# [R] Limit distribution of continuous-time Markov process

Charles C. Berry cberry at tajo.ucsd.edu
Thu Jun 5 20:38:41 CEST 2008

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);
> # (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