Symbolic Math Toolbox

Symbolic Analysis of a Markov Chain

In this example, we derive the symbolic stationary distribution of a trivial Markov chain by computing it's eigen decomposition. The stationary distribution represents the limiting, time-independent, distribution of the states for a Markov process as the number of steps or transitions increase.

This example uses the following Symbolic Math Toolbox capablities:

  • Computing eigenvalues using eig

Define (positive) transition probabilities between states { A through F } as shown in the above image.

syms a b c d e f cCA cCB positive;

Add further assumptions bounding the transition probablities. This will be helpful in selecting desirable stationary distributions later.

assumeAlso([a, b, c, e, f, cCA, cCB] < 1 & d == 1);

Define transition matrix. States { A through F } are mapped to the columns and rows { 1 through 6 }. Note the values in each row sum up to one.

P = sym(zeros(6,6));
P(1,1:2) = [a 1-a];
P(2,1:2) = [1-b b];
P(3,1:4) = [cCA cCB c (1-cCA-cCB-c)];
P(4,4) = d;
P(5,5:6) = [e 1-e];
P(6,5:6) = [1-f f];
P
 
P =
 
[     a, 1 - a, 0,                 0,     0,     0]
[ 1 - b,     b, 0,                 0,     0,     0]
[   cCA,   cCB, c, 1 - cCA - cCB - c,     0,     0]
[     0,     0, 0,                 d,     0,     0]
[     0,     0, 0,                 0,     e, 1 - e]
[     0,     0, 0,                 0, 1 - f,     f]
 

Compute all possible analytical stationary distributions of the states of the Markov chain. This is the problem of extracting eigenvectors with corresponding eigenvalues that can be equal to 1 for some value of the transition probabilities.

[V,D] = eig(P');

Analytical eigenvectors

pretty(V)
/ b - 1           (c - d) (cCB - b cCA - b cCB + c cCA)            \
| -----,   0,   - -------------------------------------, 0, -1,  0 |
| a - 1                             #1                             |
|                                                                  |
|                 (c - d) (cCA - a cCA - a cCB + c cCB)            |
|   1,     0,   - -------------------------------------, 0,  1,  0 |
|                                   #1                             |
|                                                                  |
|                                 c - d                            |
|   0,     0,             - -----------------,           0,  0,  0 |
|                           c + cCA + cCB - 1                      |
|                                                                  |
|   0,     0,                      1,                    1,  0,  0 |
|                                                                  |
|        f - 1                                                     |
|   0,   -----,                    0,                    0,  0, -1 |
|        e - 1                                                     |
|                                                                  |
\   0,     1,                      0,                    0,  0,  1 /

where

                                                   2
   #1 == (c + cCA + cCB - 1) (a + b - a c - b c + c  - 1)


Analytical eigenvalues

diag(D)
 
ans =
 
         1
         1
         c
         d
 a + b - 1
 e + f - 1
 

Find eigenvalues that are exactly equal to 1. If there is any ambiguity in determining this condition for any eigenvalue, stop with an error - this way we are sure that below list of indices is reliable when this step is successful.

ix = find(isAlways(diag(D) == 1,'Unknown','error'));
diag(D(ix,ix))
 
ans =
 
 1
 1
 d
 

Extract the analytical stationary distributions. The eigenvectors are normalized prior to display

for k = ix'
    V(:,k) = simplify(V(:,k)/norm(V(:,k)));
end
Probability = V(:,ix);
pretty(Probability)
/    b - 1                  \
| ----------,      0,     0 |
| (a - 1) #2                |
|                           |
|      1                    |
|     --,          0,     0 |
|     #2                    |
|                           |
|      0,          0,     0 |
|                           |
|      0,          0,     1 |
|                           |
|                f - 1      |
|      0,     ----------, 0 |
|             #1 (e - 1)    |
|                           |
|                  1        |
|      0,         --,     0 |
\                 #1        /

where

             /        2     \
             | (f - 1)      |
   #1 == sqrt| -------- + 1 |
             |        2     |
             \ (e - 1)      /

             /        2     \
             | (b - 1)      |
   #2 == sqrt| -------- + 1 |
             |        2     |
             \ (a - 1)      /


The probability of the steady state being A or B in the first eigenvector case is a function of the transition probabilities a and b. Visualize this dependency.

ezsurf(Probability(1), [0 1], [0 1]);
title('Probability of A');
figure(2);
ezsurf(Probability(2), [0 1], [0 1]);
title('Probability of B');

The stationary distributions confirm the following (Recall states { A through F } correspond to row indices { 1 through 6 }):

  • State C is never reached and is therefore transient i.e. the 3rd row is entirely zero.
  • Rest of the states form three groups, { A , B }, { D } and { E , F } that do not communicate with each other and are recurrent.