We have finally our transition matrix, found out that the game always succeeds, and now we should try to find how long it will take. We can’t obviously give a precise answer (512 steps or whatever) but we can estimate what is the probability that it will be within a certain amount of steps.

If you read the wikipedia article, you should know that being A a transition matrix, A² will also be a transition matrix. Its (i,j) element will give the probabiliy to go from state i to state j in two steps. In general, A^n (A to the power of n) is a transition matrix that will give the probabilities of going from a state to another in n steps.

We can therefore compute the probabilities that it will take to go from any given state to the end state (square 63) in n steps simply by taking the last column of A^n. The first element of this array in particular will give the probability of going from the start (square 0) to the end in n steps.

Lets take our A matrix, this element can be simply computed as:

countFreq 63 0 10000

(10000 being arbitrarily chosen, pick any number that won’t cause your computer to bend to its knees), which gives 0.1099 (very close to the theoretical probability of 0.111). It is clear that we are not very near to the end of the game yet.

To go further we need to be able to compute the nth power of our matrix A. Instead of reinventing the wheel we will just use the Array class, available from Data.Array. We need to modify our buildTM function so that the matrix will be constructed in accordance with the class definition:

buildTM' :: Double -> Array (Int, Int) Double
buildTM' n = array ((0,0),(63,63)) [((st,sq), countFreq sq st n) | st <- [0..63], sq <- [0..63]]

We are also not going to reinvent the matrix multiplication function, there is already a nice one given here:

matMult :: Array (Int,Int) Double -> Array (Int,Int) Double
	   -> Array (Int,Int) Double
matMult x y =  array ((0,0),(63,63))
                     [((i,j), sum [x!(i,k) * y!(k,j) |k <- range (0,63)])
                                       | i <- range (0,63),
                                         j <- range (0,63)]

Now we can compute A² simply by doing this in ghci:

let a = buildTM' 10000 -- (again feel free to experiment, this takes about 15 minutes on my PC)
let a2 = matMult a a

And the last column of it will be [a2!(i,63)|i<-[0..63]]. The first element of this column will be the probability of going from the first square to the last square in 2 steps:

a2!(0,63) -- 0.1099

Well, same result as before. Were would we expecting anything different!?
Now, increasing the exponent should start giving some results.

let a4 = matMult a2 a2
let a8 = matMult a4 a4
let a16 = matMult a8 a8
let a32 = matMult a16 a16

And now:

a32!(0,63) -- 0.7189049811868098

So the probability of ending the game in less than 32 steps is already quite high, about 72%.

let a64   = matMult a32 a32
let a128 = matMult a64 a64
let a256 = matMult a128 a128
 
a64!(0,63)  -- 0.9292194901730458
a128!(0,63) -- 0.9955122418509581
a256!(0,63) -- 0.9999819589853614

Remember that previously we found that the longest sequence in 10000 simulations was 219? As we see, the probability of ever finding a sequence longer than 256 steps is pretty low (I don’t know what you think, but a probability of 0.0018% is pretty near to “impossible” in my vocabulary). Just for curiosity, what is the probability of finding a sequence longer than 219 steps? We can write 219 as 128+64+16+8+2+1 so using our already computed matrices we easily find that:

a219!(0,63) -- 0.9999111233274545

which is indeed very close to 9999/10000 …

Do you remember the clumsy countFreq function I wrote? It turns out that the Array class provides a built-in function to do exactly this (accumArray). We can therefore cut the middle man and build our transition matrix directly using it:

countFreq :: Int -> Int -> Int -> Array Int Int
countFreq st n s = accumArray (+) 0 (0,63) [(i,1)| i <- runSim (n-1) st s]

buildTM :: Int -> Int -> Array (Int,Int) Double
buildTM n s = array ((0,0),(63,63)) [((i,j),
	fromIntegral ((countFreq i n s)!j) / fromIntegral n) |
	i <- range (0,63),
	j <- range (0,63)]

I don’t know how smart ghc is, but I expect that the countFreq in (countFreq i n s)!j) will be evaluated only once (somehow I think there should be a simpler method to combine arrays into a matrix). Both methods seem to be equivalent in terms of timing.

Assuming that 4 players are playing and that each take 15 sec per turn, it seems pretty likely that few games (if any) will last longer than half hour (since we have 4 kids, the single probabilities have to be combined). Since the probability of having 4 kids playing together in a room for half an hour without a fight breaking is pretty low, I think we can safely enunciate the Goose theorem:

It is very probable that 4 kids playing the Game of the Goose will have a fight

And this i can definitively prove from personal experience.