This short notebook aims at simulating trials of the so-called Monty-Hall problem, and thus helping to convince about the result thanks to a numerical results rather than a possibly-unclear proof.
There is $M \geq 3$ doors, and behind only one of them there is a treasure. The goal of the player is to find the treasure, following this game:
Finally, the chosen door is opened, and the player wins this round if she found the treasure.
The goal of this notebook is to numerically prove that the choice of always switching to the last door is the best one.
In [1]:
import random
In [2]:
M = 3
allocation = [False] * (M - 1) + [True] # Only 1 treasure!
assert set(allocation) == {True, False} # Check: only True and False
assert sum(allocation) == 1 # Check: only 1 treasure!
Just to check:
In [3]:
allocation
Out[3]:
We can generate a random spot for the treasure by simply shuffling (with random.shuffle()):
In [4]:
def randomAllocation():
r = allocation[:]
random.shuffle(r)
return r
Let us quickly check this function randomAllocation():
In [5]:
for _ in range(10):
print(randomAllocation())
We need to write a small function to simulate the choice of the door to show to the player, show():
In [6]:
def last(r, i):
# Select a random index corresponding of the door we keep
if r[i]: # She found the treasure, returning a random last door
return random.choice([j for (j, v) in enumerate(r) if j != i])
else: # She didn't find the treasure, returning the treasure door
# Indeed, the game only removes door that don't contain the treasure
return random.choice([j for (j, v) in enumerate(r) if j != i and v])
In [7]:
for _ in range(10):
r = randomAllocation()
i = random.randint(0, M - 1)
j = last(r, i)
print("- r =", r, "i =", i, "and last(r, i) =", j)
print(" Stay on", r[i], "or go to", r[j], "?")
We need a function to simulate the first choice of the player, and a simple choice is to select a uniform choice:
In [8]:
def firstChoice():
global M
# Uniform first choice
return random.randint(0, M - 1)
Now we can simulate a game, for a certain left-to-be-written function strategy() that decides to keep or to change the initial choice.
In [9]:
def simulate(stayOrNot):
# Random spot for the treasure
r = randomAllocation()
# Initial choice
i = firstChoice()
# Which door are remove, or equivalently which is the last one to be there?
j = last(r, i)
assert {r[i], r[j]} == {False, True} # There is still the treasure and only one
stay = stayOrNot()
if stay:
return r[i]
else:
return r[j]
We can simulate many outcome of the game for one strategy, and return the number of time it won (i.e. average number of time it found the good door, by finding r[i] = True or r[j] = True):
In [10]:
N = 10000
def simulateManyGames(stayOrNot):
global N
results = [simulate(stayOrNot) for _ in range(N)]
return sum(results)
In [11]:
def keep():
return True # True == also stay on our first choice
In [12]:
rate = simulateManyGames(keep)
print("- For", N, "simulations, the strategy 'keep' has won", rate, "of the trials...")
print(" ==> proportion = {:.2%}.".format(rate / float(N)))
$\implies$ We find a chance of winning of about $\frac{1}{M} = \frac{1}{3}$ for this strategy, which is very logical as only the initial choice matters, and due to the uniform location of the treasure behind the $M = 3$ doors, and the uniform first choice with firstChoice().
In [13]:
def change():
return False # False == never stay, ie always chose the last door
In [14]:
rate = simulateManyGames(change)
print("- For", N, "simulations, the strategy 'change' has won", rate, " of the trials...")
print(" ==> proportion = {:.2%}.".format(rate / float(N)))
$\implies$ We find a chance of winning of about $\frac{M - 1}{M} = \frac{2}{3}$ for this strategy, which is less logical.
Due to the uniform location of the treasure behind the $M = 3$ doors, and the uniform first choice with firstChoice(), we have a $\frac{1}{M}$ chance of finding the treasure from the first time.
The first case has probability $\frac{1}{M}$, probability of loosing, and the second case has probability $\frac{M - 1}{M}$, probability of loosing.
$\implies$ Conclusion : this strategy change() has a chance of winning of $\frac{2}{3}$, far better than the chance of $\frac{1}{3}$ for stay().
We proved numerically the results given and explained here on the Wikipedia page. Great!
In [15]:
def bernoulli(p=0.5):
return random.random() < p
In [16]:
rate = simulateManyGames(bernoulli)
print("- For", N, "simulations, the strategy 'bernoulli' has won", rate, " of the trials...")
print(" ==> proportion = {:.2%}.".format(rate / float(N)))
Now we can try different values for $p$, and plot the resulting chance of winning the game as a function of $p$. Hopefully, it should be monotonic, confirming the result explained above.
In [17]:
import numpy as np
import matplotlib.pyplot as plt
We generate lots of values for $p$, then a function stratBernoulli() to create the strategy described above, for some $p \in [0, 1]$.
In [18]:
values_p = np.linspace(0, 1, 500)
def stratBernoulli(p):
def stayOrNot():
return bernoulli(p=p)
return stayOrNot
Let finally do all the simulations, and store the empirical probability of winning the game when following a Bernoulli strategy.
This line takes about $4$ minutes on my laptop, it's not that quick.
In [19]:
chance_of_winning = [simulateManyGames(stratBernoulli(p)) / float(N) for p in values_p]
In [20]:
plt.figure()
plt.plot(values_p, chance_of_winning, 'r')
plt.title("Monty-Hall paradox with {} doors ({} random simulation)".format(M, N))
plt.xlabel("Probability $p$ of staying on our first choice (Bernoulli strategy)")
plt.ylabel("Probability of winning")
plt.ylim(0, 1)
plt.yticks(np.linspace(0, 1, 11))
plt.show()
In [27]:
def completeSimu():
global M
global N
allocation = [False] * (M - 1) + [True] # Only 1 treasure!
def randomAllocation():
r = allocation[:]
random.shuffle(r)
return r
def last(r, i):
# Select a random index corresponding of the door we keep
if r[i]: # She found the treasure, returning a random last door
return random.choice([j for (j, v) in enumerate(r) if j != i])
else: # She didn't find the treasure, returning the treasure door
# Indeed, the game only removes door that don't contain the treasure
return random.choice([j for (j, v) in enumerate(r) if j != i and v])
def simulate(stayOrNot):
# Random spot for the treasure
r = randomAllocation()
# Initial choice
i = firstChoice()
# Which door are remove, or equivalently which is the last one to be there?
j = last(r, i)
stay = stayOrNot()
if stay:
return r[i]
else:
return r[j]
def simulateManyGames(stayOrNot):
global N
results = [simulate(stayOrNot) for _ in range(N)]
return sum(results)
values_p = np.linspace(0, 1, 300)
chance_of_winning = [simulateManyGames(stratBernoulli(p)) / float(N) for p in values_p]
plt.figure()
plt.plot(values_p, chance_of_winning, 'r')
plt.title("Monty-Hall paradox with {} doors ({} random simulation)".format(M, N))
plt.xlabel("Probability $p$ of staying on our first choice (Bernoulli strategy)")
plt.ylabel("Probability of winning")
plt.ylim(0, 1)
plt.yticks(np.linspace(0, 1, 11))
plt.show()
In [28]:
M = 4
completeSimu()
In [29]:
M = 100
completeSimu()
$\implies$ it is now clear that the best strategy is the always change and select the other door.
As we saw it, it is clear from the numerical experiments presented above that the strategy of always changing our first choice (change()), is statistically better than the strategy of keeping our first choice (keep()).
See this repository for other Python notebook doing numerical simulations.