Continuous-Time Markov Chain populations

[jinx]

Hello,

I'm working on a CTMC three-state model to obtain time-dependent populations of each state.

S <=> E <=> G
I have built a rate matrix for this (diffusion) process.

<br /> K =<br /> \begin{pmatrix}<br /> K_{SS} &amp; K_{SE} &amp; K_{SG}\\<br /> K_{ES} &amp; K_{EE} &amp; K_{EG}\\<br /> K_{GS} &amp; K_{GE} &amp; K_{GG}<br /> \end{pmatrix}<br /> =<br /> \begin{pmatrix}<br /> -3.13E+06 &amp; 4.29E+07 &amp; 0\\<br /> 3.13E+06 &amp;-4.29E+07 &amp; 3.33E+09\\<br /> 0 &amp; 2.26E+06 &amp; -3.33E+09<br /> \end{pmatrix}<br />

The time-dependent population of state G (final state), is given as the product of (I) sum of the (elements??) of the matrix exponential of the transition matrix multiplied by time, and (II) the initial population of state S, P_S(0)=1.

<br /> P_G(t) = \left( \sum_{S}(exp(tK))\right)*P_S(0)=\left(\sum_{S}(\frac{t^n*K^n}{n!} )\right)*P_S(0)<br />

My question: I have many difficulties understanding how I could solve this matrix exponential to obtain the population of state G in practice, given the above transition matrix.In this handout, they estimate the population but for the initial state (here, denoted S) instead of final (G, what I seek). They do so by obtaining by an eigenvalue decomposition of the rate matrix, to obtain a final expression for P(t) in terms of t.
http://www.stats.ox.ac.uk/~laws/AppliedProb/handout5.pdf

 

[Stephen Tashi]

Your question is how to compute the matrix exponential \exp(tK). If this is a numerical computation, what programming language will you use?

 

[jinx]

Let me rephrase my question with more detail:

I'm working on a CTMC three-state model to obtain time-dependent populations of each state.
S &lt;=&gt; E &lt;=&gt; G
I have built a rate matrix for this (diffusion) process.

<br /> K =<br /> \begin{pmatrix}<br /> K_{SS} &amp; K_{SE} &amp; K_{SG}\\<br /> K_{ES} &amp; K_{EE} &amp; K_{EG}\\<br /> K_{GS} &amp; K_{GE} &amp; K_{GG}<br /> \end{pmatrix}<br /> =<br /> \begin{pmatrix}<br /> -3.13E+06 &amp; 4.29E+07 &amp; 0\\<br /> 3.13E+06 &amp;-4.29E+07 &amp; 1.90E+10\\<br /> 0 &amp; 1.83E+06 &amp; -1.90E+09<br /> \end{pmatrix}<br />

The time-dependent probability population of state G (final state), is given as the product of (I) sum of the (elements??) of the matrix exponential of the transition matrix multiplied by time, and (II) the initial population of state S, P_S(0)=1.

<br /> P_G(t) = \left( \sum_{S}(exp(tK))\right)*P_S(0)=\left(\sum_{S}(\frac{t^n*K^n}{n!} )\right)*P_S(0)<br />

My question: I'm trying to work out the probability (population) of state G (state C). Initially only state S is population (P_S(0) = 1, P_E(0)=0 and P_G(0)=0 )

(1) In order to obey Markov probability conservation, I have ensured that matrix K obeys the property that
each column adds to zero. Is this correct?(2) One consequence of obeying probability conservation is that upon diagonalisation, one eigenvalue
will be zero. What is the consequence of a zero eigenvalue for the physics of the problem?<br /> K =<br /> \begin{pmatrix}<br /> \lambda_{1} &amp; 0 &amp; 0\\<br /> 0 &amp; \lambda_{2} &amp; 0\\<br /> 0 &amp; 0 &amp; \lambda_{3}<br /> \end{pmatrix}<br /> =<br /> \begin{pmatrix}<br /> 0 &amp; 0 &amp; 0\\<br /> 0 &amp; -46247656 &amp; 0\\<br /> 0 &amp; 0 &amp; -19022904594<br /> \end{pmatrix}<br />

When I solve for the eigenvalues in R, it lists the eigenvalues according to INCREMENTAL size, so -19022904594, -46247656, 0
however from linear algebra exercises the value is 0, -46247656, -19022904594 so how to determine the right order of eigenvalues?

A textbook example is<br /> \begin{pmatrix}<br /> -2 &amp; 1 &amp; 1\\<br /> 1 &amp;-1 &amp; 0\\<br /> 2 &amp; 1 &amp; -3<br /> \end{pmatrix}<br />

This gives eigenvalues 0, -2, -4 [true order] but R program ranks them as -4, -2, 0.
Is there a way to determine the true order of eigenvalues?

(3) When I input this, I build a probability equation using three arbitrary diagonalisation
constants
\alpha, \beta, \gamma<br />

P_(G) = \alpha*(exp(- \lambda_{1} *t))+\beta*(exp(- \lambda_{2} *t))+\gamma*(exp(- \lambda_{3} *t))

The three constants are solved using a 3x3 system of equations depending on the populations
of state G ( as per matrix exponential algebra of Markov chain models using the forward and
backward equations)
<br /> \alpha+\beta+\gamma=1<br />
<br /> \alpha*\lambda_{1} + \beta*\lambda_{2} + \gamma*\lambda_{3} = k_GG<br />
<br /> \alpha*( \lambda_{1})^2 + \beta*(\lambda_{2})^2 + \gamma*(\lambda_{3})^2 = k_EG*k_GG + k_GG^2<br />

Now, depending on whether we choose the constraint
<br /> \alpha + \beta + \gamma =1<br />

<br /> \alpha + \beta + \gamma =0<br />

I get very different behaviour. Do we solve this relative to a state with initial population 1 or 0?I attach a link for this

http://www.stats.ox.ac.uk/~laws/AppliedProb/handout5.pdf