(This article was first published on Freakonometrics - Tag - R-english, and kindly contributed to R-bloggers)
Markov chains is a very interesting and powerful tool. Especially for parents. Because if you think about it quickly, most of the games our kids are playing at are Markovian. For instance, snakes and ladders...It is extremely easy to write down the transition matrix, one just need to define all snakes and ladders. For the one above, we have,
n=100 M=matrix(0,n+1,n+1+6) rownames(M)=0:n colnames(M)=0:(n+6) for(i in 1:6){diag(M[,(i+1):(i+1+n)])=1/6} M[,n+1]=apply(M[,(n+1):(n+1+6)],1,sum) M=M[,1:(n+1)] starting=c(4,9,17,20,28,40,51,54,62, 64,63,71,93,95,92) ending =c(14,31,7,38,84,59,67,34,19, 60,81,91,73,75,78) for(i in 1:length(starting)){ v=M[,starting[i]+1] ind=which(v>0) M[ind,starting[i]+1]=0 M[ind,ending[i]+1]=M[ind,ending[i]+1]+v[ind]}So, why is it important to have a Markov Chain ? Because, once you've noticed that you had a Markov Chain game, you can derive anything you want. For instance, you can get the distribution after some turns,
powermat=function(P,h){ Ph=P if(h>1){ for(k in 2:h){ Ph=Ph%*%P}} return(Ph)} initial=c(1,rep(0,n)) COLOR=rev(heat.colors(101)) u=1:sqrt(n) boxes=data.frame( index=1:n, ord=rep(u,each=sqrt(n)), abs=rep(c(u,rev(u)),sqrt(n)/2)) position=function(h=1){ D=initial%*%powermat(M,h) plot(0:10,0:10,col="white",axes=FALSE, xlab="",ylab="",main=paste("Position after",h,"turns")) segments(0:10,rep(0,11),0:10,rep(10,11)) segments(rep(0,11),0:10,rep(10,11),0:10) for(i in 1:n){ polygon(boxes$abs[i]-c(0,0,1,1), boxes$ord[i]-c(0,1,1,0), col=COLOR[min(1+trunc(500*D[i+1]),101)], border=NA)} text(boxes$abs-.5,boxes$ord-.5, boxes$index,cex=.7) segments(c(0,10),rep(0,2),c(0,10),rep(10,2)) segments(rep(0,2),c(0,10),rep(10,2),c(0,10))}Here, we have the following (note that I assume that once 100 is reached, the game is over)


> h=10 > (initial%*%powermat(M,h))[59:61]/ + sum((initial%*%powermat(M,h))[59:61]) [1] 0.1597003 0.5168209 0.3234788
distrib=initial%*%M game=rep(NA,1000) for(h in 1:length(game)){ game[h]=distrib[n+1] distrib=distrib%*%M} plot(1-game[1:200],type="l",lwd=2,col="red", ylab="Probability to be still playing")

> sum(1-game) [1] 32.16499
But assuming that you are playing with your daughter, and that the game is over once one player reaches the 100 cell, it is possible to get the survival distribution of the first time one of us reaches the 100 cell,
plot((1-game[1:200])^3,type="l",lwd=2,col="blue", ylab="Probability to be still playing (2 players)")

> sum((1-game)^2) [1] 23.40439And if you ask your son to join the game, the survival distribution function is
plot((1-game[1:200])^3,type="l",lwd=2,col="purple", ylab="Probability to be still playing (3 players)")

> sum((1-game)^3) [1] 20.02098
To leave a comment for the author, please follow the link and comment on his blog: Freakonometrics - Tag - R-english.
R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...


Zero Inflated Models and Generalized Linear Mixed Models with R.
Zuur, Saveliev, Ieno (2012).