Steve, Monty, and Julia

I don’t make a habit of reading obituaries, but I’m glad I clicked the link to this one in yesterday’s New York Times (that should be a gift link). Steve Selvin was a statistician who got his Ph.D. at the University of California, Berkeley, and taught there for many years. He got into the Times because of a probability puzzle he came up with in 1975: the Monty Hall Problem.

I won’t explain the problem here. Most of you already know it, and those who don’t can read the setup in the Wikipedia article. I do want to praise the Times for the cute little interactive game they included in Selvin’s obituary. It allows you to be a Let’s Make a Deal contestant and run through the problem again and again (with your computer playing the Monty Hall role) until you’re convinced that the “switching my choice” strategy is distinctly better than the “sticking with my original choice” strategy. I doubt you have the patience to play it enough times to confirm that switching wins 2/3 of the time and sticking wins 1/3 of the time.

I was less pleased with the explanation the Times gave for the probabilities:

After the first reveal, you have two doors, so you should have a 50-50 chance. Right? But you don’t. The host hasn’t randomly eliminated a door; he has deliberately revealed a goat, using his knowledge of where the car is. Your original door still has a one-in-three chance of hiding the car, while switching gives you a two-in-three chance of winning.

It’s right, of course, but it doesn’t strike me as an actual explanation—it sounds more like a restatement of the answer. To me, the way to think of it is to consider the winning strategy:

If you chose the door with a car behind it on your first guess—which will happen 1/3 of the time—you will win by sticking with that door. If you chose one of the doors with a goat behind it on your first guess—which will happen 2/3 of the time—you will win by switching to the door that Monty didn’t open. Therefore, switching gives you a 2/3 chance of winning.

Everyone’s different, so maybe you don’t need those extra words and the Times explanation is sufficient. Not for me, though.

I mentioned above that although you could play the game in the obituary enough times to get a decent estimate of the probabilities, most of us don’t have that kind of patience. The better way to run lots of simulations in a short time is to use a random number generator to have the computer do all the selecting. I’ve been meaning to try out the Julia programming language for quite a while; this seemed like a good problem to start with.

Here’s my script, switch.jl, for running 10,000 simulations using the “switch” strategy:

 1:  #!/usr/bin/env julia
 2:  
 3:  # Parameters
 4:  N = 10_000    # Number of trials
 5:  m = 0         # Number of wins
 6:  
 7:  # Loop through N trials, counting the number of wins
 8:  for i in 1:N
 9:    # Behind-the-scenes preparation
10:    carDoor = rand(1:3)
11:    goatDoors = [j for j in 1:3 if j != carDoor]
12:    
13:    # Player makes initial choice
14:    playerDoor = rand(1:3)
15:    
16:    # Monty chooses an available goat door
17:    if playerDoor == carDoor
18:      # Pick either of the two available
19:      montyDoor = goatDoors[rand(1:2)]
20:    else
21:      # Pick the one available
22:      montyDoor = [j for j in goatDoors if j != playerDoor][1]
23:    end
24:    
25:    # Player switches to the other door
26:    playerDoor = [j for j in 1:3 if (j != montyDoor) && (j != playerDoor)][1]
27:    
28:    # Increment m if the player wins
29:    if playerDoor == carDoor
30:      global m += 1
31:    end
32:  end
33:  
34:  # Print the winning fraction for this strategy
35:  println(m/N)

Because I’m new to Julia, I put in lots of comments, so I doubt there’s much need for more explanation. I will say, though, that I suspect a real Julia programmer wouldn’t use list comprehensions as often as I did. There’s probably a function for removing items from a list that’s more elegant than what I did in Lines 11, 22, and 26. I just don’t know what that is yet.

My script for running 10,000 simulations of the “stick” strategy, cleverly named stick.jl, is basically the same:

 1:  #!/usr/bin/env julia
 2:  
 3:  # Parameters
 4:  N = 10_000    # Number of trials
 5:  m = 0         # Number of wins
 6:  
 7:  # Loop through the N trials, counting the number of wins
 8:  for i in 1:N
 9:    # Behind-the-scenes preparation
10:    carDoor = rand(1:3)
11:    goatDoors = [j for j in 1:3 if j != carDoor]
12:    
13:    # Player makes initial choice
14:    playerDoor = rand(1:3)
15:    
16:    # Monty chooses an available goat door
17:    if playerDoor == carDoor
18:      # Pick either of the two available
19:      montyDoor = goatDoors[rand(1:2)]
20:    else
21:      # Pick the one available
22:      montyDoor = [j for j in goatDoors if j != playerDoor][1]
23:    end
24:    
25:    # Player sticks with original choice; do nothing
26:    playerDoor = playerDoor
27:    
28:    # Increment m if the player wins
29:    if playerDoor == carDoor
30:      global m += 1
31:    end
32:  end
33:  
34:  # Print the winning fraction for this strategy
35:  println(m/N)

The only difference is in Lines 25–26, where I’m basically doing a no-op.

These scripts give different answers each time they’re run because the random number generator starts with a different seed, but the answers are always close to 0.67 for switch.jl and 0.33 for stick.jl. As expected.

Are 10,000 simulations enough? Well, the 95% confidence interval can be calculated through this formula:

±1.96p(1−p)N

which works out to ±0.009 for both strategies. So 10,000 simulations is fine as long as we’re willing to accept the second decimal place to be off by one every now and then.