Петли для матриц умножения

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 

И я пытаюсь придумать цикл, который будет держать перемножение результатов (новые векторы) по той же матрице (это простой пример модели матриц совокупности). Мне нужны все результаты от 1-го умножения, до 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)