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-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Jonathan Knibb Sent: 01 February 2011 14:05 To: r-help@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@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}} ______________________________________________ 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.
