title: "Lecture 18: Snakes and Ladders" author: "Junaid Hasan" date: "August 6, 2021"¶
Snakes and Ladders¶
The history of the game is interesting, designed to teach children the rewards of good deeds. For more read on Wikipedia
The Game:
Played on a 10 by 10 board.
1 in the lower left corner, 100 on the upper left corner.
Snakes and Ladders at various locations.
{width=40%, height=50%}
Gameplay¶
- Players start off the board, roll a single die and move that many steps ahead.
- If they hit a snake's head, they go down to its tail.
- If they hit a ladder's foot, they go up the ladder to the top.
- Hitting the top of the ladder, or the tail of the snake does nothing.
- The first player to reach square 100 wins.
Why Monte Carlo¶
For a general system if we wish to calculate the probability of some event, we have two approaches:
- Experimenting
- Formal Modeling
In the experimentation, we simply repeat an experiment multiple times and record the relative frequencies of the results.
More likely events happen more frequently.
Less likely events happen less frequently.
The higher samples we run, the higher confidence in our results.
This is the Monte Carlo Method.
Computers are really good at doing repeated experiments.
In [1]:
import random
import statistics
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
success = False
Ladder = {1:38,
4:14,
9:31,
21:42,
28:84,
36:44,
51:67,
71:91,
80:100}
Snake = {98:78,
95:75,
93:73,
87:24,
64:60,
62:19,
56:53,
49:11,
47:26,
16:6
}
MAX_STEPS = 300
Ladder_foot = list(Ladder.keys())
Ladder_top = list(Ladder.values())
Snake_head = list(Snake.keys())
Snake_tail = list(Snake.values())
In [2]:
def playSnakesLadders(Ladder, Snake, MAX_STEPS=300, exact=True):
step = 0
position = 0
record = []
record.append(0)
while(True):
die = random.randint(1,6)
position = position + die
if exact and position > 100 :
position = position - die
step = step + 1
if position in Ladder_foot:
position = Ladder[position]
if position in Snake_head:
position = Snake[position]
record.append(position)
if position >= 100 or step >= MAX_STEPS:
break
return record
In [3]:
playSnakesLadders(Ladder, Snake, exact = False)
Out[3]:
[0, 6, 7, 11, 15, 20, 23, 29, 33, 37, 43, 44, 26, 30, 33, 37, 43, 44, 48, 11, 17, 19, 24, 25, 84, 85, 86, 24, 29, 34, 39, 44, 48, 54, 60, 63, 65, 68, 91, 75, 100]
Questions¶
- Let us run the game $10^5$ times and ask information about
- How long does a typical game last?
In [46]:
NUM_GAMES=100000
gameLength=np.zeros(NUM_GAMES)
for i in range(NUM_GAMES):
gameLength[i] = len(playSnakesLadders(Ladder,Snake, exact=False))-1
In [47]:
histogram = plt.hist(gameLength, bins = range(1,MAX_STEPS+2), density=True)
Therefore the minimum number of steps needed is 7 it seems. Around 0.2 percent of the games are ending in 7 steps.
In [48]:
gameLength.min()
Out[48]:
7.0
In [49]:
gameLength.max()
Out[49]:
274.0
In [50]:
statistics.mode(gameLength)
Out[50]:
22.0
In [51]:
statistics.median(gameLength)
Out[51]:
29.0
In [52]:
statistics.mean(gameLength)
Out[52]:
35.92964
Cumulative Probabilities¶
Suppose we want to know the probability $P(t)$ that a game finishes in $t$ steps or less.
In [11]:
np.histogram?
In [53]:
values, bin_data = np.histogram(gameLength,bins=range(1, MAX_STEPS+2), density=True)
plt.plot(bin_data[:-1], values, c='green')
Out[53]:
[<matplotlib.lines.Line2D at 0x7f84eb8d9190>]
Around 2.7% of the games finish in 20 moves (the mode).
In [54]:
bin_data
Out[54]:
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78,
79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91,
92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104,
105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117,
118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130,
131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156,
157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182,
183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195,
196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208,
209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221,
222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234,
235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247,
248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260,
261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273,
274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286,
287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299,
300, 301])
In [55]:
cumulative = np.cumsum(values)
cumulative
Out[55]:
array([0. , 0. , 0. , 0. , 0. , 0. , 0.0017 ,
0.00789, 0.01734, 0.0301 , 0.04753, 0.06924, 0.09185, 0.11584,
0.14076, 0.16583, 0.19253, 0.21988, 0.24751, 0.27514, 0.30235,
0.33054, 0.35775, 0.3842 , 0.40911, 0.43379, 0.45679, 0.47882,
0.50077, 0.52283, 0.54272, 0.56246, 0.58147, 0.59868, 0.61621,
0.63271, 0.64856, 0.66311, 0.67723, 0.69108, 0.70463, 0.71771,
0.7293 , 0.74057, 0.75188, 0.76255, 0.77261, 0.78178, 0.79109,
0.79985, 0.80806, 0.81554, 0.82353, 0.83081, 0.83831, 0.8454 ,
0.85194, 0.85796, 0.86392, 0.86914, 0.87435, 0.87953, 0.8846 ,
0.88979, 0.8945 , 0.89891, 0.90336, 0.9077 , 0.91128, 0.91509,
0.91883, 0.92222, 0.92553, 0.92867, 0.93203, 0.93519, 0.93786,
0.94056, 0.94302, 0.94558, 0.94784, 0.94996, 0.95236, 0.95428,
0.95643, 0.9584 , 0.96012, 0.96182, 0.96344, 0.96495, 0.96647,
0.96788, 0.96949, 0.97083, 0.97224, 0.97335, 0.97456, 0.97559,
0.9766 , 0.9774 , 0.9785 , 0.97937, 0.98036, 0.98126, 0.98212,
0.9829 , 0.98362, 0.98428, 0.98512, 0.98563, 0.98626, 0.98691,
0.98746, 0.98808, 0.9886 , 0.9891 , 0.98962, 0.99016, 0.99056,
0.99091, 0.99142, 0.99194, 0.99234, 0.99272, 0.99304, 0.99344,
0.99369, 0.99399, 0.9942 , 0.99454, 0.9948 , 0.99504, 0.99525,
0.99541, 0.99554, 0.99567, 0.9958 , 0.996 , 0.99618, 0.99634,
0.99646, 0.99666, 0.99683, 0.99697, 0.99709, 0.99724, 0.9974 ,
0.99748, 0.99756, 0.99762, 0.99769, 0.99778, 0.99792, 0.99803,
0.9981 , 0.99816, 0.99825, 0.99832, 0.99845, 0.99854, 0.99859,
0.99863, 0.99868, 0.99869, 0.99877, 0.99881, 0.99889, 0.99894,
0.99903, 0.9991 , 0.99912, 0.99913, 0.99914, 0.99921, 0.99927,
0.99933, 0.99935, 0.9994 , 0.99942, 0.99944, 0.99945, 0.99946,
0.99947, 0.99951, 0.99952, 0.99953, 0.99954, 0.99957, 0.99957,
0.9996 , 0.99961, 0.99964, 0.99966, 0.99966, 0.99967, 0.99968,
0.99971, 0.99972, 0.99974, 0.99975, 0.99975, 0.99975, 0.99975,
0.99975, 0.99975, 0.99976, 0.99977, 0.99979, 0.99979, 0.99979,
0.99981, 0.99982, 0.99983, 0.99985, 0.99985, 0.99986, 0.99987,
0.99987, 0.99987, 0.99987, 0.99987, 0.99988, 0.99988, 0.99988,
0.99988, 0.99988, 0.9999 , 0.9999 , 0.9999 , 0.9999 , 0.9999 ,
0.9999 , 0.99992, 0.99992, 0.99995, 0.99995, 0.99995, 0.99995,
0.99995, 0.99995, 0.99995, 0.99995, 0.99995, 0.99995, 0.99995,
0.99995, 0.99995, 0.99995, 0.99995, 0.99996, 0.99996, 0.99996,
0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996,
0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99996, 0.99997,
0.99997, 0.99997, 0.99999, 0.99999, 0.99999, 0.99999, 0.99999,
1. , 1. , 1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. , 1. ])
In [15]:
plt.plot(bin_data[:-1], cumulative, c='blue')
Out[15]:
[<matplotlib.lines.Line2D at 0x7f84eb2fc100>]
Therefore half of the games finish in 31 steps or less. This can also be observed because the median is 31.
- Question: What percent of games take 100 moves or more?
- Question: What percent of games take 50 moves or more?
In [16]:
ladderCounter = {}
snakeCounter = {}
def resetCounters():
for l in Ladder_foot:
ladderCounter[l] = 0
for s in Snake_head:
snakeCounter[s] = 0
In [17]:
resetCounters()
In [18]:
ladderCounter
Out[18]:
{1: 0, 4: 0, 9: 0, 21: 0, 28: 0, 36: 0, 51: 0, 71: 0, 80: 0}
In [19]:
snakeCounter
Out[19]:
{98: 0, 95: 0, 93: 0, 87: 0, 64: 0, 62: 0, 56: 0, 49: 0, 47: 0, 16: 0}
In [20]:
def playSnakesLadders2(Ladder, Snake, MAX_STEPS=300, exact=True):
step = 0
position = 0
record = []
record.append(0)
while(True):
die = random.randint(1,6)
position = position + die
if exact and position > 100 :
position = position - die
step = step + 1
if position in Ladder_foot:
ladderCounter[position] += 1
position = Ladder[position]
if position in Snake_head:
snakeCounter[position] += 1
position = Snake[position]
record.append(position)
if position >= 100 or step >= MAX_STEPS:
break
return record
In [21]:
NUM_GAMES = 10000
resetCounters()
for i in range(NUM_GAMES):
playSnakesLadders2(Ladder,Snake, exact=False)
In [22]:
1/6
Out[22]:
0.16666666666666666
In [23]:
(1/6) + (2/36) + (1/216)
Out[23]:
0.22685185185185183
In [24]:
16000/100000
Out[24]:
0.16
In [25]:
plt.bar(list(ladderCounter.keys()), list(ladderCounter.values()))
Out[25]:
<BarContainer object of 9 artists>
In [26]:
plt.bar(list(snakeCounter.keys()), list(snakeCounter.values()))
Out[26]:
<BarContainer object of 10 artists>
Answer:¶
- It seems that the ladder at 28 and 36 are used the most.
- The snake at 48 is used the most.
More questions¶
- Why is the ladder at position 1 used the least?
- Why is ladder at position 4 used the second least?
- What about other ladders?
Answer:¶
- A probable answer is that the ladder at position 1 can only be reached if the throw is a 1.
- Let us see how many fraction of games the ladder at position 1 was used.
In [27]:
ladderCounter[1]/NUM_GAMES
Out[27]:
0.1676
- On the first throw, what is the chance of throwing a 1?
- Answer: $\frac{1}{6}$.
- Therefore we expect with $\frac{1}{6}$ probability ladder at position 1 is used. Indeed
In [28]:
1/6.0
Out[28]:
0.16666666666666666
Limitations of Monte Carlo approach¶
Monte Carlo works well, but sometimes it is simply not possible to repeat an experiment multiple times. Sometimes we only get one shot.
Moreover, random numbers we use may not be perfectly random. There is an whole area devoted to generating random numbers in the best possible way.
How?¶
- When a player is at a particular square, we do not care how they arrived at that square.
- The only thing that matters is what they get when they roll the die.
- In other words, the past does not influence what will happen at the next stage.
Markov Chain Model¶
- If we model the game as a Markov Chain.
- What are the states?
- What is a transition matrix?
Answer:¶
- Suppose each box represents a state.
- From one box $n$ we can go to 6 boxes ahead $n+1, n+2, \ldots, n+6$ with probability $\frac{1}{6}$ each.
- However, we must encorporate snakes and ladders.
- If there is a snake head or a ladder foot at state $k$ then we get instantly teleported from $k$ to another box $l$.
- In other words instead of the $(n,k)$ entry as $\frac{1}{6}$ we must write the $(n,l)$ entry as $\frac{1}{6}$.
Caveats¶
- There are two minor caveats:
- Consider the state $52$.
{width=50%, height=5%}- There is a snake that moves from 56 to 53.
- Therefore from 52 the probability of reaching 53 is $\frac{2}{6}$.
- Consider the final row
{width=50%, height=5%}- If we are playing the inexact game, then from state 97 the probability looks like
- $p(78 | 97) = \frac{1}{6}$ because of the snake at 98.
- $p(99 | 97) = \frac{1}{6}$.
- $p(100| 97) = \frac{4}{6}$ because throwing 3,4,5,6 all result in win.
- However, if we play exact game, then from state 97
- $p(78 | 97) = \frac{1}{6}$ because of snake at 98.
- $p(99 | 97) = \frac{1}{6}$.
- $p(100| 97) = \frac{1}{6}$.
- $p(97 | 97) = \frac{3}{6}$.
Questions¶
What should be the size of the transition matrix?
A guess is 101 since we start outside the board state 0.
However, we can actually use 82 states only. $$101 - \text{ num of ladders} - \text{ num of snakes}$$
This is because the states at the head of a snake, or a foot of a ladder are not necessary since we immediately teleport from those states.
82 vs 101 can be better because matrix multiplication is $O(n^3)$ and $\left ( \frac{82}{100} \right)^3 \approx 0.54$
Let us use numpy to create the matrices. I use numpy because of performance reasons. Sympy will give exact computation but is slow.
In [29]:
slmatrix=np.zeros(shape = (82,82))
In [30]:
slmatrix
Out[30]:
array([[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.]])
In [31]:
box2state = {}
countSL = 0
for i in range(0, 101):
if i in Ladder_foot or i in Snake_head:
countSL +=1
box2state[i] = i - countSL
for i in range(0,101):
if i in Ladder_foot:
box2state[i] = box2state[Ladder[i]]
if i in Snake_head:
box2state[i] = box2state[Snake[i]]
The first ladder is between box 1 and 38. Moreover, box 38 is state 31 (there are 7 snakes or ladders before box 38). Therefore we want box2state[1] = 31
In [32]:
box2state
Out[32]:
{0: 0,
1: 31,
2: 1,
3: 2,
4: 11,
5: 3,
6: 4,
7: 5,
8: 6,
9: 25,
10: 7,
11: 8,
12: 9,
13: 10,
14: 11,
15: 12,
16: 4,
17: 13,
18: 14,
19: 15,
20: 16,
21: 35,
22: 17,
23: 18,
24: 19,
25: 20,
26: 21,
27: 22,
28: 69,
29: 23,
30: 24,
31: 25,
32: 26,
33: 27,
34: 28,
35: 29,
36: 37,
37: 30,
38: 31,
39: 32,
40: 33,
41: 34,
42: 35,
43: 36,
44: 37,
45: 38,
46: 39,
47: 21,
48: 40,
49: 8,
50: 41,
51: 54,
52: 42,
53: 43,
54: 44,
55: 45,
56: 43,
57: 46,
58: 47,
59: 48,
60: 49,
61: 50,
62: 15,
63: 51,
64: 49,
65: 52,
66: 53,
67: 54,
68: 55,
69: 56,
70: 57,
71: 75,
72: 58,
73: 59,
74: 60,
75: 61,
76: 62,
77: 63,
78: 64,
79: 65,
80: 81,
81: 66,
82: 67,
83: 68,
84: 69,
85: 70,
86: 71,
87: 19,
88: 72,
89: 73,
90: 74,
91: 75,
92: 76,
93: 59,
94: 77,
95: 61,
96: 78,
97: 79,
98: 64,
99: 80,
100: 81}
Create the transition matrix¶
- Let us now create the transition matrix.
In [61]:
exact = True
slmatrix=np.zeros(shape = (82,82))
for i in range(0,101):
if i in Ladder_foot or i in Snake_head:
continue
if i <=94 :
for j in range(1,7):
slmatrix[box2state[i],box2state[i+j]] += 1/6.0
elif i>94 and exact==False:
for j in range(1,7):
if i+j > 100:
slmatrix[box2state[i],box2state[100]] += 1/6.0
else:
slmatrix[box2state[i],box2state[i+j]] += 1/6.0
else:
for j in range(1,7):
if i+j > 100:
slmatrix[box2state[i],box2state[i]] += 1/6.0
else:
slmatrix[box2state[i],box2state[i+j]] += 1/6.0
In [62]:
slmatrix[80]
Out[62]:
array([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. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0.83333333, 0.16666667])
In [63]:
gamestart = np.zeros(82)
gamestart[0] = 1
gamestart
Out[63]:
array([1., 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., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
The prob vector stores the probability of being in a state after moves.
In [64]:
moves = 7
P_n = np.linalg.matrix_power(slmatrix,moves)
prob = np.matmul(gamestart, P_n)
prob
Out[64]:
array([0. , 0. , 0. , 0. , 0.01974023,
0.00477252, 0.00609425, 0.00769819, 0.02860654, 0.0145176 ,
0.0130637 , 0.01448903, 0.01603581, 0.01781121, 0.01730038,
0.02263017, 0.02106553, 0.01972594, 0.02436628, 0.0306963 ,
0.02628458, 0.05069373, 0.02715621, 0.02614883, 0.02430556,
0.02965321, 0.02253372, 0.01855781, 0.01769333, 0.02029392,
0.01654664, 0.01402463, 0.0127672 , 0.0119563 , 0.01133473,
0.03378987, 0.0179291 , 0.03833376, 0.0223944 , 0.02382687,
0.0233089 , 0.01593579, 0.00887346, 0.01617513, 0.00959505,
0.00797325, 0.00712663, 0.00816972, 0.00785894, 0.00920568,
0.00506544, 0.00405093, 0.00208976, 0.001261 , 0.01286723,
0.00212191, 0.002497 , 0.00273277, 0.00315786, 0.00717664,
0.00265775, 0.00558699, 0.00251843, 0.00219693, 0.00303284,
0.00187543, 0.00082519, 0.00042867, 0.00026792, 0.03094993,
0.0042081 , 0.00474037, 0.00539052, 0.00610497, 0.00688372,
0.00678727, 0.00384374, 0.00350437, 0.00249343, 0.00169324,
0.00041795, 0.00151106])
Printing the board¶
However to print the board, we must convert the states back to squares. In particular 82 states will have to be converted back to 100 squares (plus the 0 box)
This converted probability is the newprob
In [65]:
newprob = np.zeros((10,10))
newprob
Out[65]:
array([[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., 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.]])
In [66]:
newprob = np.zeros((10,10))
countS = 0
countL = 0
for i in range(0,10):
for j in range(1,11):
box = 10*i +j
if box in Ladder_foot:
countL +=1
continue
if box in Snake_head:
countS +=1
continue
state = box - countL - countS
newprob[i,j-1] = prob[state]
newprob
Out[66]:
array([[0. , 0. , 0. , 0. , 0. ,
0.01974023, 0.00477252, 0.00609425, 0. , 0.00769819],
[0.02860654, 0.0145176 , 0.0130637 , 0.01448903, 0.01603581,
0. , 0.01781121, 0.01730038, 0.02263017, 0.02106553],
[0. , 0.01972594, 0.02436628, 0.0306963 , 0.02628458,
0.05069373, 0.02715621, 0. , 0.02614883, 0.02430556],
[0.02965321, 0.02253372, 0.01855781, 0.01769333, 0.02029392,
0. , 0.01654664, 0.01402463, 0.0127672 , 0.0119563 ],
[0.01133473, 0.03378987, 0.0179291 , 0.03833376, 0.0223944 ,
0.02382687, 0. , 0.0233089 , 0. , 0.01593579],
[0. , 0.00887346, 0.01617513, 0.00959505, 0.00797325,
0. , 0.00712663, 0.00816972, 0.00785894, 0.00920568],
[0.00506544, 0. , 0.00405093, 0. , 0.00208976,
0.001261 , 0.01286723, 0.00212191, 0.002497 , 0.00273277],
[0. , 0.00315786, 0.00717664, 0.00265775, 0.00558699,
0.00251843, 0.00219693, 0.00303284, 0.00187543, 0. ],
[0.00082519, 0.00042867, 0.00026792, 0.03094993, 0.0042081 ,
0.00474037, 0. , 0.00539052, 0.00610497, 0.00688372],
[0.00678727, 0.00384374, 0. , 0.00350437, 0. ,
0.00249343, 0.00169324, 0. , 0.00041795, 0.00151106]])
The code below generates newprob from prob. In other words is generates the box data from the state data.
In [67]:
#generate box data from state data
newprob = np.zeros((10,10))
newbox = np.zeros((10,10))
countS = 0
countL = 0
for i in range(0,10):
for j in range(1,11):
box = 10*i +j
newbox[i,j-1] = box
if box in Ladder_foot:
countL +=1
continue
if box in Snake_head:
countS +=1
continue
state = box - countL - countS
newprob[i,j-1] = prob[state]
In [68]:
newbox
Out[68]:
array([[ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.],
[ 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.],
[ 21., 22., 23., 24., 25., 26., 27., 28., 29., 30.],
[ 31., 32., 33., 34., 35., 36., 37., 38., 39., 40.],
[ 41., 42., 43., 44., 45., 46., 47., 48., 49., 50.],
[ 51., 52., 53., 54., 55., 56., 57., 58., 59., 60.],
[ 61., 62., 63., 64., 65., 66., 67., 68., 69., 70.],
[ 71., 72., 73., 74., 75., 76., 77., 78., 79., 80.],
[ 81., 82., 83., 84., 85., 86., 87., 88., 89., 90.],
[ 91., 92., 93., 94., 95., 96., 97., 98., 99., 100.]])
Heatmap¶
Seaborn package allows us to convert the box probabilities newprob into an heatmap.
In [69]:
import seaborn as sns
from matplotlib import colors
for i in range(10):
if i%2 == 1:
newprob[i,:] = np.flip(newprob[i,:])
newbox[i,:] = np.flip(newbox[i,:])
plot = sns.heatmap(newprob,annot=newbox,norm=colors.LogNorm(vmin=1E-3, vmax=1))
plot.invert_yaxis()
plot.axis('off')
fig = plot.get_figure()
fig.savefig(f"img-{moves}.png")
Another way to print the board¶
Suppose we have a 4 by 4 board for simplicity.
Note that additionally we must flip every even row.
In [70]:
a = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
a
Out[70]:
array([[ 1, 2, 3, 4],
[ 5, 6, 7, 8],
[ 9, 10, 11, 12],
[13, 14, 15, 16]])
The code below flips the even rows and it invertes the rows so that the rows start from the bottom.
In [71]:
for i in range(4):
if i%2 == 1:
a[i,:] = np.flip(a[i,:])
for i in range(4//2):
tmp = np.copy(a[3-i,:])
a[3-i,:] = a[i,:]
a[i,:] = tmp
a
Out[71]:
array([[16, 15, 14, 13],
[ 9, 10, 11, 12],
[ 8, 7, 6, 5],
[ 1, 2, 3, 4]])
The code below inputs a picture containing the board and adds it to the plot as a background with transparency.
- The first section is about generating the box data from the state data.
- The second section is about flipping the even rows and starting the rows from the bottom.
- The third section is about plotting.
In [73]:
# 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)
#generate box data from state data
newprob = np.zeros((10,10))
countS = 0
countL = 0
for i in range(0,10):
for j in range(1,11):
box = 10*i +j
if box in Ladder_foot:
countL +=1
continue
if box in Snake_head:
countS +=1
continue
state = box - countL - countS
newprob[i,j-1] = prob[state]
# print board
fig, ax = plt.subplots()
board = plt.imread('snakes-ladders2.gif')
# Compute & reshape the probability vector
for i in range(10):
if i%2 == 1:
newprob[i,:] = np.flip(newprob[i,:])
for i in range(10//2):
tmp = np.copy(newprob[9-i,:])
newprob[9-i,:] = newprob[i,:]
newprob[i,:] = tmp
# Show result over the image of the board
ax.imshow(board, alpha=0.9)
im = ax.imshow(newprob, 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 {moves}")
Out[73]:
Text(0.5, 1.0, 'Turn 7')
Exact computation¶
- Before this we have been using
numpy. Let us attempt to use sympy to do some exact computation.
In [74]:
import sympy as sp
Since most entries in each row are zeros, the slmatrix2 is a sparse matrix (a more efficient way to store matrix data).
In [75]:
exact = False
slmatrix3=sp.zeros(82)
slmatrix2 = sp.SparseMatrix(slmatrix3)
for i in range(0,101):
if i in Ladder_foot or i in Snake_head:
continue
if i <=94:
for j in range(1,7):
slmatrix2[box2state[i],box2state[i+j]] += sp.Rational(1,6)
elif i>94 and exact==False:
for j in range(1,7):
if i+j > 100:
slmatrix2[box2state[i],box2state[100]] += sp.Rational(1,6)
else:
slmatrix2[box2state[i],box2state[i+j]] += sp.Rational(1,6)
else:
for j in range(1,7):
if i+j > 100:
slmatrix2[box2state[i],box2state[i]] += sp.Rational(1,6)
else:
slmatrix2[box2state[i],box2state[i+j]] += sp.Rational(1,6)
Caution¶
Run the cell below with caution. The print command for sympy tries to print an exact fraction using mathjax (latex on the browser). It can be slow.
In [76]:
gamestart2 = sp.zeros(1,82)
gamestart2[0] = 1
gamestart2
Out[76]:
$\displaystyle \left[\begin{array}{cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc}1 & 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 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$
Let us evaluate the exact probabilities and store in prob2 after moves
In [79]:
moves = 7
P_n2 = slmatrix2**moves
prob2 = gamestart2*P_n2
Caution¶
Again run the cell below with caution. The print can be super slow.
In [80]:
prob2
Out[80]:
$\displaystyle \left[\begin{array}{cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc}0 & 0 & 0 & 0 & \frac{307}{15552} & \frac{167}{34992} & \frac{853}{139968} & \frac{2155}{279936} & \frac{1001}{34992} & \frac{127}{8748} & \frac{1219}{93312} & \frac{169}{11664} & \frac{4489}{279936} & \frac{277}{15552} & \frac{4843}{279936} & \frac{6335}{279936} & \frac{5897}{279936} & \frac{2761}{139968} & \frac{6821}{279936} & \frac{8593}{279936} & \frac{3679}{139968} & \frac{14191}{279936} & \frac{1267}{46656} & \frac{305}{11664} & \frac{7}{288} & \frac{2767}{93312} & \frac{1577}{69984} & \frac{5195}{279936} & \frac{1651}{93312} & \frac{5681}{279936} & \frac{193}{11664} & \frac{1963}{139968} & \frac{1787}{139968} & \frac{3347}{279936} & \frac{3173}{279936} & \frac{1051}{31104} & \frac{1673}{93312} & \frac{3577}{93312} & \frac{6269}{279936} & \frac{3335}{139968} & \frac{725}{31104} & \frac{1487}{93312} & \frac{23}{2592} & \frac{283}{17496} & \frac{1343}{139968} & \frac{31}{3888} & \frac{665}{93312} & \frac{2287}{279936} & \frac{275}{34992} & \frac{859}{93312} & \frac{709}{139968} & \frac{7}{1728} & \frac{65}{31104} & \frac{353}{279936} & \frac{1801}{139968} & \frac{11}{5184} & \frac{233}{93312} & \frac{85}{31104} & \frac{221}{69984} & \frac{2009}{279936} & \frac{31}{11664} & \frac{391}{69984} & \frac{235}{93312} & \frac{205}{93312} & \frac{283}{93312} & \frac{175}{93312} & \frac{77}{93312} & \frac{5}{11664} & \frac{25}{93312} & \frac{361}{11664} & \frac{589}{139968} & \frac{1327}{279936} & \frac{503}{93312} & \frac{1709}{279936} & \frac{1927}{279936} & \frac{475}{69984} & \frac{269}{69984} & \frac{109}{31104} & \frac{79}{34992} & \frac{143}{93312} & \frac{13}{31104} & \frac{89}{46656}\end{array}\right]$
Computing stationary state exactly.¶
Let us see if we can calculate the eigenvectors exactly.
In [81]:
transpose = slmatrix2.T
**THE CELL BELOW CAN BE SUPER SLOW*** It tries to solve a 82 degree polynomial exactly.
Computing eigenvectors exactly is very hard. Therefore we will stop the cell below after a while.
In [82]:
transpose.eigenvects?
Approx computation¶
- Let us attempt to find eigenvectors as an approximate. Therefore we use
slmatrixwhich stores an approximate data, instead ofslmatrix2which had the exact data.
In [83]:
transpose_approx = slmatrix.T
In [84]:
eigenvalues, eigenvectors = np.linalg.eig(transpose_approx)
Approx eigenvectors, but they are inaccurate. Eigenvectors are highly sensitive to change in data, so numpy does not help. Sympy is super slow, so we do not have a way of computing eigenvectors.
Let us print a few eigenvectors.
In [85]:
eigenvectors[0]
Out[85]:
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])
We can observe that numpy cannot give us usable eigenvectors. And sympy takes too much time, because it is an 82 by 82 matrix.
However, we know that the stationary state is $[0,0, \ldots, 0,1]$ because in the long run all games end at the finish box 100.
In [89]:
np.linalg.matrix_power(slmatrix, 300)[:,-1]
Out[89]:
array([0.99999164, 0.99999148, 0.99999166, 0.99999163, 0.9999917 ,
0.99999177, 0.99999186, 0.9999918 , 0.99999196, 0.9999921 ,
0.99999222, 0.99999233, 0.9999925 , 0.99999286, 0.9999929 ,
0.99999293, 0.99999294, 0.99999338, 0.99999324, 0.99999315,
0.99999309, 0.99999304, 0.999993 , 0.99999249, 0.99999264,
0.99999269, 0.99999275, 0.99999281, 0.99999289, 0.99999295,
0.99999299, 0.99999307, 0.99999321, 0.99999334, 0.99999329,
0.99999347, 0.99999325, 0.99999351, 0.999994 , 0.99999411,
0.99999451, 0.99999502, 0.99999476, 0.99999485, 0.99999495,
0.99999507, 0.999995 , 0.99999509, 0.99999535, 0.99999558,
0.99999576, 0.99999651, 0.99999689, 0.99999686, 0.99999683,
0.99999688, 0.99999689, 0.99999687, 0.99999667, 0.99999666,
0.99999716, 0.99999693, 0.99999677, 0.99999668, 0.99999663,
0.99999662, 0.99999562, 0.99999589, 0.99999612, 0.99999634,
0.99999655, 0.99999671, 0.99999748, 0.9999974 , 0.99999761,
0.99999778, 0.99999761, 0.99999846, 0.99999878, 0.99999878,
1. , 1. ])
This shows that at a high enough power $n$ the last column of the matrix $P^n$ is $$\begin{bmatrix}1 \\ 1 \\ \vdots \\ 1 \end{bmatrix}.$$
Therefore the stationary distribution is $[0, 0, \ldots, 1]$ where you reach the end goal eventually.
In [ ]: