look at function ?mtx.exp() in the Malmig package, e.g., A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE) n.zero <- c(1,0,0) ################## library(Malmig)
all.equal(A %*% (A %*% n.zero), mtx.exp(A, 2) %*% n.zero) all.equal(A %*% (A %*% (A %*% n.zero)), mtx.exp(A, 3) %*% n.zero) mtx.exp(A, 15) %*% n.zero I hope it helps. Best, Dimitris ---- Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/16/336899 Fax: +32/16/337015 Web: http://www.med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm ----- Original Message ----- From: "Rau, Roland" <[EMAIL PROTECTED]> To: "R-Help" <[email protected]> Sent: Wednesday, August 17, 2005 12:15 PM Subject: [R] power of a matrix > Dear all, > > I have a population with three age-classes, at time t=0 the > population > is: > n.zero <- c(1,0,0) > I have a transition matrix A which denotes "fertility" and > "survival": > > A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE) > > To obtain the population at t=1, I calculate: > A %*% n.zero > > To obtain the population t=2, I calculate: > A %*% (A %*% n.zero) > > ... and so on ... > > I thought now to obtain the population at time x, I should simply > do: > A^x %*% n.zero > > But this, of course, does not work since the following two > statements > are not equivalent for matrices: > A %*% A > A * A > > Is there something like a "powermatrix"-function? > > Thanks, > Roland > > P.S. The example is taken from: > Example 3.1 ("A linear, time-invariant model") from Keyfitz/Caswell: > "Applied Mathematical Demography", 3rd Edition, 2005, p. 52f > > > +++++ > This mail has been sent through the MPI for Demographic > Rese...{{dropped}} > > ______________________________________________ > [email protected] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
