Jesper.Lindberg
Jesper.Lindberg

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)

MichaelChirico
MichaelChirico

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

Ben Wynne-Morris
Ben Wynne-Morris

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

agila
agila

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

yarnabrina
yarnabrina

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

GabCaz
GabCaz

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