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@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
______________________________________________
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.
