CODEWITHSUNDEEP
rfor-loopmatrix
Homer Jay Simpson
Homer Jay Simpson

Reputation: 1280

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

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 :enter image description here

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?

Upvotes: 1

Views: 64

Answers (1)

Martin Gal
Martin Gal

Reputation: 16978

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

Upvotes: 3

Related Questions