[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);
> 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
More information about the R-help
mailing list