Markov Process Steady State (Eigenvector) Determination

Exercize to calculate the steady state (eigenvector) of a Markov process (reflecting a matrix made from your student ID) in three different ways : 
  1. Pump any starting probability distribution vector into the transition matrix, and iteratively feedback the output back into the transition.
  2. Raise the transtion matrix to a high enough power such that the columns have fixed value.
  3. Solve analytically, using the side condition that the sum of the distribution is unity.

Preliminary initialization of transition matrix.

Get the student ID from an escape to the (Unix) operating system.

sID = <<"!whoami"

sID = "s1231234";

Strip off the leading "s".
Parse as number.
Put extra digits before and after.

digits = 
Append[
Prepend[
IntegerDigits[
ToExpression[
StringDrop[sID, 1]
]
],
9],
9]

Compose a matrix with rows normalized (sum to unity).

(The //MatrixForm operator cannot be part of the assignment, since it changes the matrix form, so invoke it after assignment to the variable.)

matrix0 = {
(({digits[[1]],digits[[2]],digits[[3]]})/(digits[[1]] + digits[[2]] + digits[[3]])),
(({digits[[4]], digits[[5]], digits[[6]]})/(digits[[4]] + digits[[5]] + digits[[6]])),
(({digits[[7]], digits[[8]], digits[[9]]})/(digits[[7]] + digits[[8]] + digits[[9]]))
}

%//MatrixForm

1. Pump an (arbitrary) probability distribution into this transition matrix, feeding the output back in repeatedly until a "fixed point" is apparently reached.

Define a vector as a list:

p = {pa,pb, pc}

To derive future state given current state, regular multiplication isn't what we want:

p matrix0//MatrixForm

But the dot product is (collapsing the columns as sums), using left-multiply:

p . matrix0

Right-multiply should not be used, as it yields something different:

matrix0 . p

Now iteratively multiply an arbitrary probability distribution times this transition matrix, feeding back the results as successive inputs. (“//N” means to convert the output to a real number, instead of the default rational [fractional] representation.)

{1,0,0}.matrix0 //N

("%" means the results of the last calculation.)

%.matrix0 //N

%.matrix0 //N

%.matrix0 //N

%.matrix0 //N

%.matrix0 //N

2. The steady state of an ergodic process does not depend on initial conditions, so let' s try another technique that doesn’t use a start vector.

matrix1 = {{a, b}, {c, d}}

Raising a matrix to the second power is the same as taking the "dot product" with itself.

matrix1 //MatrixForm

MatrixPower[matrix1, 2] //MatrixForm

matrix1 . matrix1 //MatrixForm

What power must your transition matrix ("matrix0") be raised to in order to have its columns stabilized?

(Please ignore the warnings about the deprecated package.)

<<BarCharts`;

StackedBarChart[
Transpose[
Table[{1,0,0} . MatrixPower[matrix0,i],{i,0,10}]
],
BarStyle->{Red, Green, Blue}
]

Since system is ergodic, steady state solution is independent of original conditions.

StackedBarChart[
Transpose[
Table[{0,1,0} . MatrixPower[matrix0,i],{i,0,10}]
],
BarStyle->{Red, Green, Blue}
]

StackedBarChart[
Transpose[
Table[{0,0,1} . MatrixPower[matrix0,i],{i,0,10}]
],
BarStyle->{Red, Green, Blue}
]

Alternatively, dispense with original state, and just exponentiate transition matrix (of ergodic process) until it stabilizes.

MatrixPower[matrix0, 2]//N //MatrixForm

MatrixPower[matrix0, 5]//N //MatrixForm

MatrixPower[matrix0, 10] //N //MatrixForm

MatrixPower[matrix0, 20]//N //MatrixForm

StackedBarChart[
Transpose[
Table[MatrixPower[matrix0,i][[1]],{i,0,10}]
],
BarStyle->{Red, Green, Blue}
]

Note how processes with a recurrence (non-unity loop) never settle:

matrix2 = {{0,1},{1,0}}

matrix2 //MatrixForm

StackedBarChart[
Transpose[
Table[{1,0} . MatrixPower[matrix2,i],{i,0,9}]
],
BarStyle->{Red, Green}
]

matrix3 ={{0,1,0},{0,0,1},{1,0,0}};

matrix3 //MatrixForm

GraphPlot[matrix3,
VertexLabeling->True, DirectedEdges->True
]

MatrixPower[matrix3,1] //MatrixForm

MatrixPower[matrix3,2] //MatrixForm

MatrixPower[matrix3,3] //MatrixForm

MatrixPower[matrix3,4] //MatrixForm

StackedBarChart[
Transpose[
Table[{1,0,0} . MatrixPower[matrix3,i],{i,0,9}]
],
BarStyle->{Red, Green, Blue}
]

But a leak in just one node will allow equilibrium:

matrix4 = 
{
{0,1/2,1/2},
{0,0,1},
{1,0,0}
};

%//MatrixForm

GraphPlot[matrix4,
VertexLabeling->True, DirectedEdges->True
]

MatrixPower[matrix4,45] //MatrixForm //N

Apply this "MatrixPower" technique to our own transition matrix, observing twenty iterations.

StackedBarChart[
Transpose[
Table[{1,0,0} . MatrixPower[matrix4,i],{i,0,19}]
],
BarStyle->{Red, Green, Blue}
]

3. Alternatively, solve the system analytically.

Since the matrix rows are not independent (as they always sum to unity), the determinant of the characteristic equation system vanishes.

Det[matrix0 - IdentityMatrix[3]]

Find the equilibrium solution (Oishi, p. 50):

Unconstrained system is under-determined (since the matrix is singular),
so has no unique solution.

Solve[
{
p . matrix0== p
},
p
]

We need to include extra constraint that p is probability distribution (which adds up to unity):

Solve[
{
p . matrix0== p,
pa + pb + pc == 1
},
p
]

Actually, only 2/3 of the matrix is needed, since there is a linear dependence.

Solve[
{
(p  . matrix0)[[1]]== pa,
(p  . matrix0)[[2]]== pb,
pa + pb + pc == 1
}, 
p
]

This result should be the same as those obtained iteratively above.

% //N

StackedBarChart[
Transpose[
p /.%
],
BarStyle->{Red, Green, Blue}
]