Re: [R] Limit distribution of continuous-time Markov process

From: Charles C. Berry <cberry_at_tajo.ucsd.edu>
Date: Thu, 05 Jun 2008 11:38:41 -0700

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

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. Received on Thu 05 Jun 2008 - 20:13:48 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 05 Jun 2008 - 20:32:34 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive