Markov Chains

Matrix Powers

The key to understanding the properties of Markov chains is to understand the powers P^n of the transition matrix P.  Here is the Sunny / Rainy weather chain.

P = {{7/10,3/10},{8/10,2/10}};
{P.P//MatrixForm,P.P.P//MatrixForm,MatrixPower[P,10]//N//MatrixForm}

P^n  converges to a limit matrix which gives the long term probabilities of being Sunny or Rainy.

Here it the ant-on-a-grid chain.  Once again, it converges to a limit.

P = {{0,1/3,1/3,1/3},{1/2,0,1/2,0},{1/3,1/3,0,1/3},{1/2,0,1/2,0}};
{P.P//MatrixForm,P.P.P//MatrixForm,MatrixPower[P,30]//N//MatrixForm}

An initial distribution is given by a row vector u0.  For the ant chain, u0 = {1,0,0,0} means the ant starts at vertex number 1 with probabiliy 1.  The distribution at time n is obtained by multiplying repeatedly by P on the right or by multiplying u0 by P^n.  This vector also converges to a limit.

u0={1,0,0,0};
u1 = u0.P
u2=u1.P
u10= u0.MatrixPower[P,10]//N
u100= u0.MatrixPower[P,100]//N

Absorbing Chains

If the chain has absorbing states, the matrix powers reveal the probabilities of being absorved at various states given that we start at a certain state.

Here is the basketball shoot-out chain where players A,B,C try to sink the first basket.  The matrix in the upper right corner gives the probabilities for who will win, given who starts.

P={{0,8/10,0,2/10,0,0},{0,0,7/10,0,3/10,0},{5/10,0,0,0,0,5/10},{0,0,0,1,0,0},{0,0,0,0,1,0},{0,0,0,0,0,1}};
P//MatrixForm
MatrixPower[P,100]//N//MatrixForm

These limiting probabilities can also be found using the fundamental matrix N0  and the block R (you can’t call the fundamental matrix N since that is a Mathematica command for giving a numerical approximation)

Q = P[[{1,2,3},{1,2,3}]]; R = P[[{1,2,3},{4,5,6}]];
{Q//MatrixForm,R//MatrixForm}

N0 = Inverse[IdentityMatrix[3]-Q];
{N0//MatrixForm,N0.R//MatrixForm,N0.{1,1,1}}

Ergodic Chains -- Finding the Equilibrium Distribution

For a regular Markov chain, the matrix powers P^n converge to a matrix with equal rows.  The row vector which appears in the rows is the unique equilibrium (fixed) distribution, and it is the limiting distribution for any initial distribution.  There are several ways to find it.

You can approximate the limit by calculating a large matrix power.

P = {{0,1/3,1/3,1/3},{1/2,0,1/2,0},{1/3,1/3,0,1/3},{1/2,0,1/2,0}};
{P//MatrixForm,MatrixPower[P,100]//N//MatrixForm}

To find an exact answer, solve the equation P^T V= V or (P^T - I)V = 0  to get the column vector V = w^T.  The anwer is not unique and it has to be normalized to make a probability vector.  You can just row-reduce the matrix and then solve by hand to get a solution with w4 arbitrary.  The Mathematica Solve gives the answer with w1 arbitrary.  Either way choose the constant to make it a probability vector.

A = Transpose[P]-IdentityMatrix[4];
{A//MatrixForm,RowReduce[A]//MatrixForm}
{w1,w2,w3,w4}/.Solve[A.{w1,w2,w3,w4}=={0,0,0,0}][[1]]

One way to do the normalization automatically is to replace the last row of P^T-I with a vector of all 1’s representing the normalization condition.  If the resulting matrix is called B then we want BV = (0,0,..., 1)^T. The answer is unique.  Here we first define B = P^T-I then alter the last row.

B = Transpose[P]-IdentityMatrix[4];
B[[4]]= {1,1,1,1};
B//MatrixForm
{w1,w2,w3,w4}/.Solve[B.{w1,w2,w3,w4}=={0,0,0,1}][[1]]

A simple alternative (if matrix inverses are available) is to read off the last column of the inverse of B.

Inverse[B]//MatrixForm

The equilibrium vector is also an eigenvector corresonding to eigenvalue lambda = 1.  The eigenvector needs normalization though.

Eigenvalues[Transpose[P]]
Eigenvectors[Transpose[P]]

Here is the two-step weather chain.

P={{1/2,0,1/2,0},{2/5,0,3/5,0},{0,3/5,0,2/5},{0,4/5,0,1/5}};
B = Transpose[P]-IdentityMatrix[4];
B[[4]]= {1,1,1,1};
B//MatrixForm
{w1,w2,w3,w4}/.Solve[B.{w1,w2,w3,w4}=={0,0,0,1}][[1]]
Inverse[B]//MatrixForm

Markov Chains -- Mean Passage Times and Recurrence

Modify the ant on a grid chain to make state 3 absorbing and find the expected times to reach it.

Q={{0,1/3,1/3},{1/2,0,0},{1/2,0,0}}
R= {{1/3},{1/2},{1/2}}
N0 = Inverse[-Q+IdentityMatrix[3]]
N0.{1,1,1}

Similarly, we could make state 2 absorbing.

Q={{0,1/3,1/3},{1/3,0,1/3},{1/2,1/2,0}}
R= {{1/3},{1/3},{0}}
N0 = Inverse[-Q+IdentityMatrix[3]]
N0.{1,1,1}

By symmetry we can find the whole matrix of mean passage times.

M= {{0,4,5/2,4},{9/4,0,9/4,5},{5/2,4,0,4},{9/4,5,9/4,0}};

Find the matrix M+D = PM+C which has the mean passage and recurrence times.

C0=Table[1,{i,1,4},{j,1,4}]
P = {{0,1/3,1/3,1/3},{1/2,0,1/2,0},{1/3,1/3,0,1/3},{1/2,0,1/2,0}};P.M+C0//MatrixForm