Markov Chains

Andy Long
Spring, 2018

Some interesting sports examples from Olinick’s book An Introduction to Mathematical Models in the Social and Life Sciences.

These examples have “absorbing states”, so we have some new and interesting questions to ask about them, like “How long before one falls into an absorbing state, if one starts in state S0?”

Best of Three Games:

p. 296, two teams in best of three series: overall winner is first to two games. States come in two kinds: there are transient states 0-0, 1-0, 0-1,1-1; and there are absorbing states 2-0, 1-2, 2-1, 0-2.

Clear[p,q];
identity=IdentityMatrix[4];
zeros=ConstantArray[0,{4,4}];
rrrr=ConstantArray[0,{4,4}];
qqqq=ConstantArray[0,{4,4}];
rrrr[[2,1]]=p;
rrrr[[3,4]]=q;
rrrr[[4,2]]=p;
rrrr[[4,3]]=q;
qqqq[[1,2]]=p;
qqqq[[1,3]]=q;
qqqq[[2,4]]=p;
qqqq[[3,4]]=q;
(* R represents transitions from transient states to absorbing states: *)
MatrixForm[rrrr]
(* Q represents transitions from transient states to transient states: *)
MatrixForm[qqqq]
(* and we put it all together (with the identity matrix for absorbing states to absorbing states) 
 to create the markov chain, taking states to states: *)
TransitionMatrix=Join[
Join[identity,zeros,2],
Join[rrrr,qqqq,2]
];
MatrixForm[TransitionMatrix]

Some illuminating matrices related to the transition matrix:

ImQ=identity-qqqq;
enn=Inverse[ImQ]; (* The fundamental matrix of the Markov chain *)
MatrixForm[enn]
bee=enn.rrrr;
MatrixForm[bee]

https://www.wolfram.com/mathematica/new-in-9/markov-chains-and-queues/analyze-a-tennis-game.html

proc=With[{p=3/5,q=2/5},DiscreteMarkovProcess[8,({{1,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0},{0,0,1,0,0,0,0,0},{0,0,0,1,0,0,0,0},{0,0,0,0,0,p,q,0},{p,0,0,0,0,0,0,q},{0,0,0,q,0,0,0,p},{0,p,q,0,0,0,0,0}})]];
proc[[1]]=5; (* this says we start in state 5 *)
Graph[proc,VertexSize->{"Scaled",0.07},VertexLabels->Thread[Range[8]->Map[Placed[Style[#,{Bold}],Center]&,states={"2-0","1-2","2-1","0-2","0-0","1-0","0-1","1-1"}]],ImageSize->Large,ImagePadding->25]

PDF[proc[∞],{1,2,3,4,5,6,7,8}]//N
Mean[FirstPassageTimeDistribution[proc,{6,2}]]
MarkovProcessProperties[proc,"TransientTotalVisitMean"]
MarkovProcessProperties[proc,"TransientVisitMean"][[All]]
MarkovProcessProperties[proc,"TransientVisitVariance"][[All]]
p=3/5
q=1-p
MatrixForm[enn]
MatrixForm[enn.ConstantArray[1,4]]
MatrixForm[bee]
N[%]
Clear[p,q]

MarkovProcessProperties[proc]

f[p_]:=2(1+p-p^2) (* expected length of the series, the sum of the first row of enn *)
Plot[f[p],{p,0,1}]
soln=N[Solve[f[p]==(67*2 + 45*3)/(67+45),p]]
 (* Data from USLTA Official Encyclopedia of Tennis, 112 tourneys, 67 two sets, 45 three set *)
prob=p/.soln[[2]]
twoSets=112*(prob^2+(1-prob)^2)
threeSets=112-twoSets

Tennis Chains.

This is a really interesting example, using real data:

Clear[p,q];
identity=IdentityMatrix[6];
zeros=ConstantArray[0,{6,9}];
rrrr=ConstantArray[0,{9,6}];
qqqq=ConstantArray[0,{9,9}];
rrrr[[4,1]]=p;
rrrr[[6,6]]=q;
rrrr[[7,2]]=p;
rrrr[[8,5]]=q;
rrrr[[9,3]]=p;
rrrr[[9,4]]=q;
qqqq[[1,2]]=p;
qqqq[[1,3]]=q;
qqqq[[2,4]]=p;
qqqq[[2,5]]=q;
qqqq[[3,5]]=p;
qqqq[[3,6]]=q;
qqqq[[4,7]]=q;
qqqq[[5,7]]=p;
qqqq[[5,8]]=q;
qqqq[[6,8]]=p;
qqqq[[7,9]]=q;
qqqq[[8,9]]=p;
TransitionMatrix=Join[
Join[identity,zeros,2],
Join[rrrr,qqqq,2]
];
MatrixForm[TransitionMatrix]

Some illuminating matrices related to the transition matrix:

ImQ=IdentityMatrix[9]-qqqq;
enn=Inverse[ImQ];
MatrixForm[enn]
bee=enn.rrrr;
MatrixForm[bee]

p=prob  (* computed from the previous example, and data *)
q=1-p
proc=With[{p=p,q=q},DiscreteMarkovProcess[15,({{1,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,1,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,1,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,1,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,1,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,p,q,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,p,q,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,p,q,0,0,0},{p,0,0,0,0,0,0,0,0,0,0,0,q,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,p,q,0},{0,0,0,0,0,q,0,0,0,0,0,0,0,p,0},{0,p,0,0,0,0,0,0,0,0,0,0,0,0,q},{0,0,0,0,q,0,0,0,0,0,0,0,0,0,p},{0,0,p,q,0,0,0,0,0,0,0,0,0,0,0}})]];
proc[[1]]=7; (* this says we start in state 5 *)
Graph[proc,VertexSize->{"Scaled",0.07},VertexLabels->Thread[Range[15]->Map[Placed[Style[#,{Bold}],Center]&,states={"3-0","3-1","3-2","2-3","1-3","0-3","0-0","1-0","0-1","2-0","1-1","0-2","2-1","1-2","2-2"}]],ImageSize->Large,ImagePadding->25]

PDF[proc[∞],Table[i,{i,1,15}]]//N
MarkovProcessProperties[proc,"TransientTotalVisitMean"]
MarkovProcessProperties[proc,"TransientVisitMean"][[All]]
MarkovProcessProperties[proc,"TransientVisitVariance"][[All]]
MatrixForm[enn]
MatrixForm[enn.ConstantArray[1,9]]
MatrixForm[bee]

In 89 singles championships of the USLTA from 1881 to 1969, there were 
38 three-set matches,
28 four-set matches, and
23 five-set matches.

Thanks, https://mathematica.stackexchange.com/questions/5579/performing-a-chi-square-goodness-of-fit-test

chi[u_]:=Module[{dim,dof,rs,cs,n,full,exp,chis,restbl,tbl,exptbl,pv,res},dim=Dimensions[u];
dof=Times@@(dim-1);
rs=Map[Plus@@#&,u];
cs=Map[Plus@@#&,Transpose[u]];
n=Total[Flatten[u]];
full=Append[MapThread[Append[#1,#2]&,{u,rs}],Join[cs,{n}]];
exp=Outer[Times,rs,cs]/n;
chis=Total[Flatten[(u-exp)^2/exp]];
restbl=Grid[{{"Degrees of Freedom",dof},{"Chi Square Statistic",N[chis,3]},{"p-value",pv}},Alignment->{{Left,"."}},Frame->All];
tbl=Grid[full,Background->{None,None,{{{dim[[1]]+1,dim[[1]]+1},{1,dim[[2]]+1}}->Pink,{{1,dim[[1]]+1},{dim[[2]]+1,dim[[2]]+1}}->Pink}}];
exptbl=Grid[N[exp],Alignment->{".","."}];
pv=N[SurvivalFunction[ChiSquareDistribution[dof],chis]];
res={"ChiSquareStatistic"->N[chis,2],"p-value"->pv,"result"->restbl,"table"->tbl,"expectedtable"->exptbl,"fullresults"->Column[{"Data Table",tbl,"Expected Values Table",exptbl,"Hypothesis Testing",restbl},Frame->All,Background->{Yellow,None,Yellow,None,Yellow,None}],"Properties"->{"ChiSquareStatistic","p-value","result","table","expectedtable","fullresults"}};
#/.res&]
chisqt[u_,r_]:=chi[u][r]

x=Transpose[{89*{bee[[1,1]]+bee[[1,6]],bee[[1,2]]+bee[[1,5]],bee[[1,3]]+bee[[1,4]]},
{38,28,23}}]
chisqt[x,"table"] 
chisqt[x,"expectedtable"] 
chisqt[x,"result"]
chisqt[x,"fullresults"] y
chisqt[x,"Properties"]

A test against known results -- it checks out....

chisqt[Transpose[{{116,74,59},{84,126,41}}],"fullresults"]