Solved – Transition probability matrix rows not summing to 1

by

I'm implementing a method of sampling from a CTMC from here.

I'm trying to calculate the transition probability matrix but the rows are not summing to 1, except at 0.

I've diagonalised my rate matrix $Q$ to give me a matrix $D$ such that $Q=UDU^{-1}$ where $U$ is the matrix with it's columns the eigenvectors of $Q$ and $D$ the diagonal matrix with diagonal entries the eigenvalues of $Q$.

I should then be able to calculate the transition probability matrix of the CTMC using the equation $$P(t)=e^{Qt}=Ue^{tQ}U^{-1}.$$

I've implemented this using R, with the following code:

eig <- eigen(rateM) U <- eig$vectors invU <- solve(U)  Pt <- U%*%diag(exp(eig$values*t))%*%invU

where

> rateM            [,1]       [,2]       [,3]       [,4] [1,] -1.1224490  0.6122449  0.3061224  0.2040816 [2,]  0.4081633 -0.9183673  0.3061224  0.2040816 [3,]  0.2040816  0.3061224 -0.9183673  0.2040816 [4,]  0.2040816  0.3061224  0.6122449 -1.1224490 > t [1] 3

The rows of the rate matrix all sum to zero, $UU^{-1}=I$ as it should, but the larger t gets, the further from 1 Pt gets. E.g.

> sum((U%*%diag(exp(eig$values*0))%*%invU)[1,])     [1] 1     > sum((U%*%diag(exp(eig$values*1))%*%invU)[1,]) [1] 0.9773404 > sum((U%*%diag(exp(eig$values*2))%*%invU)[1,])     [1] 0.9316814     > sum((U%*%diag(exp(eig$values*3))%*%invU)[1,]) [1] 0.8799868 > sum((U%*%diag(exp(eig$values*10))%*%invU)[1,]) [1] 0.5701228

Is this owing to some sort of numerical error, or am I doing something incorrectly?

Contents hide
1 Best Answer
2 Similar Posts:

Best Answer

Converting comment to answer, since the comment was on the money, thereby allowing this to be closed. Credit to whuber for his initial comment.

The 3rd row of rateM sums to -0.2040817, not to zero, and is therefore erroneous. Therefore the transition probability matrix will be erroneous.

Similar Posts:

Rate this post

Leave a Comment