(This article was first published on Freakonometrics » R-english, and kindly contributed to R-bloggers)
A couple of days ago, we had a quick chat on Karl Broman‘s blog, about snakes and ladders (see http://kbroman.wordpres...
(This article was first published on Freakonometrics » R-english, and kindly contributed to R-bloggers)
A couple of days ago, we had a quick chat on Karl Broman‘s blog, about snakes and ladders (see http://kbroman.wordpress.com/…) with Karl and Corey (see http://bayesianbiologist.com/….), and the use of Markov Chain. I do believe that this application is truly awesome: the example is understandable by anyone, and computations (almost any kind, from what we’ve tried) are easy to perform. At the same time, some French students asked me specific details regarding some old lectures notes on Markov chains, and on some introductory example I used as a possible motivation: the stepping stone algorithm. In the notes, I just mentioned the idea of this popular generic algorithm (introduced in Sawyer (1976)) and I use simulations to show - visually - how it works. Again, it was just to motivate the course which actually did focus on the theory of Markov Chains. But those student wanted more, like how did I get the transition matrix, for instance. And that is actually not a simple question, from a computational perspective. I mean, I can easily generate this Markov Chain, but writing explicitly the transition, that was another story. Which took me a bit longer. In a very specific case…
But let us get back to the roots, and to the stepping stone algorithm. At least, one of them (the one I used in my notes) because it looks like there are several algorithm. We do consider a grid, say , with some colors inside, say possible colors. Each cell of the grid has a given color. Then, at some stage, we select randomly one cell in the grid, and it will take the color of one of its neighbor (some kind of absorption, or mutation). This is, more or less, what is also detailed in some lecture notes by James Propp (see also e Sato (1983) or Zähle et al. (2005) for more theoretical details about that Markov chain). This is extremely simple to generate (that’s what I did in my notes, with very big grids, and a lot of colors). But what if we want to write the transition matrix ?
First of all, we need to define the state space. Basically, we do have cells, each of them has one color, chosen among . Which gives us possible states…. And that can be large. I mean, if we consider the smallest possible grid (that might be interesting), say , and only colors, then we talk about possible states. That is large, not huge. But we should keep in mind that we have to compute a transition matrix, that would be a matrix with elements. More generally, we talk about writing down matrices with elements. If we want black and white grids, that would mean a matrix with which mean 4 billion elements ! And if we consider an red-green-blue grid, we have to explicit a matrix with i.e almost 400 million elements. So, let’s face it: we can only work with bi-color grids.
So let’s try… The good thing is that it can be related to work I’ve been doing recently on binomial recombining trees (binomial being related to bi-color). First of all, our grid will be describes as follows
> h=3
> M=matrix(1:(h^2),h,h)
> M
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
with two colors
> color=c("red","blue")
Then, we should look for neighbors, or derive an neighborhood matrix,
> d=function(i,j) dist(rbind(c((i-1)%/%h,(i-1)%%h),
+ c((j-1)%/%h,(j-1)%%h)))
> Neighb=matrix(Vectorize(d)(rep(1:(h^2),each=h^2),
+ rep(1:(h^2),h^2)),h^2,h^2)
> trunc(Neighb*100)/100
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.00 1.00 2.00 1.00 1.41 2.23 2.00 2.23 2.82
[2,] 1.00 0.00 1.00 1.41 1.00 1.41 2.23 2.00 2.23
[3,] 2.00 1.00 0.00 2.23 1.41 1.00 2.82 2.23 2.00
[4,] 1.00 1.41 2.23 0.00 1.00 2.00 1.00 1.41 2.23
[5,] 1.41 1.00 1.41 1.00 0.00 1.00 1.41 1.00 1.41
[6,] 2.23 1.41 1.00 2.00 1.00 0.00 2.23 1.41 1.00
[7,] 2.00 2.23 2.82 1.
