- Mon 16 November 2020
- misc
- Akindele Aroge
This notebook is a detailed walkthrough of the chutes and ladders simulation which was first presented by Jake here. In this notebook, I have further included additional interesting questions.
![[img: Chutes and ladders game board]](http://jakevdp.github.io/images/ChutesAndLadders-board.gif)
Simulation¶
One way to quick gain some insight into the stochastic characeristics of the game is to simply "play" several games and analyse the statistic of interest. The way to "play" seveal games here ill be to leverage computational power to simulate several games.
SPECIAL_TRANSITIONS = {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}
import numpy as np
from random import Random
def play_game(rSeed=None, die_sides=6):
"""Simulate a complete chutes and ladders game and
return the number of turns to completion
"""
#rand = np.random.RandomState(rSeed)
rand = Random(rSeed)
position = 0
turns = 0
while position < 100:
turns += 1
roll = rand.randint(1, die_sides)
if (position + roll) <= 100:
position += roll
position = SPECIAL_TRANSITIONS.get(position, position)
return turns
The above function returns the number of turns for one instance of chutes and ladders
play_game()
We would now simulate a decent number of games to give us a sense of the distribution of the number of turns it requires to complete a game
import matplotlib.pyplot as plt
%matplotlib inline
n_games = 10000
n_turns = [play_game() for game in range(n_games)]
rv_max = np.max(n_turns)
plt.hist(n_turns, bins=range(300), density=True);
plt.xlabel('Number of Turns to win')
plt.ylabel('Fraction of games')
plt.title('Game lengths of Simulated Schutes and Ladders');
While this is interesting, the above approach only approximates the actual distribution of game lengths with a 100000 samples.
Interestingly, it is possible to derive the exact distribution with probability theory and harnessing the memoryless property of the game!
The Game as a Markov Process¶
According to the wikipedia article here:
"A stochastic process has the Markov property if the conditional probability distribution of future states of the process (conditional on both past and present states) depends only upon the present state, not on the sequence of events that preceded it. A process with this property is called a Markov process."
Indeed, the conditionl distribution of the next state in the schutes and ladder game is only dependent on the current position. The implication is that the entire game may be modelled as a Markov process - and that is what this next section does
Let's take a concrete example -At each position, we could roll a die which will open up say 6 (is we have a six-sided die) options. Also, the probabilities for each of this option will typically be uniform if we have a fair die (We will explore what happens with when the die is biased). Concretely, assuming we are currently on square 2, the possibilities are 3, 4, 5, 6, 7, 8.
Mathematically: from a position x, we can advance to some set of new positions y with probabilities that sum to 1
$$ \sum_{i=1}^n P(y_i|x) = 1 $$For each state, we could store this possibilities in a vector. E.g from state 2 ($x = 2$), we will have:
2: [0, 0, 0, 1/6, 1/6, 1/6, 1/6, 1/6, 1/6, ...]
Note that this vector completely captures to stochastic bechaviour of the game if you are on position 2 in this case. We can do the same for every other box to give us a so called transition matrix!.
In a similar manner, the special chutes and ladder transitions within the game will also be captured in a transition matrix. In this case, the vectors will hold a 1 for transitions where there is a chute or a ladder.
def trans_mat(die_sides=6):
# create the basic transition matrix:
mat = np.zeros((101, 101))
for i in range(101):
mat[i+1:i+1+die_sides, i] = 1.0 / die_sides
# fix rolls off the board:
# typically, if the roll leads off the board, we don't move
# therefore Assign the rest of the probability to not moving
mat[range(101), range(101)] = 1.0 - mat.sum(0)
# create another trans matrix for schutes and ladders
cl_mat = np.zeros((101, 101))
idx = [SPECIAL_TRANSITIONS.get(i, i) for i in range(101)]
cl_mat[idx, range(101)] = 1
return np.matmul(cl_mat, mat)
mat = trans_mat()
# plt.matshow(mat)
# plt.colorbar()
With the transition matrix, playing the game now involves evolving the state vector which starts as:
v_0 = [1, 0, 0, 0, ...] (there's a probability of 1 we are at position 0 when we start). The state vector is evolved my multiplying by the transition matrix.
def get_game_state(n):
"""computet the state vector after n turns"""
mat = trans_mat()
v_0 = [1, *np.zeros(100)]
return np.linalg.matrix_power(mat, n) @ v_0
We can now return to generating the actual distribbution: For that, we will be interested in the last value in the state vector which is the probability of being at the finishing spot!.
probs = [get_game_state(i)[-1] for i in range(rv_max)]
plt.plot(np.arange(1, rv_max), np.diff(probs), label='prob distribution')
plt.hist(n_turns, bins=range(rv_max), density=True, label='sample distribution');
plt.xlabel('Number of Turns to win')
plt.ylabel('Fraction of games')
plt.title('Game lengths of Simulated Schutes and Ladders');
plt.legend()
Indeed the probability distribution looks fairly similar to the sample ditribution that was estimated earlier
Game Lengths¶
Another way of representing the distribution is the cumulative density function. From this, we could directly read the probability that a game last for as long as n turns or less
plt.plot(np.arange(rv_max), probs, color='black')
plt.hist(n_turns, bins=range(rv_max), density=True, align='mid', cumulative=True)
plt.axhline(y=0.9, color='r', alpha=0.5)
plt.axvline(x=71.5, color='r', alpha=0.5)
plt.xlabel('Number of turns')
plt.ylabel('Cumulative fraction of games completed')
plt.title('Cumulative distribution of number of turn to game completion')
plt.grid()
plt.xlim(0, rv_max)
plt.ylim(0, 1)
plt.grid(False)
extraticks = [72]
plt.xticks(list(plt.xticks()[0]) + extraticks);
plt.yticks(list(plt.yticks()[0]) + [0.9]);
For two players (as is typically the case) we can compute this by reasoning probabilistically from the simplge player case as follows:
$$P(T\le 72) = 0.9$$for two players with turns $T_1$ and $T_2$ (and independent player assumption), the probability may be computed as : $$P(min(T_1, T_2) \le 72) \\ = 1 - P(min(T_1, T_2) > 72) \\ = 1 - [P(T_1 > 72) \cap P(T_2 > 72)]\\ $$
1 - (0.1 * 0.1)
That is 99% ends by 72 moves. As expected, this probability is higher since more player simply implies more chances
Quickest game¶
We will explore the game lengths further. It may be interesting to know the least amount of time the game will take. Within our workflow, that translates to the first time the probability of being in state 100 turns out a non-zero value as we increase the number of turns. As shown in the line below, the minimum number of turns to end the game is 7.
np.nonzero(probs)[0][0]
But what constitutes the 7 turns?. This question is effective asking the shotest path to 100 starting from point 0. This may be computed as follows
from scipy.sparse.csgraph import shortest_path
lengths, predecessors = shortest_path(mat.T, indices=0, directed=True,
unweighted=True, return_predecessors=True)
lengths
The length variable above hold the minimum number of turns required to get to the specific position. For example, from point 0, it takes 0 turn(s) to get to point 0, while it takes inf (infinite) turns to get to point 1 (this simply means it is not possible to reach point 1 from point zero (See the game image to verify that).
From the information encoded in lengths, it is possible to extract the relevant 7 turns as follows;
path = [100]
while path[0] > 0:
path.insert(0, predecessors[path[0]])
path
There we have it: the lucky game is what the above represents - roll a 1, roll 1, roll 6, and so on.
Typical Game¶
Typical game may be considered as the mean, median or mode depending on the relevant interpretation. Here we calculate them all.
First the mean: we simply use the expectation formula:
$$\sum_{turns = 1}^{maximum \ turns} turns * prob_{turns}$$np.dot(np.diff(probs), np.arange(1, len(probs)))
So on average, games will end after 32 turns.
The median will be more in tune to our sense of a typical game since it is not impacted by the rare instances of very long games. For that we get a value of 32 turns:
np.searchsorted(probs, 0.5)
The most common game length we experience as we play a large number of games will be given by the mode:
np.argmax(np.diff(probs)) + 1 # to account for zero indexing: 0 is 1 turn
More Information to be extracted¶
Entropy¶
Recall the state vector discussed previously: as the game evolves, the jprobability values for the states change. But what how does this cahnge happen? In particular, how are the probabilities distrited across the the stated? How spreadout or how concentrated is the distribution.
These type of questions are formalized in the entropy measure.
from scipy import stats
turns = np.arange(101)
entropy = [stats.entropy(get_game_state(turn)) for turn in turns]
plt.plot(turns, entropy)
plt.xlabel('number of turns')
plt.ylabel('entropy of distribution')
plt.title('chutes and ladders');
np.argmax(entropy)
The above graph shows the trend of the entroy as the game proceeds. That entropy is maximized after 11 turns implies that that the game is least predictable after 11 turns. As the number of turns increase however, the entropy begins to decrease; implying that there's not many more states we can find ourselves on the board.
Eigenvectors and Stationary States¶
We have been woring with a (transition) matrix and a (state) vector. For a matrix and vector manipulation like this, we may be interested in knowing how the matrix transforms the vector. One widely studied aspect of transformations are eigen vectors and their eigen values.
For the transition matrix $T$, the eigen state ($s$) is that which is only stretched or contracted by the transition matrix as follows:
$$Ts = \lambda s$$where \lambda is just a scalar value.
In the special case where lambda is 1, this implies that the state s does not change and this is what is typically described as an absorbing state. For the game as we know it, we expect this only at when we are on the box 100 and we can check this state.
evals, evecs = np.linalg.eig(mat)
evals[:10]
evecs[:, 0]
Indeed we see that the wigen vector with the eigen value of 1 the vector representing the state where we are on the state with probability of 1. This is good news as other state with eigen value 1 will mean the game could get stuck in some cases
Conclusion¶
It's really fascinating how much could be revealed in an otherwise simple game using powerful mathematical ideas. Of course it's possible to glean more information still and I hope to include more sections as the ideas pop. Also, this exploration follows the presentation by Jake VanderPlas closely but still a great execise to work through.