A visual guide to Markov Chains

Introduction

Mathematica has built in functions to deal with Markov Chains. 
Consider the following transition matrix:

Q = (({{(1/2), (1/3), 0, (1/5), 0}, {0, (2/3), (1/3), (1/5), (1/2)}, {(1/2), 0, (1/3), (1/5), 0}, {0, 0, (1/3), (1/5), 0}, {0, 0, 0, (1/5), (1/2)}}));

And suppose we want to start our state in

p0 = {0,1,0,0,0};

With these two ingredients, we may use the function “DiscreteMarkovProcess” to do all sorts of manipulations of Markov Chains.

mark= DiscreteMarkovProcess[p0,Q\[Transpose]]

But be careful: Mathematica describes the chain such that the rows add up to 1, not the columns. So when we feed it the matrix, we must feed Qᵀ

We can now visualize the matrix using the following function

Clear[MarkovGraph];
MarkovGraph[mark_]:=Module[{g},
g = Graph[mark];
Scan[(PropertyValue[{g,#},EdgeLabels]=PropertyValue[{g,#},"Probability"])&,EdgeList[g]];
(*Show[g,ImageSize->600]*)
g
];

MarkovGraph[mark]

Running a simulation

Now let us run a simulation of our Markov Graph.
To do that we use the function “RandomFunction”

data=RandomFunction[mark,{0,50}]

Mathematica returns a “TemporalData” object. 
This is useful because now you can extract properties of your data using built-in functions. 
If you want you can extract the full data using “Normal”

First@Normal@data

This is a list of tuples (t,(Subscript[X, t])). We can plot it if we want:

But Mathematica’s functions understand TemporalData objects. 
So we can use it, for instance, to make plots.

ListPlot[data,Filling->Axis,FrameTicks->{{{1,2,3,4,5},None},{Automatic,None}},PlotRange->{0,5.4},
GridLines->{None,Range[1,5]},
Frame->True,
FrameLabel->{"t","(Subscript[X, t])"},
BaseStyle->{FontFamily->"Times",18},
ImageSize->500
]

Or we can visualize this as follows:

Animate[Graph[mark,
PlotLabel->"t = "<>ToString@i,
VertexStyle->{(First@Normal@data)[[i,2]]->Darker@Green}],{i,1,Length@First@Normal@data,1}]

Each time we repeat the experiment we get a new Stochastic Trajectory:

Table[
ListPlot[RandomFunction[mark,{0,30}],Filling->Axis,Ticks->{Automatic,{1,2,3,4,5}},PlotRange->{0,5.4},
GridLines->{None,Range[1,5]},
Frame->True,
FrameLabel->{"t","(Subscript[X, t])"},
BaseStyle->{FontFamily->"Times",18},
ImageSize->230
],{6}]

Repeating the experiment many times

Now we repeat the experiment several times, always starting at the same initial condition:

data2 = RandomFunction[mark,{0,500},(10^4)]

We can use this, for instance, to compute means:

mean=TimeSeriesThread[Mean,data2]

ListLinePlot[mean,PlotRange->{1,5},FrameLabel->{"t","⟨(Subscript[X, t])⟩"}]

Or variances:

ListLinePlot[TimeSeriesThread[Variance,data2],PlotRange->All]

At t = 0 we started with the probability vector (Subscript[p, 0]) = (1,0,0,0,0). 
This means we had 100% certainty the system was at 1. 
Hence, at t = 0 the mean should be 1 and the variance should be 0.

We can also visualize the distribution as a function of time:

DistributionChart[Table[data2["SliceData",t],{t,0,10}],
FrameLabel->{"t","(Subscript[X, t])"},
BaseStyle->{FontFamily->"Times",18},
FrameTicks->{{{1,2,3,4,5},None},{Automatic,None}}
]

It helps to compare this chart with our graph:

MarkovGraph[mark]

Let us look at the distribution after a long time has passed:

Histogram[data2["SliceData",50],{0.5,5.5,1},"Probability"]

This is the Steady-State.
We can also compare it with the analytical calculation.

Steady-state

We first compute the eigenvalues of Q

N[Eigenvalues[Q]]

We see that the chain has only one eigenvalue = 1
The corresponding eigenvector is the first eigenvector

ps = First@Eigenvectors[Q]

But this is not yet normalized:

Total@ps

So we normalize it as

ps = (ps/(Total@ps))

We may now compare this with our simulations

Show[
Histogram[data2["SliceData",50],{-0.5,5.5,1},"PDF",PlotRange->{{0.5,5.5},All}],
ListPlot[ps,Filling->Axis,Ticks->{Automatic,{1,2,3,4,5}}]
]

It is already pretty close.

Unused functions

tmp=Tally[data2["SliceData",10]][[All,2]];
tot = Total@tmp;

tally[slice_]:=Module[{tal},
tal = Tally[slice];
Do[If[Not@MemberQ[tal[[All,1]],i],AppendTo[tal,{i,0}]],{i,1,5,1}];
Sort@tal
]

Graphics[{
Flatten@Table[
tmp=tally[data2["SliceData",t]][[All,2]];
tot = Total@tmp;
{PointSize[0.1 ((tmp[[i]])/tot)],Point[{t,i}]},{t,0,10,1},{i,1,5}]
},
Frame->True,
FrameLabel->{"t","(Subscript[X, t])"},
BaseStyle->{FontFamily->"Times",18},
FrameTicks->{{{1,2,3,4,5},None},{Automatic,None}}
]

Other examples of Markov Chains

Gambler’s ruin

Two gamblers playing until they run out of money.
 p = probability of winning
 n = total amount of money the two have

GamblersRuin[p_,n_]:=SparseArray[{{1,1}->1,{n+1,n+1}->1,{i_,j_}/;1<i<n+1&&j==i+1->p,{i_,j_}/;1<i<n+1&&j==i-1->1-p},{n+1,n+1}]\[Transpose];
ntot = 7;
p = 0.5;
GamblersRuin[p,ntot]//MatrixForm

gamb = DiscreteMarkovProcess[4,GamblersRuin[p,ntot]\[Transpose]]

Show[MarkovGraph[gamb],ImageSize->200]

SeedRandom[2347];
ListPlot[{#[[1]],#[[2]]-1}&/@RandomFunction[gamb,{0,30}]["Path"],Filling->Axis]

We can compute the “First Passage Time” , which is the average time it takes for the Gambler to go broke.

Mean[FirstPassageTimeDistribution[gamb,ntot+1]]

Or the probability of winning/loosing

DiscretePlot[PDF[gamb[Infinity],k],{k,1,ntot+1},PlotRange->All]

PDF[gamb[Infinity],ntot+1]
PDF[gamb[Infinity],1]

A horse on a Chess board

gr=KnightTourGraph[8,8,EdgeStyle->LightGray]

Mathematica understands a random walk on a graph by simply giving it the graph

knightmark = DiscreteMarkovProcess[1,gr]

Let’s do a simulation

sim=RandomFunction[knightmark,{0,100}];
HighlightGraph[gr,Subgraph[gr,sim["ValueList"]]]

We can also get, for instance, the stationary distribution

steadystate=PDF[StationaryDistribution[knightmark],Range[64]]

Recall that (Subscript[p, n]^( * )) = 1/(Subscript[τ, n]), where (Subscript[τ, n]) is the average time the system spends in each state. 
Thus

MatrixForm[1/Partition[steadystate,8]]

The horse spends more time on the corners, instead of the center.