1

I have this matrix:

            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

And this vector:

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

And I am trying to come up with a loop that will keep multiplicating the results (new vectors) by the same matrix (it is a simple exemple of population matrices model). I need all the results from the 1st multiplication, until 100th, so I can put them into a graph. Any ideas?

> 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."))
Matthew Lundberg
  • 42,009
  • 6
  • 90
  • 112
TWest
  • 765
  • 2
  • 6
  • 27

2 Answers2

5

This should return a list of the vectors:

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

Note that the first element of the list is the original vector, and there are 100 subsequent elements.

Matthew Lundberg
  • 42,009
  • 6
  • 90
  • 112
3

A bit of maths will help here, you want to calculate

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

Calculating the power of a matrix is standard matrix algebra (eigenvalues / diagonalization of a matrix)

The expm package has the %^% operator to do this

 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)
mnel
  • 113,303
  • 27
  • 265
  • 254