Reputation: 341
Fill the diagonal and diagonal+1 in a matrix
I want a matrix with 0 everywhere and the diagonal and diagonal +1 has 0.5 values.
I create the matrix with the following code:
n = 10
transProbs = matrix(0, nrow = n, ncol = n)
Then, filling the diagonal with:
diag(transProbs) = 0.5
The matrix now looks the following:
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[2,] 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[3,] 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[4,] 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0
[5,] 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0
[6,] 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
[7,] 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0
[8,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0
[9,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0
[10,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
However, I want it to be:
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[2,] 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[3,] 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0
[4,] 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0
[5,] 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0
[6,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0
[7,] 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0
[8,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0
[9,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5
[10,] 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
Edit:
This matrix will be used in library(HMM), initHMM as the transProbs matrix.
My desired output for emissionProbs is:
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2
[2,] 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2
[3,] 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0
[4,] 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0
[5,] 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0
[6,] 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0
[7,] 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0
[8,] 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2
[9,] 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2
[10,] 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2
Note that it is the diag +/-2 that is filled with 0.2. In the first matrix it is the diag +1 that is filled with 0.5. This means that in the end, the probabilities might "overlap" and get in the bottom left corner.
Upvotes: 0
Views: 1655
Answers (5)
Reputation: 34773
I would use matrix indexing (allows you to replace rows according to "coordinates"); if you look at the source code for diag<- (print(`diag<-`)), you'll see it does this for the simpler diagonal-only case.
NN = nrow(transProbs)
idx = seq_len(NN)
transProbs[cbind(idx, idx)] = .5 # replace diagonal
transProbs[cbind(idx[-NN], idx[-NN] + 1L)] = .5 # replace off-diagonal
transProbs
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,] 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
# [2,] 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
# [3,] 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0
# [4,] 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0
# [5,] 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0
# [6,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0
# [7,] 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0
# [8,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0
# [9,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5
# [10,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
You could also do this in one [<- call but it's a bit uglier/harder to read:
transProbs[cbind(c(idx, idx[-NN]), c(idx, idx[-NN] + 1L))] = .5
This might be slightly more efficient but because only the first [<- call copies transProbs (see ?tracemem and .Internal(inspect(transProbs))), I thinkt he difference should be small.
Upvotes: 0
Reputation: 141
First create a function to shift a vector to the left and right (note - there is almost certainly a library or function that already does this but I couldn't find it!)
shiftSeq <- function(n, shift){
#return vector 1:n, but start shifted
# e.g. shiftSeq(5,shift=1) returns c(2,3,4,5,1)
# e.g. shiftSeq(5,shift=-1) returns c(5,1,2,3,4)
if(shift>=1){
res <- c((shift+1):n, 1:(shift))
} else if(shift==0){
res <- 1:n
} else{
res <- c((n+1+shift):n, 1:(n+shift))
}
return(res)
}
> shiftSeq(5,shift=1)
[1] 2 3 4 5 1
We'll use this shiftSeq function inside another function (below). The idea is to use apply with shiftSeq to shift each of the columns in the 'building block' diagonal matrix up and down, which we do several times, each time accumulating this shifted matrix in the result matrix.
The key is to set rowShift and colShift arguments correctly...
createTranProb <- function(n, prob, rowShift, colShift){
# create transition probability matrix of size nxn
# - prob is non-zero prob
# - rowShift is number of rows to move prob down
# - colShift is number of cols to move prob to right
shifts = setdiff(c(-rowShift:colShift), 0)
matDiag <- diag(n)*prob
matRes <- matDiag
for(i in shifts){
matRes <- matRes +
apply(matDiag, 2,
function(x) x[shiftSeq(n,i)])
}
return(matRes)
}
It works for the prob=0.5 case:
> createTranProb(10, 0.5, rowShift=0, colShift=1)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[2,] 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[3,] 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0
[4,] 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0
[5,] 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0
[6,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0
[7,] 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0
[8,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0
[9,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5
[10,] 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
As it does for prob=0.2 if we set rowShift=2 and colShift=2:
> createTranProb(10, 0.2, 2, 2)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2
[2,] 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2
[3,] 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0
[4,] 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0
[5,] 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0
[6,] 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0
[7,] 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0
[8,] 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2
[9,] 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2
[10,] 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2
Just for fun, I've added one with prob=0.33333:
> createTranProb(10, 0.33333, 1, 1)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.33333 0.33333 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.33333
[2,] 0.33333 0.33333 0.33333 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
[3,] 0.00000 0.33333 0.33333 0.33333 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
[4,] 0.00000 0.00000 0.33333 0.33333 0.33333 0.00000 0.00000 0.00000 0.00000 0.00000
[5,] 0.00000 0.00000 0.00000 0.33333 0.33333 0.33333 0.00000 0.00000 0.00000 0.00000
[6,] 0.00000 0.00000 0.00000 0.00000 0.33333 0.33333 0.33333 0.00000 0.00000 0.00000
[7,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.33333 0.33333 0.33333 0.00000 0.00000
[8,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.33333 0.33333 0.33333 0.00000
[9,] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.33333 0.33333 0.33333
[10,] 0.33333 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.33333 0.33333
Upvotes: 0
Reputation: 3492
Another solution using the fact that matrices are vectors with a dim attribute.
n <- 10
m <- 10
transProbs = matrix(0.0, nrow = n, ncol = m)
diag(transProbs) <- 0.5
transProbs[(1:(m - 1)) * (n + 1)] <- 0.5
transProbs
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> [2,] 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> [3,] 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0
#> [4,] 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0
#> [5,] 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0
#> [6,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0
#> [7,] 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0
#> [8,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0
#> [9,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5
#> [10,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
Created on 2019-09-21 by the reprex package (v0.3.0)
The sequence (1:(m - 1)) * (n + 1) selects all matrix elements one off from the diagonal.
Edit
You can achieve want you ask in the edit taking the rest of the modulus of the same sequence (plus a shift) wrt the total number of elements in the sum, i.e.
n <- 10
m <- 10
transProbs = matrix(0.0, nrow = n, ncol = m)
diag(transProbs) <- 0.2
transProbs[((1:m) * (n + 1)) %% (n * m)] <- 0.2
transProbs[((1:m) * (n + 1) + m) %% (n * m)] <- 0.2
transProbs[((1:m) * (n + 1) + 7 * m) %% (n * m)] <- 0.2
transProbs[((1:m) * (n + 1) + 8 * m) %% (n * m)] <- 0.2
transProbs
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2
#> [2,] 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2
#> [3,] 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0
#> [4,] 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0
#> [5,] 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0
#> [6,] 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0 0.0
#> [7,] 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2 0.0
#> [8,] 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2 0.2
#> [9,] 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2 0.2
#> [10,] 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.2
Created on 2019-09-21 by the reprex package (v0.3.0)
You can determine the value of the shift (i.e. +m, +7m and +8m) looking at column where the sequence starts and subtracting 2. For example, to generate the sequence which starts in the third column you have to sum (3 - 2)*m which is simply m.
I hope its clear.
Upvotes: 1
Reputation: 1676
I don't like this solution, but it does the job:
element_on_diagonal <- 0.5
element_above_and_below_diaginal <- 0.2
a <- diag(x = element_on_diagonal,
nrow = 10)
for(i in seq_len(length.out = ncol(x = a)))
{
temp <- sapply(X = setdiff(x = seq(from = (i - 2),
to = (i + 2)),
y = i),
FUN = function(j) if (j %in% 1:10) j else if (j != 0) j %% 10 else 10)
a[temp, i] <- element_above_and_below_diaginal
}
a
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 0.5 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2
#> [2,] 0.2 0.5 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2
#> [3,] 0.2 0.2 0.5 0.2 0.2 0.0 0.0 0.0 0.0 0.0
#> [4,] 0.0 0.2 0.2 0.5 0.2 0.2 0.0 0.0 0.0 0.0
#> [5,] 0.0 0.0 0.2 0.2 0.5 0.2 0.2 0.0 0.0 0.0
#> [6,] 0.0 0.0 0.0 0.2 0.2 0.5 0.2 0.2 0.0 0.0
#> [7,] 0.0 0.0 0.0 0.0 0.2 0.2 0.5 0.2 0.2 0.0
#> [8,] 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.5 0.2 0.2
#> [9,] 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.5 0.2
#> [10,] 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.2 0.2 0.5
Upvotes: 0
Reputation: 119
diag(transProbs[,-1]) = 0.5 will do it
In my terminal, the output is:
transProbs
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[2,] 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
[3,] 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0
[4,] 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0
[5,] 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0
[6,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0
[7,] 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0
[8,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0
[9,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5
[10,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
Afterwards, you can add the last "out of nowhere" 0.5 with:
transProbs[10, 1] = 0.5
Upvotes: 4
Related Questions
- Fill a matrix in R
- How to fill off-diagonals and ignore diagonals in matrix in R?
- Set values along a diagonal in a matrix
- Filling a matrix in R
- Adding a value to each element of a matrix diagonal in R
- Fill matrix in R
- Fill a matrix with sequence of value
- Looping through diagonal+1 of a matrix
- Insert number above the diagonal in R matrix
- How to fill in a matrix given diagonal and off-diagonal elements in r?