From: <gschultz_at_scriptpro.com>

Date: Thu, 05 Jun 2008 08:41:37 -0500

Sum <- matrix(c(rep(0, 9)), ncol=3, byrow=T); for (i in 0:60)

{ # Naive, Taylor series summation:

}

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 - 15:46:12 GMT

Date: Thu, 05 Jun 2008 08:41:37 -0500

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? 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 # Diagonalizedform

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. Received on Thu 05 Jun 2008 - 15:46:12 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 - 21:31:04 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.
*