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:

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

We will start by having you pick one door, but

**do not open it**.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!)

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) it doesn't matter?

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 [1]:

```
# 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
car_is_behind = random.choice(doors.keys())
doors[car_is_behind] = 'car'
return doors
```

In [2]:

```
print create_doors()
print create_doors()
print create_doors()
```

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 [3]:

```
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 :(
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 [4]:

```
# ok, let's try this out!
print always_switch()
print always_switch()
print always_switch()
print never_switch()
print never_switch()
print never_switch()
```

OK, a couple of notes here.

First, this *is* a true implementation of the Monty Hall scenario, right? Even though you're randomly choosing doors, it's no different from what you'd normally be doing, right -- guessing?

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

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

In [5]:

```
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'
```

We can do the same with never switching:

In [6]:

```
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'
```

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

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

In [7]:

```
# 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
car_is_behind = random.choice(doors.keys())
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 [8]:

```
def always_switch2(M):
# set up the problem
doors_d = create_multi_doors(M)
# pick a door
doors = 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 = 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 [9]:

```
print always_switch2(3)
```

In [10]:

```
print always_switch2(100)
print always_switch2(100)
print always_switch2(100)
print never_switch2(100)
print never_switch2(100)
print never_switch2(100)
```

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

In [11]:

```
won_count = 0
N = 1000
M = 3
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'
```

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

What's the actual mathematical formula underlying the results?

In [12]:

```
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)
```

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 [13]:

```
import matplotlib.mlab as mlab
# 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 = hist(always_list, facecolor='green', alpha=0.5, range=(0, 1), bins=50, normed=1)
always_avg, always_stddev = calc_stddev(always_list)
y1 = mlab.normpdf(bins1, always_avg, always_stddev)
plot(bins1, y1, 'g--')
# get the actual and fitted distribution for the 'never switch' decision
n2, bins2, patches2 = hist(never_list, facecolor='red', alpha=0.5, range=(0, 1), bins=50, normed=1)
never_avg, never_stddev = calc_stddev(never_list)
y2 = mlab.normpdf(bins2, never_avg, never_stddev) # fit
plot(bins2, y2, 'r--')
# label the plot
xlabel('probability of winning')
ylabel('distribution of p for multiple trials')
title('Monty Hall problem (%d doors)' % M)
```

Out[13]:

In [14]:

```
plot(bins1, y1, 'g', label='never switch')
plot(bins2, y2, 'r', label='always switch')
axis(ymax=12)
legend()
xlabel('probability of winning')
ylabel('distribution of p for multiple trials')
title('Monty Hall problem (%d doors)' % M)
print 'average win rate, always switching:', always_avg
print 'average win rate, never switch:', never_avg
```

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