矩阵乘法循环

Question

我有这个矩阵：

            V1         V2         V3        V4         V5         V6         V7 V8
[1,] 0.8399983 0.01558029 0.00000000 0.0000000 0.00000000 0.00000000 0.00000000  0
[2,] 0.0000000 0.89022017 0.02570281 0.0000000 0.00000000 0.00000000 0.00000000  0
[3,] 0.0000000 0.00000000 0.87910624 0.0242963 0.00000000 0.00000000 0.00000000  0
[4,] 0.0000000 0.00000000 0.00000000 0.0000000 0.03428571 0.00000000 0.00000000  0
[5,] 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 0.02988506 0.00000000  0
[6,] 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 0.73438228 0.01666667  0
[7,] 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.00000000  0
[8,] 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.00000000  0

还有这个向量：

     [,1]
[1,]  908
[2,]  516
[3,]  269
[4,]   85
[5,]   32
[6,]   13
[7,]    2
[8,]    3

我正在尝试提出一个循环，该循环将继续将结果（新向量）乘以相同的矩阵（这是人口矩阵模型的一个简单示例）。我需要第一次乘法的所有结果，直到第 100 次，所以我可以将它们放入图表中。有什么想法吗？

> dput(mat)
structure(c(0.8399983, 0, 0, 0, 0, 0, 0, 0, 0.01558029, 0.89022017, 
0, 0, 0, 0, 0, 0, 0, 0.02570281, 0.87910624, 0, 0, 0, 0, 0, 0, 
0, 0.0242963, 0, 0, 0, 0, 0, 0, 0, 0, 0.03428571, 0, 0, 0, 0, 
0, 0, 0, 0, 0.02988506, 0.73438228, 0, 0, 0, 0, 0, 0, 0, 0.01666667, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(8L, 8L), .Dimnames = list(
    c("[1,]", "[2,]", "[3,]", "[4,]", "[5,]", "[6,]", "[7,]", 
    "[8,]"), c("V1", "V2", "V3", "V4", "V5", "V6", "V7", "V8"
    )))
> dput(vec)
structure(c(908L, 516L, 269L, 85L, 32L, 13L, 2L, 3L), .Dim = c(8L, 
1L), .Dimnames = list(c("[1,]", "[2,]", "[3,]", "[4,]", "[5,]", 
"[6,]", "[7,]", "[8,]"), "X..1."))

Answer 1

这应该返回一个向量列表：

Reduce(f=function(v, x) mat %*% v, x=1:100, init=vec, accumulate=TRUE)

注意列表的第一个元素是原始向量，后面有100个元素。

Answer 2

这里有一点数学帮助，你想计算一下

p_1 =  M \times p
p_2 =  M  \times p_1 = M \times M  \times p = M^2 \times p
....
p_n = M^n \times p

计算矩阵的幂是标准矩阵代数（矩阵的特征值/对角化）

expm 包有 %^% 运算符来执行此操作

 library(expm)

 # a function to calculate 
 # the population at time n for transition matrix M and initial population p


pop <- function(n, M, p){
     (M%^%n) %*% p
   }

# use sapply to get the populations (each column is a time point

sapply(1:100, pop, M = M, p = p)