How can I fix this wrong result in matrix for loop in R?-CodePudding

I want to enter in R the following matrix (transition probability) but something I make wrong and the output is not right. The matrix is in the picture :

my effort is the following :

n=100
A = matrix(c(0),n 1,n 1);A
A[1,2] = 1/4 
A[1,1] = 3/4
A[2,1] = 1/4 
A[2,2] = 1/2 
A[2,1] = 1/4
for (k in 2:n) {
  A[k,k 1] = 1/4
  A[k,k]   = 1/2 
  A[k,k-1] = 1/6 
  A[k,k-2] = 1/12
  A[n,n]   = 3/4 
  A[n,n-1] = 1/6 
  A[n,n-2] = 1/12}
A


pA = A
for (i in 10000){
  pA = pA %*% A }
pA

but the resulting columns (first 3 columns) must be: 0.2087, 0.1652, 0.1307

Any help?

CodePudding user response：

You could try a slightly different approach here:

n <- 100
A <- matrix(0, n   1, n   1)

# create index for the diagonal and the "other" diagonales
A_diag <- matrix(c(1:(n   1), 1:(n   1)), ncol = 2)
A_diag_1 <- matrix(c(1:(n   1), 1   1:(n   1)), ncol = 2)
A_diag_minus_1 <- matrix(c(1:(n   1), - 1   1:(n   1)), ncol = 2)
A_diag_minus_2 <- matrix(c(1:(n   1), - 2   1:(n   1)), ncol = 2)

# populate those diagonales
A[A_diag] <- 1/2
A[A_diag_1[A_diag_1[,2] <= n   1,]] <- 1/4
A[A_diag_minus_1[A_diag_minus_1[,2] >= 1,]] <- 1/6
A[A_diag_minus_2[A_diag_minus_2[,2] >= 1,]] <- 1/12

# create starting values
A[1, 1] <- 3/4
A[2, 1] <- 1/4

# calculate pA
pA <- A
for (i in 1:10000){ pA <- pA %*% A }

This returns

pA[1, 1:3]
#> [1] 0.2087122 0.1651514 0.1306823

For the calulation of pA you could use the expm package to speed things up:

library(expm)

pB <- A %^% 10000

pB[1,1:3]
#> [1] 0.2087122 0.1651514 0.1306823