One of the simplest non-trivial examples of something interesting to do with Markov chains (that’s pretty easy to visualise and understand) is working out how many moves (on average) it takes to complete a game of snakes and ladders? This selection of an example is somewhat intentional because it’s the same example as the wikipedia entry for absorbing Markov chains, and because it’s just hard enough to actually justify solving the problem properly.

A full description of absorbing Markov chains is available on wikipedia, which is a near verbatim transcription from Kemeny and Snell’s book (also available online). In short, a (Discrete Time, Finite State) Markov chain is a system that describes the evolution of probability of going from one state to another in an instantaneous transition. Transitions between states are modelled using a stochastic matrix (a matrix where the rows must sum to one).

The wikipedia article on snakes and ladders quotes the solution as 39.6, for the Milton Bradley version of Chutes and Ladders. That’s a useful baseline for comparison, because the board layout is easy to find online. It’s also worth noting that the page for absorbing Markov chains quotes the solution as 39.2.

Each of the states in the chain represents a location on the board, and we’ll have one extra state (state 0) that represents the initial state off the board. That means we can conveniently number those state 0 - 100 (inclusive), where 0 is the initial state off the board, and 100 is the terminal state.

Ignoring (initially) the chutes and ladders, a transition from a state results in you advancing by the number of squares on a six sided dice. From any state $x$, you can reach states $x+1, x+2, …, x+6$, unless such a transition would take you off the board (beyond state 100), where the transition is invalid (which means you stay where you are). The probability of each of these transitions is ¹⁄₆, except for the edge case where the move is invalid - in those cases a small adjustment is needed to keep the player in that states.

def snakes_and_ladders_matrix():
    n = 101
    A = np.zeros((n, n))
    for i in xrange(6):
        A += np.diag(np.ones(n-1-i),i+1)

    # work out the cells that are reachable for every initial cell 
    b = np.sum(A, axis=1)

    # update the transition matrix to have a self-edge for every missing edge 
    A += np.diag(6 - b)

    # normalize
    A /= 6.
    
    return A

A = snakes_and_ladders_matrix()

We can verify that the resulting matrix is a valid (right) stochastic matrix trivially - and it’s always a good idea to enforce constraints like this when working with Markov Chains.

np.testing.assert_allclose( 1, np.sum(A, axis=1) )

So for the game with no snakes or ladders, this is enough to show how many moves it take to finish the game. The equation for mean number of steps can be written in terms of the fundamental matrix $N$ and the identity matrix $I$:

$$ N = (I-Q)^{-1} $$

The equation is then

$$ t = N \mathbf{1} $$

Where $\mathbf{1}$ is the column vector of ones. While it’s possible to solve the equation like this, but only just, it’s not an effective way to use a computer to solve the equation. Solutions to linear equations on a computer are solved in the form:

$$ A x = b $$

But our unknown is on the wrong side of the equation. Somewhat conveniently, multiplying through by $I-Q$ gives:

$$ (I-Q) t = \mathbf{1} $$

This means we can solve the equation:

Q = A[:-1,:-1]
I = np.identity(Q.shape[0])
o = np.ones((Q.shape[0],1))
result = np.linalg.solve(I-Q, o)
print 'mean turns to complete (without snakes and ladders)', result[0, 0]

Which gives 33.3, which is close enough to give confidence that we’re in the ballpark when we add the snakes and ladders.

ladders = np.array([
    [1, 38],  [4, 14],  [9, 31],  [21,42],
    [28,84],  [36,44],  [51,67],  [71,91],
    [80,100]
]) 

snakes = np.array([
    [98, 78],  [95, 75],  [93, 73],  [87, 24],
    [64, 60],  [62, 19],  [56, 53],  [49, 11],
    [48, 26],  [16, 6],
]) 

def snakes_and_ladders_matrix():
    n = 101
    A = np.zeros((n, n))
    for i in xrange(6):
        A += np.diag(np.ones(n-1-i),i+1)

    # work out the cells that are reachable for every initial cell 
    b = np.sum(A, axis=1)

    # update the transition matrix to have a self-edge for every missing edge 
    A += np.diag(6 - b)

    edges = np.vstack([ladders, snakes])

    for src, dst in edges:
        b = A[:, src].copy()
        A[:, dst] += b
        A[:, src] = 0
    
    A /= 6.
    return A


A = snakes_and_ladders_matrix()
assert_allclose( 1, np.sum(A, axis=1) )

And running the code gives the solution

Q = A[:-1,:-1]
I = np.identity(Q.shape[0])
o = np.ones((Q.shape[0],1))
result = np.linalg.solve(I-Q, o)
print 'mean turns to complete', result[0, 0]

Gives 39.598