Simulating Chutes & Ladders in Python

Edit 12/19/2017: added a new subsection on analyzing Chutes & Ladders as an Absorbing Markov Chain.

[img: Chutes and Ladders animated simulation]

This weekend I found myself in a particularly drawn-out game of Chutes and Ladders with my four-year-old. If you've not had the pleasure of playing it, Chutes and Ladders (also sometimes known as Snakes and Ladders) is a classic kids board game wherein players roll a six-sided die to advance forward through 100 squares, using "ladders" to jump ahead, and avoiding "chutes" that send you backward. It's basically a glorified random walk with visual aids to help you build a narrative. Thrilling. But she's having fun practicing counting, learning to win and lose gracefully, and developing the requisite skills to be a passionate sports fan, so I play along.

On the approximately twenty third game of the morning, as we found ourselves in a near endless cycle of climbing ladders and sliding down chutes, never quite reaching that final square to end the game, I started wondering how much longer the game could last: what is the expected length of a game? How heavy are the tails of the game length distribution? How succinctly could I answer those questions in Python? And then, at some point, it clicked: Chutes and Ladders is memoryless — the effect of a roll depends only on where you are, not where you've been — and so it can be modeled as a Markov process! By the time we (finally) hit square 100, I basically had this blog post written, at least in my head.

When I tweeted about this, people pointed me to a number of similar treatments of Chutes & Ladders, so I'm under no illusion that this idea is original. Think of this as a blog post version of a dad joke: my primary goal is not originality, but self-entertainment, and if anyone else finds it entertaining that's just an added bonus.

Direct Simulation

The most straightforward way to get a handle on the dynamics of the game is through direct simulation: if we simulate enough games, we'll obtain a distribution of game lengths that will approximate the "true" distribution. The first step in this is to examine the game board and somehow encode the positions of the chutes and ladders on the grid:

[img: Chutes and ladders game board]

(Image source: uncyclopedia)

We'll use a Python dictionary to store these positions:

In [1]:
# Mapping of start : end spaces of chutes & ladders
CHUTES_LADDERS = {1:38, 4:14, 9:31, 16:6, 21:42, 28:84, 36:44,
                  47:26, 49:11, 51:67, 56:53, 62:19, 64:60,
                  71:91, 80:100, 87:24, 93:73, 95:75, 98:78}

With this in place, we can simulate the game in a few lines of Python:

In [2]:
from random import Random


def simulate_cl_game(rseed=None, max_roll=6):
    """
    Simulate a full Chutes & Ladders game
    and return the number of turns to finish
    """
    rand = Random(rseed)
    position = 0
    turns = 0
    while position < 100:
        turns += 1
        roll = rand.randint(1, max_roll)
        
        # if the roll takes us past square 100, we don't move
        if position + roll > 100:
            continue
            
        # otherwise, move the position according to the roll
        position += roll
        
        # go up/down any chute/ladder
        position = CHUTES_LADDERS.get(position, position)
    return turns

Calling the function tells us how many moves were required to finish the particular game:

In [3]:
simulate_cl_game()
Out[3]:
29

If we simulate many games, the result will be a distribution of how many turns are required to get to the end:

In [4]:
# Jupyter notebook plotting setup & imports
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn')
In [5]:
sim_games = [simulate_cl_game() for i in range(10000)]

plt.hist(sim_games, bins=range(200), normed=True)
plt.xlabel('Number of Turns to win')
plt.ylabel('Fraction of games')
plt.title('Simulated Lengths of Chutes & Ladders Games');

This gives us some insight, but the problem here is that the result is just an estimate of the "true" distribution; to make our estimate more precise will require many more simulations, and this can get computationally expensive.

Fortunately, there's another approach.

Chutes and Ladders as a Markov Process

Instead of brute force simulation, we might think about the game probabilistically. On any given turn, there are six equally probable options: rolling a 1, 2, 3, 4, 5, or 6. Depending on which space you start on, these lead to six well-defined results. For example, the first turn, the possibilities are the squares 38, 2, 3, 14, 5, or 6, each with equal probability. We could encode this set of probabilities as a vector of length 101, with 1/6 in each associated index (here the zeroth element represents the start of the game, off of the board):

0: [0, 0, 1/6, 1/6, 0, 1/6, 1/6, 0, 0, ...]

Each entry in this vector is the probability of going from square zero to the corresponding square. This vector completely describes the first turn of the game.

Similarly, we could construct the vector describing turns from square 2:

2: [0, 0, 0, 1/6, 0, 1/6, 1/6, 1/6, 1/6, 0, 0, ...]

Again, this completely describes any turn of the game that starts at square two.

The key insight is that Chutes and Ladders is memoryless; if you're on, say, square 19, it doesn't matter if you got there by sliding down from square 62, or by rolling a five from square 14 – the probabilities on the next turn are exactly the same.

This kind of situation: a sequence of probabilistic transitions from one state to another with no memory of previous states, is known as a Markov Chain or Markov Process, and can be entirely described by a N x N transition matrix where N is the number of states in the system (101 states for Chutes and Ladders). The columns of the matrix contain the probability vectors we discussed above, such that element (n, m) is the probability of transitioning to element n when starting from element m.

Let's create a Markov transition matrix for Chutes and Ladders:

In [6]:
import numpy as np

def cl_markov_matrix(max_roll=6, jump_at_end=True):
    """
    Create a Markov transition matrix
    
    If jump_at_end is True, then apply ladder/chute jumps at the end of each turn.
    If False, then apply them at the beginning of the next turn.
    """  
    # Create the basic transition matrix:
    mat = np.zeros((101, 101))
    for i in range(101):
        mat[i + 1:i + 1 + max_roll, i] = 1. / max_roll
        
    # We could alternatively use scipy.linalg.circulent as follows:
    # mat = circulant([0, *np.ones(max_rolls) / 6, *np.zeros(100)])[:101, :101]

    # rolls off the end of the board don't change the state;
    # add these probabilities to the diagonal
    mat[range(101), range(101)] += 1 - mat.sum(0)

    # account for the presence of chutes and ladders
    # we'll do this via  another transition matrix
    cl_mat = np.zeros((101, 101))
    ind = [CHUTES_LADDERS.get(i, i) for i in range(101)]
    cl_mat[ind, range(101)] = 1
    if jump_at_end:
        return cl_mat @ mat
    else:
        return mat @ cl_mat

mat = cl_markov_matrix()
plt.matshow(mat)
plt.grid(False)

The above matrix encodes everything you need to know about the game of Chutes and Ladders. (Note that there are two ways to construct this matrix: applying the chutes and ladders jumps at the end of each turn, as typically played, or applying them at the beginning of the next turn. We make this second option available for later).

The game starts with a state vector describing 100% probability of being in the zeroth state:

v_0 = [1, 0, 0, 0, 0...]

If you think about it for a while, you can convince yourself that the way to evolve this (probabilistically) to the next state is to left-multiply the vector by our transition matrix:

In [7]:
np.set_printoptions(suppress=True, precision=2)

v_0 = [1, *np.zeros(100)]
mat @ v_0
Out[7]:
array([ 0.  ,  0.  ,  0.17,  0.17,  0.  ,  0.17,  0.17,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.17,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.17,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ])

The result is precisely the vector we described above, encoding the probabilities of game state after one move.

To perform multiple moves, we could repeatedly multiply this vector by our transition matrix, but repeating such an operation can compound floating-point inaccuracies; a more effective (and more efficient) approach is to use a matrix power. With this in mind, lets compute the probability of finishing the game in n moves for several values of n, and compare it to the simulation result from above:

In [8]:
def cl_probability(n):
    """Compute the state vector after n turns"""
    mat = cl_markov_matrix()
    v_0 = [1, *np.zeros(100)]
    return np.linalg.matrix_power(mat, n) @ v_0
In [9]:
probs = [cl_probability(i)[-1] for i in range(200)]

plt.hist(sim_games, bins=range(200), normed=True,
         align='mid', alpha=0.5);
plt.plot(np.arange(1, 200), np.diff(probs), color='black')
plt.title('Chutes and Ladders Completion Rate')
plt.xlabel('number of turns')
plt.ylabel('fraction of games completed');

By eye it looks like our direct simulation was quite close to the Markov model, which indicates that we're on the right track in our thinking.

The cumulative distribution is also interesting:

In [10]:
plt.hist(sim_games, bins=range(200), normed=True,
         align='mid', cumulative=True, alpha=0.5);
plt.plot(np.arange(200), probs, color='black');
plt.title('Chutes and Ladders Cumulative Completion Rate')
plt.xlabel('number of turns')
plt.ylabel('cumulative fraction of games completed');

From this, we can see that 90% of single-player games finish within about 72 moves, though it is possible for games to be much longer.

With two players, this translates to a 1% chance that the game will go 72 moves without either of the players winning. Assuming roughly 20 seconds per round, that is about 24 minutes of play time, though from personal experience I can say it feels roughly twenty times that long.

Minimum Game Length

But what if you get really lucky... how quickly can the game end? From our Markov chain results, we see that there is a nonzero probability of finishing the game in just seven moves, and this will happen about once in every 660 games:

In [11]:
probs[1:8]
Out[11]:
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0015110596707818924]

The Markov chain approach doesn't immediately tell us how this happens, though that could be determined with a bit of book-keeping to track which probabilities contribute to the result in question.

Typical Game Length

We've talked about the extremes, but how long should you expect the typical game to last? The answer to that question depends on what you mean by typical. Some common approaches to defining typical are the mean, median, or mode of game lengths.

The mean game length is the expected value of the number of turns, averaged over all possibilities and weighted by their relative probability:

In [12]:
turns = np.arange(1, len(probs))
np.dot(np.diff(probs), turns)
Out[12]:
39.106460290714438

We see that, on average, a player will finish the game in about 39 turns.

This value, though, is weighted by game length, so the small possibility of having a very long game contributes to this estimate in an outsized way. Often the median game length is a more useful statistic:

In [13]:
np.searchsorted(probs, 0.5)
Out[13]:
32

This tells us that approximately 50% of games will be finished in fewer than 32 moves, and 50% will be finished in more than 32 moves.

But both the mean and the median are offset from the peak (or mode) of the probability distribution, which we can compute as follows:

In [14]:
np.argmax(np.diff(probs)) + 1
Out[14]:
22

This says that you'll finish the game in 22 moves more often than any other specific number of moves.

Probably the most useful types of statistics for skewed probability distributions like this one, though, are quantiles. For example:

In [15]:
np.searchsorted(probs, [0.025, 0.975])
Out[15]:
array([ 11, 106])

This tells us that 95% of the time, games will last between 11 and 106 turns. Or, to tighten the constraints a bit:

In [16]:
np.searchsorted(probs, [0.25, 0.75])
Out[16]:
array([22, 50])

half of all games will last between 22 and 50 moves.

Learning More from the Markov Process

Aside from summary statistics, the Markov transition matrix for Chutes and Ladders can tell us a number of other interesting things, which we'll explore here.

Shannon Entropy

One interesting thing we can explore is the entropy of the distribution: you can think of the entropy as a measure of how "spread out" the distribution is in probability space; a high entropy state is one in which we have very little knowledge of where an individual draw will lie.

We can compute and visualize the entropy using the entropy function from SciPy's stats module. Let's plot the entropy as a function of the number of turns played:

In [17]:
from scipy import stats

turns = np.arange(101)
entropy = [stats.entropy(cl_probability(turn)) for turn in turns]

plt.plot(turns, entropy)
plt.xlabel('number of turns')
plt.ylabel('entropy of distribution')
plt.title('Chutes and Ladders Entropy');

We see that the entropy is maximized after eleven turns:

In [18]:
np.argmax(entropy)
Out[18]:
11

What this tells you is that when you complete 11 turns, you should expect positions of various players to be very spread out. Prior to this, the probability distribution is clumped near the beginning of the board. After 11 moves, you can take heart that you're unavoidably on the path to the game's end.

Eigenvectors and Stationary States

The Markov transition matrix can be examined directly to learn about the dynamics of the game. For example, the eigenvectors of a transition matrix with eigenvalue equal to one tell us about stationary states; that is, states which are unchanged after applying the transition matrix.

We can compute these eigenvectors this way:

In [19]:
evals, evecs = np.linalg.eig(mat)
evals[:10]
Out[19]:
array([ 1.00+0.j  ,  0.83+0.j  ,  0.96+0.j  ,  0.40+0.65j,  0.40-0.65j,
        0.76+0.j  ,  0.60+0.32j,  0.60-0.32j,  0.24+0.59j,  0.24-0.59j])

There is only one stationary state in Chutes and Ladders (with eigenvalue equal to one); let's see what it is:

In [20]:
evecs[:, 0]
Out[20]:
array([ 0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j,
        0.+0.j,  0.+0.j,  1.+0.j])

Graciously, the single stationary state is a completed game. What this means is that you are very unlikely to find yourself in an infinitely looping, unending game. I doubt the game designers did this stationary state calculation when designing the game, but I bless them for their accidental success in this regard.

Fundamental Matrix

As we have set up the transition matrix, once you reach state 100 you remain there with 100% probability. In mathematical terms, the end of the game is an "absorbing state", because once you are there you cannot transition to any other state. What this means is that Chutes and Ladders is an absorbing Markov Chain, and we can do a few more interesting things with the transition matrix.

This all starts with computing the fundamental matrix N of the absorbing Markov chain, which is a simple transformation of the Markov transition matrix:

In [21]:
M = cl_markov_matrix(jump_at_end=False)
N = np.linalg.inv(np.eye(100) - M[:100, :100])

(Note that we used the version of our process where jumps happen at the beginning of the turn, so that we can more easily reason about the effects of ladders and chutes).

The (n,m) entry of the fundamental matrix tells us the expected number of times that we will visit state n when starting at state m on our way to the absorbing state.

Let's plot these expected values when starting at state 0, the beginning of the game

In [22]:
expected_visits = N[:, 0]
plt.plot(expected_visits)
plt.xlabel('Square Number')
plt.ylabel('Expected Number of Visits in a Game');

The profile is pretty interesting. The most striking feature is that you're expected to visit square 99 approximately 1.5 times per game, on average! This is due to the fact that once you land on square 99, you have to roll a 1 to exit it; each additional turn that you don't roll a 1 is counted as another visit.

Let's look at the top few other most visited states:

In [23]:
col = N[:, 0]
np.argsort(col)[::-1][:5]
Out[23]:
array([99,  0, 47, 97, 46])

The next most visited state is state 0, which you visit exactly once at the start of each game. After that is square 47, which is visited about 3/4 times per game on average. Interestingly, this square is the top of a chute, and thus is the most used of all chutes in the game.

We can also look at the least visited squares on the board:

In [24]:
np.argsort(col)[:5]
Out[24]:
array([67,  1,  2,  3, 68])

Surprisingly, the least likely square is 67, and is visited even less often than square 1, which you have the possibility of hitting only once per game, if you roll a 1 on the first roll! Looking at the graph above, it is clear that squares 55-75 or so are visited relatively rarely, due to the positions of chutes and of ladders within the game board.

The Big Chute and Big Ladder

One thing my daughter would probably like to know: how often will you to hit the "big slide" or the "big ladder" during a game? These are at squares 87, and 28, respectively:

In [25]:
col[87]
Out[25]:
0.52512870634459763
In [26]:
col[28]
Out[26]:
0.57577654067476058

We see that a player is expected to slide down the big chute a bit over once in every two games on average, and climb the big ladder slightly more often: just under six times in every ten games on average.

Time to Absorption

The fundamental matrix can also be used to compute other quantities; for example the matrix product with a column of ones gives the expected "time to absorption"; that is, the expected number of turns until finishing the game from any given square. Let's use the "jump-at-the-end" version of the matrix again:

In [27]:
np.set_printoptions(precision=1)
np.dot(N.T, np.ones(100))
Out[27]:
array([ 39.7,  35.4,  40.2,  39.7,  37.6,  39.8,  39.6,  39.3,  39.1,
        36.7,  39.2,  38.7,  38.3,  37.9,  37.6,  37.1,  39.6,  35.9,
        35.8,  35.7,  35.6,  34. ,  34. ,  34.5,  34.8,  35.1,  35.3,
        35.4,  22.9,  37.4,  36.9,  36.7,  36.5,  36.3,  36. ,  35.8,
        33.9,  35.7,  35.4,  34.9,  34.5,  34.6,  34. ,  34.7,  33.9,
        32. ,  31.7,  35.3,  30.5,  38.7,  28.9,  20.9,  29.9,  29.5,
        29.1,  28.6,  29.5,  28.6,  28.3,  27.2,  26.3,  25.5,  35.7,
        22.6,  26.3,  20.8,  20.8,  20.9,  20.5,  20.3,  20.3,  16.2,
        21. ,  21. ,  18.1,  19.2,  20. ,  20.6,  20.9,  21. ,   1. ,
        25.9,  24.8,  23.8,  22.9,  21.9,  21.2,  34.8,  17.9,  18.1,
        17.1,  16.2,  16.9,  21. ,  12.3,  19.2,  11. ,  11. ,  20.9,   6. ])

From the first square, the expected time to completion is 39.7 moves: this is half a move than the 39.2 moves we computed previously, because we're using a different form of the matrix (here we count climbing the ladder from square 80 to square 100 as one extra move)

This vector can tell us what the best first roll is: as my daughter has intuited, starting out by rolling a 1 and climbing the first ladder is best: it decreases your expected time to completion by five moves. The worst first roll is a 2: this actually increases your expected time to completion by 0.5 turns, on average: rolling a 2 to start is akin to moving away from the goal.

By the time you get to the 99th square, your espected time to completion is 6 moves, because you have a 1/6 chance of rolling the 1 required to win. Two squares back, on the 97th square, your expected time to completion goes up to 11 moves, because of the chute present on the 98th square.

The fundamental matrix lets you quantititatively explore a number of other properties of the game; for example, we could adjust the matrix to make square 80 an absorber as well, and ask how probable we are to complete the game there versus landing on square 100 directly. Or we could compute quantities such as the variance of the number of visits and number of steps above. For more discussion of applications of the fundamental matrix of an absorbing Markov chain, see the wikipedia article. A nice undergraduate-level introduction to these concepts can be found in Chapter 11 of Grinstead and Snell's “Introduction to Probability” (pdf available here).

Visualizing the Board

Just for fun, let's use Python's tools to visualize these probability states in an intuitive way. We'll use the above scan of the original Chutes and Ladders game as our background, and plot the probability distribution on top of this using a custom colormap with changing transparency

In [28]:
from matplotlib import colors

# Make a blue colorbar with increasing opacity
c = np.zeros((100, 4))
c[:, -1] = np.linspace(0, 1, 100)  # transparency gradient
c[:, 2] = 0.5  # make the map dark blue
TransparencyMap = colors.ListedColormap(c)


def show_board(turn):
    fig, ax = plt.subplots()
    board = plt.imread('ChutesAndLadders-board.gif')
    
    # Compute & reshape the probability vector
    prob = cl_probability(turn)
    prob = prob[1:].reshape(10, 10)[::-1]
    prob[::2] = prob[::2, ::-1]
    
    # Show result over the image of the board
    ax.imshow(board, alpha=0.8)
    im = ax.imshow(prob, extent=[10, 800, 810, 10],
                   norm=colors.LogNorm(vmin=1E-3, vmax=1),
                   cmap=TransparencyMap)
    fig.colorbar(im, ax=ax, label='Fraction of games')
    ax.axis('off')
    ax.set_title(f"Turn {turn}")
    
    return fig

Using this, we can visualize the probability distribution after any number of turns:

In [29]:
show_board(1);

After one turn, the distribution is quite tight: there are six possible states we can find ourselves in.

In [30]:
show_board(2);

After two turns, things become more spread out, though they are still clumped around a few more likely values.

In [31]:
show_board(3);

After 3 turns, the distribution is getting even further apart.

In [32]:
show_board(50);