如何在 R 中的矩阵 for 循环中修复此错误结果？

Question

我想在 R 中输入以下矩阵（转移概率），但我做错了，输出不正确。 矩阵在图中：

我的努力如下：

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

但结果列（前 3 列）必须为：0.2087、0.1652、0.1307

有什么帮助吗？

Answer 1

您可以在这里尝试稍微不同的方法：

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 }

返回

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

对于pA 的计算，您可以使用expm 包来加快速度：

library(expm)

pB <- A %^% 10000

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