# coding: utf-8
# # Numerical simulations of the Monty-Hall "paradox"
# This short notebook aims at simulating trials of the so-called [Monty-Hall problem](https://en.wikipedia.org/wiki/Monty_Hall_problem), and thus helping to convince about the result thanks to a numerical results rather than a possibly-unclear proof.
# ----
# ## Definition of the problem
#
# 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:
#
# 1. The player first chose a door, but does not open it yet,
# 2. All remaining doors but one are opened, and none of them contain the treasure. The player sees $M - 2$ bad doors,
# 3. So there is just 2 doors, the one first chosen, and the last one. She knows the treasure is behind one of them,
# 4. And she has to decide if she wants to **stay on her initial choice**, or **switches to the last door**.
#
# 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.
# ## Starting our numerical simulations
# We start by importing some modules, then define a function, `randomAllocation()`, to generate a random allocation of treasures behind the `M` doors.
#
# > *Note:* all the code below is generic for any $M \geq 3$, but `M = 3` is used to keep small and clear visualizations.
# 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
# We can generate a random spot for the treasure by simply shuffling (with [`random.shuffle()`](https://docs.python.org/3/library/random.html#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)
# ----
# ## Comparing two strategies, on many randomized trials
# We will simulate the two strategies, `keep()` vs `change()`, on $N = 10000$ randomized games.
# ### Keeping our first choice, `keep()`
# 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()`.
# ### Changing our first choice, `change()`
# 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.
#
# 1. If we found it, the last door does not contain the treasure, hence **we loose** as we switch to it,
# 2. However, if we did not find it, the last door has to contain the treasure, hence **we win** as we switch to it deterministically (i.e. always).
#
# 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](https://en.wikipedia.org/wiki/Monty_Hall_problem#Direct_calculation). Great!
# ### Bernoulli choice
# We can try a randomized strategy, a simple one can be to follow the decision of a (biased) coin:
#
# 1. Toss a coin (with fixed head probability equals to $p \in [0, 1]$)
# 2. If head, switch to the last door.
# 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()
# ---
# ## Examples with $M = 100$ doors
# 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()
# ### Last plot, for $M = 100$
# It clearly shows the linear behavior we expected.
# In[29]:
M = 100
completeSimu()
# $\implies$ it is now clear that the best strategy is the **always change and select the other door**.
# ----
# ## Conclusion
# 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](https://github.com/Naereen/notebooks/tree/master/simus/).