First, a coin toss simulation to refresh our memories before diving into Monty Hall. All materials here belong to the Harvard Data Science class. I am only using these for my self-study purposes.
Here's a way to quickly simulate 500 coin "fair" coin tosses (where the probabily of getting Heads is 50%, or 0.5)
import matplotlib.pyplot as plt
import numpy as np
x = np.random.binomial(500, .5)
print("number of heads:", x)
Repeat this simulation 500 times, and use the plt.hist() function to plot a histogram of the number of Heads (1s) in each simulation
#your code here
#this line prepares IPython for working with matplotlib
%matplotlib inline
heads = []
# 3 ways to run the simulations
# loop
# for i in range(500):
# heads.append(np.random.binomial(500, .5))
# "list comprehension"
# heads = [np.random.binomial(500, .5) for i in range(500)]
# pure numpy
heads = np.random.binomial(500, .5, size=500)
histogram = plt.hist(heads, bins=10)
Here's a fun and perhaps surprising statistical riddle, and a good way to get some practice writing python functions
In a gameshow, contestants try to guess which of 3 closed doors contain a cash prize (goats are behind the other two doors). Of course, the odds of choosing the correct door are 1 in 3. As a twist, the host of the show occasionally opens a door after a contestant makes his or her choice. This door is always one of the two the contestant did not pick, and is also always one of the goat doors (note that it is always possible to do this, since there are two goat doors). At this point, the contestant has the option of keeping his or her original choice, or swtiching to the other unopened door. The question is: is there any benefit to switching doors? The answer surprises many people who haven't heard the question before.
We can answer the problem by running simulations in Python. We'll do it in several parts.
First, write a function called simulate_prizedoor
. This function will simulate the location of the prize in many games -- see the detailed specification below:
"""
Function
--------
simulate_prizedoor
Generate a random array of 0s, 1s, and 2s, representing
hiding a prize between door 0, door 1, and door 2
Parameters
----------
nsim : int
The number of simulations to run
Returns
-------
sims : array
Random array of 0s, 1s, and 2s
Example
-------
>>> print simulate_prizedoor(3)
array([0, 0, 2])
"""
def simulate_prizedoor(nsim):
#compute here
return answer
#your code here
def simulate_prizedoor(nsim):
return np.random.randint(0,3, (nsim))
Next, write a function that simulates the contestant's guesses for nsim
simulations. Call this function simulate_guess
. The specs:
"""
Function
--------
simulate_guess
Return any strategy for guessing which door a prize is behind. This
could be a random strategy, one that always guesses 2, whatever.
Parameters
----------
nsim : int
The number of simulations to generate guesses for
Returns
-------
guesses : array
An array of guesses. Each guess is a 0, 1, or 2
Example
-------
>>> print simulate_guess(5)
array([0, 0, 0, 0, 0])
"""
#your code here
#def simulate_guess(nsim):
# return np.random.randint(0,3,(nsim))
def simulate_guess(nsim):
return np.zeros(nsim, dtype=np.int)
Next, write a function, goat_door
, to simulate randomly revealing one of the goat doors that a contestant didn't pick.
"""
Function
--------
goat_door
Simulate the opening of a "goat door" that doesn't contain the prize,
and is different from the contestants guess
Parameters
----------
prizedoors : array
The door that the prize is behind in each simulation
guesses : array
THe door that the contestant guessed in each simulation
Returns
-------
goats : array
The goat door that is opened for each simulation. Each item is 0, 1, or 2, and is different
from both prizedoors and guesses
Examples
--------
>>> print goat_door(np.array([0, 1, 2]), np.array([1, 1, 1]))
>>> array([2, 2, 0])
"""
#your code here
def goat_door(prizedoors, guesses):
#strategy: generate random answers, and
#keep updating until they satisfy the rule
#that they aren't a prizedoor or a guess
result = np.random.randint(0, 3, prizedoors.size)
while True:
bad = (result == prizedoors) | (result == guesses)
if not bad.any():
return result
result[bad] = np.random.randint(0, 3, bad.sum())
Write a function, switch_guess
, that represents the strategy of always switching a guess after the goat door is opened.
"""
Function
--------
switch_guess
The strategy that always switches a guess after the goat door is opened
Parameters
----------
guesses : array
Array of original guesses, for each simulation
goatdoors : array
Array of revealed goat doors for each simulation
Returns
-------
The new door after switching. Should be different from both guesses and goatdoors
Examples
--------
>>> print switch_guess(np.array([0, 1, 2]), np.array([1, 2, 1]))
>>> array([2, 0, 0])
"""
#your code here
def switch_guess(guesses, goatdoors):
result = np.zeros(guesses.size)
switch = {(0, 1): 2, (0, 2): 1, (1, 0): 2, (1, 2): 1, (2, 0): 1, (2, 1): 0}
for i in [0, 1, 2]:
for j in [0, 1, 2]:
mask = (guesses == i) & (goatdoors == j)
if not mask.any():
continue
result = np.where(mask, np.ones_like(result) * switch[(i, j)], result)
return result
Last function: write a win_percentage
function that takes an array of guesses
and prizedoors
, and returns the percent of correct guesses
"""
Function
--------
win_percentage
Calculate the percent of times that a simulation of guesses is correct
Parameters
-----------
guesses : array
Guesses for each simulation
prizedoors : array
Location of prize for each simulation
Returns
--------
percentage : number between 0 and 100
The win percentage
Examples
---------
>>> print win_percentage(np.array([0, 1, 2]), np.array([0, 0, 0]))
33.333
"""
#your code here
def win_percentage(guesses, prizedoors):
return 100 * (guesses == prizedoors).mean()
Now, put it together. Simulate 10000 games where contestant keeps his original guess, and 10000 games where the contestant switches his door after a goat door is revealed. Compute the percentage of time the contestant wins under either strategy. Is one strategy better than the other?
#your code here
nsim = 10000
#keep guesses
print("Win percentage when keeping original door")
print(win_percentage(simulate_prizedoor(nsim), simulate_guess(nsim)))
#switch
pd = simulate_prizedoor(nsim)
guess = simulate_guess(nsim)
goats = goat_door(pd, guess)
guess = switch_guess(guess, goats)
print("Win percentage when switching doors")
print(win_percentage(pd, guess).mean())
Many people find this answer counter-intuitive (famously, PhD mathematicians have incorrectly claimed the result must be wrong. Clearly, none of them knew Python).
One of the best ways to build intuition about why opening a Goat door affects the odds is to re-run the experiment with 100 doors and one prize. If the game show host opens 98 goat doors after you make your initial selection, would you want to keep your first pick or switch? Can you generalize your simulation code to handle the case of n
doors?