Reminder: use Shift-ENTER to run each cell; feel free to edit, etc, before and after running!

The Monty Hall problem

You, lucky contestant, are on a game show! You stand a chance to win a car (with a trade-in value of $50k)!

The host, one Monty Hall, describes the game to you as follows:

  1. Here before you, you see three doors -- A, B, and C! Behind one of them is the car.

  2. I will start by having you pick one door, but I do not open it.

  3. Then, I will reveal what is behind one of the two doors you did not choose (needless to say, it will not be the car!)

  4. At that point, you may either keep your first choice, OR you may switch your choice to another door.

If you can identify which door has the car behind it, you get to keep the car (or the money).

The question for you is this: as a winning strategy, should you:

(a) always stay with your original door choice

(b) always switch

or

(c) does it not matter?

Monty Hall Monte Carlo

You, being an accomplished computational scientist, decide that rather than thinking hard about this, you're just going to simulate it using the Monte Carlo method (http://en.wikipedia.org/wiki/Monte_Carlo_method). This lets you explore the more complex emergent properties of this system by only worrying about the simple rules defining the system. (http://en.wikipedia.org/wiki/Agent-based_model)

Let's start by setting up the problem.

In [34]:
# grab the pseudo-random # generator module in Python: http://blog.doughellmann.com/2010/10/pymotw-random-pseudorandom-number.html
import random

# now, build a function to create 3 doors, with one car behind them.
# this function will return a dictionary with keys A, B, and C, containing associated values 'empty' or 'car'.
def create_doors():
    doors = {}
    doors['A'] = 'empty'
    doors['B'] = 'empty'
    doors['C'] = 'empty'
    
    # randomly choose *one* of the three doors
    keylist = list(doors.keys())
    car_is_behind = random.choice(keylist)
    doors[car_is_behind] = 'car'
    
    return doors

Let's take a look at the output of this function:

In [35]:
print(create_doors())
print(create_doors())
print(create_doors())
{'A': 'empty', 'B': 'empty', 'C': 'car'}
{'A': 'empty', 'B': 'car', 'C': 'empty'}
{'A': 'car', 'B': 'empty', 'C': 'empty'}

Simulating some choices with two functions

Next, here are two functions, one which simulates choice (a) -- staying with your original door -- and one which simulates choice (b) -- switch. If the function returns true, you've won! If the function returns false, you've lost.

In [36]:
def always_switch():
    # set up the problem
    doors_d = create_doors()
    
    # pick a door
    doors = ['A', 'B', 'C']
    my_choice = random.choice(doors)
    
    # remove it from Monty's consideration -- he will never choose this one
    doors.remove(my_choice)
    assert len(doors) == 2, "you should only have two doors left..."
    
    # now Monty Hall picks a door:
    while 1:
        monty_choice = random.choice(doors)
        if doors_d[monty_choice] != 'car': # he'll never reveal the car!
            break
    
    doors.remove(monty_choice)
    
    # now, because you always switch, you're left with monty's non-choice:
    assert len(doors) == 1, "you should only have one door left..."
    
    my_choice = doors[0]
    
    if doors_d[my_choice] == 'car':
        return True                    # you win!
    return False                       # you lose :(
In [37]:
def never_switch():
    # set up the problem
    doors_d = create_doors()
    
    # pick a door
    doors = ['A', 'B', 'C']
    my_choice = random.choice(doors)
    doors.remove(my_choice)
    assert len(doors) == 2, "you should only have two doors left..."
    
    # now Monty Hall picks a door:
    while 1:
        monty_choice = random.choice(doors)
        if doors_d[monty_choice] != 'car': # he'll never reveal the car!
            break
    
    doors.remove(monty_choice)
    
    # now, because you never switch, you're left with your original choice:
    assert len(doors) == 1, "you should only have one door left..."
    
    # you stick with your original choice:
    if doors_d[my_choice] == 'car':
        return True                    # you win!
    return False                       # you lose :(
In [38]:
# ok, let's try this out!
print(always_switch())
print(always_switch())
print(always_switch())

print(never_switch())
print(never_switch())
print(never_switch())
True
False
True
False
False
True

OK, a couple of things:

Question: Is this a true implementation of the Monty Hall scenario? Why or hwy not?

Second, note that every time you run this, you will get a different set of answers -- because we're using a random number generator.

The First Question: is it better to always switch? Or should you never switch? Or is it the same?

To answer this question, let's run 'always_switch()' a bunch of times, and calculate how often we win.

In [39]:
won_count = 0
N = 1000
for i in range(N):
    if always_switch():
        won_count += 1
        
print('With always_switch(), I won', won_count, 'of', N, 'tries')
With always_switch(), I won 667 of 1000 tries

We can do the same with never switching:

In [40]:
won_count = 0
N = 1000
for i in range(N):
    if never_switch():
        won_count += 1
        
print('With never_switch(), I won', won_count, 'of', N, 'tries')
With never_switch(), I won 319 of 1000 tries

OK, so according to our simulation, we should always switch. ...can we get this answer some other way?

Another question: what happens if you increase the number of doors, and Monty opens n-1 doors?

Since all our doors are virtual, we can change their number quite easily. (Well, it involves a little more programming...)

In [41]:
# write a function to create M doors, not just 3.
def create_multi_doors(M):
    doors = {}
    for i in range(M):
        doors[i] = 'empty'
        
    # randomly choose *one* of the three doors
    keylist = list(doors.keys())
    car_is_behind = random.choice(keylist)
    doors[car_is_behind] = 'car'
    
    return doors

We'll also need to update our always_switch/never_switch functions to use create_multi_doors():

In [42]:
def always_switch2(M):
    # set up the problem
    doors_d = create_multi_doors(M)
    
    # pick a door
    doors = list(doors_d.keys())
    my_choice = random.choice(doors)
    
    doors.remove(my_choice)
    
    while len(doors) > 1:
        door = random.choice(doors)
        if doors_d[door] != 'car':
            doors.remove(door)
    
    assert len(doors) == 1, doors
    
    # now pick the one that's left:
    my_choice = doors[0]
    
    if doors_d[my_choice] == 'car':
        return True                    # you win!
    return False                       # you lose :(

def never_switch2(M):
    # set up the problem
    doors_d = create_multi_doors(M)
    
    # pick a door
    doors = list(doors_d.keys())
    my_choice = random.choice(doors)
    
    doors.remove(my_choice)
    
    while len(doors) > 1:
        door = random.choice(doors)
        if doors_d[door] != 'car':
            doors.remove(door)
    
    assert len(doors) == 1, doors
    
    # ...but we're sticking with our original choice.    
    if doors_d[my_choice] == 'car':
        return True                    # you win!
    return False                       # you lose :(
In [43]:
print(always_switch2(3))
True
In [44]:
print(always_switch2(100))
print(always_switch2(100))
print(always_switch2(100))

print(never_switch2(100))
print(never_switch2(100))
print(never_switch2(100))
True
True
True
False
False
False

OK, the functions don't obviously break... what happens if we run them a bunch?

In [45]:
won_count = 0
N = 1000 # num trials
M = 3 # doors

for i in range(N):
    if always_switch2(M):
        won_count += 1
        
print('With always_switch2() and', M, 'doors, I won', won_count, 'of', N, 'tries')

won_count = 0

for i in range(N):
    if never_switch2(M):
        won_count += 1
        
print('With never_switch2() and', M, 'doors, I won', won_count, 'of', N, 'tries')
With always_switch2() and 3 doors, I won 641 of 1000 tries
With never_switch2() and 3 doors, I won 359 of 1000 tries

A few questions to explore --

How do we know always_switch2() and never_switch2() work properly?

What's the actual mathematical formula underlying the results?

Visualization!

In [46]:
import math

# calculate the avg & standard deviation
def calc_stddev(series):
    avg = sum(series) / float(len(series))
    devs = [ i - avg for i in series ]
    devs = [ i**2 for i in devs ]
    stddev = math.sqrt(sum(devs) / float(len(series)))
    return avg, stddev


# run the given function with parameter M (# doors), N times; return average
def try_N_times(fn, M, N):
    won = 0
    for i in range(N):
        if fn(M):
            won += 1
    return won / float(N)

Explore the simulation

What does N represent and what happens as you change N in different ways?

(Answer/Notes here. Double-click to edit.)

What does n_trials represent and what happens as you change n_trials in different ways?

(Answer/Notes here. Double-click to edit.)

What does M represent and what happens as you change M in different ways?

(Answer/Notes here. Double-click to edit.)

In [47]:
%matplotlib inline
import matplotlib.mlab as mlab
import pylab

# Parameters we've made easy to change
#################
n_trials = 500  # How many times to gather statics from groups of N samples
M = 3           # How many doors to use (originally 3)
N = 100         # How many samples from which to gather statistics (like mean of N coin flips)
#################

# Try M doors N times; then run n_trials times, and track.
always_list = [ try_N_times(always_switch2, M, N) for _ in range(n_trials) ]
never_list = [ try_N_times(never_switch2, M, N) for _ in range(n_trials) ]

# get the actual and fitted distribution for the 'always switch' decision
n1, bins1, patches1 = pylab.hist(always_list, facecolor='green', alpha=0.5, range=(0, 1), bins=50, normed=1, label='always')
always_avg, always_stddev = calc_stddev(always_list)
y1 = mlab.normpdf(bins1, always_avg, always_stddev)
pylab.plot(bins1, y1, 'g--')

# get the actual and fitted distribution for the 'never switch' decision
n2, bins2, patches2 = pylab.hist(never_list, facecolor='red', alpha=0.5, range=(0, 1), bins=50, normed=1, label='never')
never_avg, never_stddev = calc_stddev(never_list)
y2 = mlab.normpdf(bins2, never_avg, never_stddev) # fit
pylab.plot(bins2, y2, 'r--')

# label the plot
pylab.xlabel('probability of winning')
pylab.ylabel('distribution of p for multiple trials')
pylab.title('Monty Hall problem (%d doors)' % M)
pylab.legend(loc='upper center')
Out[47]:
<matplotlib.legend.Legend at 0x106993978>
In [48]:
pylab.plot(bins1, y1, 'g', label='never switch')
pylab.plot(bins2, y2, 'r', label='always switch')
pylab.axis(ymax=12)
pylab.legend()
pylab.xlabel('probability of winning')
pylab.ylabel('distribution of p for multiple trials')
pylab.title('Monty Hall problem (%d doors)' % M)
print('average win rate, always switching:', always_avg)
print('average win rate, never switch:', never_avg)
average win rate, always switching: 0.6697599999999995
average win rate, never switch: 0.33708000000000027

Bonus Question: What are the assumptions we're making in this simulation compared to the real game show?