Main Content

Markov Chain Analysis and Stationary Distribution

This example shows how to derive the symbolic stationary distribution of a trivial Markov chain by computing its 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.

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 probabilities. This will be helpful in selecting desirable stationary distributions later.

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

Define the 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 = 

(a1-a00001-bb0000cCAcCBc1-cCA-cCB-c00000d000000e1-e00001-ff)[a, 1 - a, sym(0), sym(0), sym(0), sym(0); 1 - b, b, sym(0), sym(0), sym(0), sym(0); cCA, cCB, c, 1 - cCA - cCB - c, sym(0), sym(0); sym(0), sym(0), sym(0), d, sym(0), sym(0); sym(0), sym(0), sym(0), sym(0), e, 1 - e; sym(0), sym(0), sym(0), sym(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

V
V = 

(b-1a-10-c-dcCB-bcCA-bcCB+ccCAσ10-1010-c-dcCA-acCA-acCB+ccCBσ101000-c-dc+cCA+cCB-10000011000f-1e-1000-1010001)where  σ1=c+cCA+cCB-1a+b-ac-bc+c2-1[(b - 1)/(a - 1), sym(0), -((c - d)*(cCB - b*cCA - b*cCB + c*cCA))/((c + cCA + cCB - 1)*(a + b - a*c - b*c + c^2 - 1)), sym(0), -sym(1), sym(0); sym(1), sym(0), -((c - d)*(cCA - a*cCA - a*cCB + c*cCB))/((c + cCA + cCB - 1)*(a + b - a*c - b*c + c^2 - 1)), sym(0), sym(1), sym(0); sym(0), sym(0), -(c - d)/(c + cCA + cCB - 1), sym(0), sym(0), sym(0); sym(0), sym(0), sym(1), sym(1), sym(0), sym(0); sym(0), (f - 1)/(e - 1), sym(0), sym(0), sym(0), -sym(1); sym(0), sym(1), sym(0), sym(0), sym(0), sym(1)]

Analytical eigenvalues

diag(D)
ans = 

(11cda+b-1e+f-1)[sym(1); sym(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 = 

(11d)[sym(1); sym(1); d]

Extract the analytical stationary distributions. The eigenvectors are normalized with the 1-norm or sum(abs(X)) prior to display.

for k = ix'
    V(:,k) = simplify(V(:,k)/norm(V(:,k)),1);
end
Probability = V(:,ix)
Probability = 

(b-1a-1σ2001σ2000000010f-1σ1e-1001σ10)where  σ1=f-12e-12+1  σ2=b-12a-12+1[(b - 1)/((a - 1)*sqrt((b - 1)^2/(a - 1)^2 + 1)), sym(0), sym(0); 1/sqrt((b - 1)^2/(a - 1)^2 + 1), sym(0), sym(0); sym(0), sym(0), sym(0); sym(0), sym(0), sym(1); sym(0), (f - 1)/(sqrt((f - 1)^2/(e - 1)^2 + 1)*(e - 1)), sym(0); sym(0), 1/sqrt((f - 1)^2/(e - 1)^2 + 1), sym(0)]

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.

fsurf(Probability(1), [0 1 0 1]);
xlabel a
ylabel b
title('Probability of A');

figure(2);
fsurf(Probability(2), [0 1 0 1]);
xlabel a
ylabel b
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 third row is entirely zero.

  • The rest of the states form three groups, { A , B }, { D } and { E , F } that do not communicate with each other and are recurrent.