[R] Limit distribution of continuous-time Markov process

[Charles C. Berry]

On Thu, 5 Jun 2008, gschultz at scriptpro.com wrote:

>
> I have (below) an attempt at an R script to find the limit distribution
> of
> 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?

The exposition there seemed pretty clear.

Maybe you need to brush up on algebra - particularly matrix 
decompositions. I'd go to Rao's book (Linear Statistical Inference and Its 
Applications,1973) and look at the early chapter (sorry I can't say which 
one - I am at home now and do not have a copy at hand) that covers matrix 
decomposition/diagonalization and work a few exercises to get up to speed.

But there are lots of texts that would cover this stuff, so pick one and 
have at it.

The point of that exposition is that you can get the matrix exponential of 
the transition matrix from the matrix exponential of the diagonal matrix 
of eigenvalues and the eigenvectors of the rate matrix.

So, eigen() does most of the work.

Once you realize what $\lim_{ t \to \infty }{ E^{tD }$ is (using the 
notation of the link you provided), you will see the computation is 
trivial given the results of eigen().

HTH,

Chuck



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
> form
> 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);
> Answer;
> # (Right answer for this example is a matrix with all values = 1/3.)
>
>
> Grant D. Schultz
>
> ______________________________________________
> 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.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901