Re: Limit distrib. of contin.-time Markov process (in R)

On May 22, 9:53 am, gschu...@xxxxxxxxxxxxx wrote:
I am attempting to write an R script to find the limit distribution of
a continuous-time Markov process, using the formulae outlined athttp://www..uwm.edu/~ziyu/ctc.pdf, page 5.

1.  Is there a better exposition of a practical algorithm for doing
this?  I have not found an R package that does this specifically, nor
anything on the web.
2.  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 very much.

---------------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 is the matrix with all values = 1/3.

You probably should consider moving this question
to the R-help mailing list. See the link "Mailing lists" on the R
home.
page.
.