>>>>> "EV" == Eik Vettorazzi <e.vettora...@uke.uni-hamburg.de>
>>>>> on Tue, 1 Feb 2011 18:02:42 +0100 writes:
EV> yes there are. But Christofer doesn't need exp(A), but
EV> A^n.
EV> But there is a matpow-function %^% in this package, which is a little
EV> bit slower, I think:
EV> library(expm)
EV> states<-1000
EV> tm <- matrix(runif(states^2),nrow=states) # random transition matrix
EV> for illustration
EV> tm <- t(apply(tm,1,function (x) x/sum(x))) # make its rows sum to 1
EV> p0<-pm <- c(0.5,0.5,rep(0,states-2)) # half of cases start in state 1,
EV> half in state 2
EV> n<-10000
EV> system.time({dd<-eigen(tm,symmetric=F)
EV> as.real(p0%*% dd$vectors%*% diag(dd$values^n)%*%solve(dd$vectors))})
EV> User System elapsed
EV> 15.20 0.09 15.57
EV> system.time(p0%*%(tm%^%n))
EV> User System elapsed
EV> 38.61 0.00 39.62
Indeed, thank you for doing the experiment.
Note however that the eigen() method is only available in *some* cases for
the matrix power, e.g. only when the matrix is non-singular,
whereas the `%^%` in package expm should work reliably in all
cases.
Also for (much) smaller powers than n=10'000
the cpu time needed is more comparable.
Martin Maechler, ETH Zurich
EV> 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
EV> --
EV> Eik Vettorazzi
EV> Institut für Medizinische Biometrie und Epidemiologie
EV> Universitätsklinikum Hamburg-Eppendorf
EV> Martinistr. 52
EV> 20246 Hamburg
EV> T ++49/40/7410-58243
EV> F ++49/40/7410-57790
______________________________________________
R-help@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.
