[R] better way to iterate matrix multiplication?

Martyn Byng Martyn.Byng at nag.co.uk
Tue Feb 1 16:22:25 CET 2011


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}}



More information about the R-help mailing list