Lets say you have some Markov transition matrix, M.
We know that at steady state, there is some row vector P, such that P*M = P. We can recover that vector from the eigenvector of M' that corresponds to a unit eigenvalue.
So easy ,peasy. But suppose that M was some large symbolic matrix, with symbolic coefficients?
The problem is, if M is a transition matrix of size larger than 4, the eigenvalue/eigenvector computation becomes equivalent to solving for the roots of a polynomial of degree higher than 4. We provably cannot solve a general problem like that for an analytical solution.
But we don't need to solve for the eigenvalues, because if this is a Markov transition matrix, then all we need is to solve for a specific eigenvector. First, a simple matrix with all numerical values.
M = [.1 .2 0 .7 0;.3 .4 .3 0 0;0 0 .5 0 .5;.2 .2 .2 .2 .2;0 .1 .2 .3 .4];
You can verify this is a valid transition matrix. As expected, that eigenvalue of 1 is there.
eig(M')
ans =
1
-0.21443
0.49014
0.13762
0.18667
[V,D] = eig(M');
P = V(:,1)';P = P./sum(P)
P =
0.0898 0.14276 0.28045 0.18996 0.29703
P*M
ans =
0.0898 0.14276 0.28045 0.18996 0.29703
So indeed, we have extracted the steady state transition probabilities in the vector P.
But what if M is not a simple numerical matrix?
M = [.1 .2 0 .7 0;.3 .4 .3 0 0;0 0 k 0 1-k;.2 .2 .2 .2 .2;0 .1 .2 .3 .4]
M =
[ 1/10, 1/5, 0, 7/10, 0]
[ 3/10, 2/5, 3/10, 0, 0]
[ 0, 0, k, 0, 1 - k]
[ 1/5, 1/5, 1/5, 1/5, 1/5]
[ 0, 1/10, 1/5, 3/10, 2/5]
As it turns out, the symbolic toolbox can actually solve for the eigenvalues of M here. But we don't need to go as far as eig. As long as we know that M is a valid transition matrix, then we need only solve the linear system:
P = P*M
subject to the constraint that sum(P)==1. We need the constraint that sum(P)==1 because the matrix problem we have formulated is singular. That is, an all zero vector P will satisfy the above problem. So we need to supply another piece of information. That sum(P)==1 is a good way. There are several ways we might accomplish the solution. We could use solve from the symbolic toolbox. Or we could just use slash.
In the case of the original matrix M, use of slash will reduce to this:
M = [.1 .2 0 .7 0;.3 .4 .3 0 0;0 0 k 0 1-k;.2 .2 .2 .2 .2;0 .1 .2 .3 .4]
M =
[ 1/10, 1/5, 0, 7/10, 0]
[ 3/10, 2/5, 3/10, 0, 0]
[ 0, 0, k, 0, 1 - k]
[ 1/5, 1/5, 1/5, 1/5, 1/5]
[ 0, 1/10, 1/5, 3/10, 2/5]
n = size(M,1);
P = [zeros(1,n),1]/[M-eye(n),ones(n,1)]
P =
[ (390*(k - 1))/(3125*k - 3734), (620*(k - 1))/(3125*k - 3734), -609/(3125*k - 3734), (825*(k - 1))/(3125*k - 3734), (1290*(k - 1))/(3125*k - 3734)]
When k=0.5, we get the vector we had before.
double(subs(P,k,.5))
ans =
0.0898 0.14276 0.28045 0.18996 0.29703
At the extremes of k, we see these limiting cases:
subs(P,k,0)
ans =
[ 195/1867, 310/1867, 609/3734, 825/3734, 645/1867]
subs(P,k,1)
ans =
[ 0, 0, 1, 0, 0]
A problem becomes when the transition matrix gets larger and more general. Then symbolic solution as I did above may become intractable. Suppose I gave you some almost completely general 20x20 matrix M, the symbolic solve will become so complicated that you will run out of memory. In fact, I'd bet that might happen even for a completely general 10x10 matrix. Solution of symbolic linear systems get very big VERY fast.
As far as the problems you show, just create the corresponding transition matrix, M. Then it becomes virtually one line of code.
