[R] better way to iterate matrix multiplication?

[Eik Vettorazzi]

if you have a homogeneous mc (= a constant transition matrix), your
state at time 10 is given by (chapman-kolmogorov)
p10=p0 %*% tm^(10)
so you need a matrix power function.

You can use the eigendecomposition and some linear algebra

A^n=(VDV^{-1})^n)=VD^nV^{-1}

dd<-eigen(tm,symmetric=F)
as.real(p0%*% dd$vectors%*% diag(dd$values^n)%*%solve(dd$vectors))

but for your toy example, there is no difference in time consumed.

hth

Am 01.02.2011 15:04, schrieb Jonathan Knibb:
> I'm simulating a Markov process using a vector of proportions. Each
> value in the vector represents the proportion of the population who are
> in a particular state (so the vector sums to 1). I have a square matrix
> of transition probabilities, and for each tick of the Markov clock the
> vector is multiplied by the transition matrix.
> 
> To illustrate the sort of thing I mean:
> 
> pm <- c(0.5,0.5,0)  # half of cases start in state 1, half in state 2
> tm <- matrix(runif(9),nrow=3)  # random transition matrix for illustration
> tm <- t(apply(tm,1,function (x) x/sum(x)))  # make its rows sum to 1
> total.months = 10
> for(month in 1:total.months) {pm <- pm %*% tm}  # slow!
> pm  # now contains the proportion of cases in each state after 10 months
> 
> My question is: is there a quicker, more R-idiomatic way of doing it,
> avoiding 'for'? I've been trying to use apply() to fill a matrix with
> the vectors, but I can't get this to act iteratively.
> 
> Suggestions gratefully received!
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


-- 
Eik Vettorazzi
Institut für Medizinische Biometrie und Epidemiologie
Universitätsklinikum Hamburg-Eppendorf

Martinistr. 52
20246 Hamburg

T ++49/40/7410-58243
F ++49/40/7410-57790