[R] Limit distribution of continuous-time Markov process

gschultz at scriptpro.com gschultz at scriptpro.com
Thu Jun 5 15:41:37 CEST 2008

I have (below) an attempt at an R script to find the limit distribution
a continuous-time Markov process, using the formulae outlined at 
http://www.uwm.edu/~ziyu/ctc.pdf, page 5.

First, is there a better exposition of a practical algorithm for doing 
this? I have not found an R package that does this specifically, nor 
anything on the web.

Second, the script below will give the right answer, _if_ I "normalize" 
the rate matrix, so that the average rate is near 1.0, and _if_ I 
tweak the multiplier below (see **), and then watch for the Answer to 
converge to a matrix where the rows to sum to 1.0. (This multiplier 
is "t" in the PDF whose URL is above.) Is there a known way to get 
this to converge?

Thank you.

---------------R script:--------------
# The rate matrix:
Q <- matrix(c(-1, 1, 0, 1, -2, 1, 0, 1, -1), ncol=3, byrow=T);
M <- eigen(Q)$vectors # diagonalizer matrix
D <- ginv(eigen(Q)$vectors) %*% Q %*% eigen(Q)$vectors # Diagonalized 
Sum <- matrix(c(rep(0, 9)), ncol=3, byrow=T);
for (i in 0:60)
{ # Naive, Taylor series summation:
	Temp <- D;
	diag(Temp) <- (4 * diag(D)) ^ i; # ** =4
	Sum <- Sum + Temp / factorial(i);

Answer <- M %*% Sum %*% ginv(M);
# (Right answer for this example is a matrix with all values = 1/3.)

Grant D. Schultz

More information about the R-help mailing list