On Thu, 5 Jun 2008, [EMAIL PROTECTED] 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
______________________________________________
[email protected] 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:[EMAIL PROTECTED] UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
______________________________________________
[email protected] 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.
