[R] better way to iterate matrix multiplication?

[Eik Vettorazzi]

yes there are. But Christofer doesn't need exp(A), but A^n.
But there is a matpow-function %^% in this package, which is a little
bit slower, I think:

library(expm)
states<-1000
tm <- matrix(runif(states^2),nrow=states)  # random transition matrix
for illustration
tm <- t(apply(tm,1,function (x) x/sum(x)))  # make its rows sum to 1
p0<-pm <- c(0.5,0.5,rep(0,states-2))  # half of cases start in state 1,
half in state 2
n<-10000

system.time({dd<-eigen(tm,symmetric=F)
as.real(p0%*% dd$vectors%*% diag(dd$values^n)%*%solve(dd$vectors))})
       User      System     elapsed
      15.20        0.09       15.57

system.time(p0%*%(tm%^%n))
       User      System     elapsed
      38.61        0.00       39.62



Am 01.02.2011 17:16, schrieb Ben Bolker:
> Eik Vettorazzi <E.Vettorazzi <at> uke.uni-hamburg.de> writes:
> 
>>
>> 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.
> 
>   There are matrix exponential functions in the Matrix and expm
> packages ... don't know about their speed
> 
> ______________________________________________
> 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