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!
Hi,
If you only want the final matrix, i.e. in this case the pm at 10
months, then you might be better off looking at something like the
square-and-multiply algorithm
(http://en.wikipedia.org/wiki/Exponentiation_by_squaring) rather than a
brute force multiplication.
Martyn
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of Jonathan Knibb
Sent: 01 February 2011 14:05
To: r-help at r-project.org
Subject: [R] better way to iterate matrix multiplication?
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.
________________________________________________________________________
This e-mail has been scanned for all viruses by Star.\ _...{{dropped:12}}
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