Monty Hall problem, Statistics, Python and solving problems differently

I have heard this problem many times and with lots of variations. A participant in a game show must choose one out of three doors. In one of them there is a valuable reward while in the other 2, nothing. The participant chooses one door. The host of the game show (who knows what is behind each door) must reveal an empty door (1 of the remaining ones) for the participant. Then he asks the participant:

Do you want to change the door you chose?

The answer is yes. You want to change the door you chose to the third door (the one you didn't open and the host didn't reveal).

But why?

I have seen couple of videos about it, and frankly, I still didn't understand completely.

In the same time, I have just finished the book 'Think Stats' version 2 by Allen B. Downey. (I added a link for PDF below). This book is about understanding statistics with coding and not about formulas. To understand pvalue, regression and correlation, instead of using formulas, we can just understand what it means and perform tests on data which generate the desired result.

So, I wrote a Python script that simulates the Monty Hall problem and another code that analyzes the results.

The amazing thing is that while I was writing the code, I was thinking about the problem differently. I didn't fall for the 'now we have a new problem with 50-50 chance for each door' trap. I start realizing that some of the lines of the code are not contributing any information. We don't need all this information. When my code was completed I realized the essence of this problem. The host didn't give the participant any information. He just revealed a door that didn't have the reward in it. But the participant already knew there is a door he didn't choose with no prize. So, the chance of switching door is 67% and not 50% because no new information was given to the participant. It is as if the participant can choose 1 door with 33% chance of win or 2 doors with 67% chance of win.

This realization came to me through writing code and not thinking about statistics. I just thought 'this line doesn't give any meaningful information. That line doesn't change anything'. In my mind, I reduced the problem to the original one.

The code that simulates the Monty Hall problem is monty_hall_problem.py (added a link to GitHub below):

import random


def simulate_monty_hall(n=1000):


    def monty_hall():
        doors = [True, False, False] # generate doors
        random.shuffle(doors)
        c1 = random.randint(0,2)
        chosen_door = doors[c1] # choose first door
        doors.pop(c1) 
        falses = [] # doors the host can show the participent
        for idx, val in enumerate(doors):
            if val == False:
                falses.append(idx)
        false_door = random.choice(falses) # the host show random false door to the participent
        doors.pop(false_door)
        alternative_door = doors[0] # the remaining door
        return chosen_door, alternative_door


    results = []
    for _ in range(n):
        results.append(monty_hall())


    counter_for_chosen = 0
    for result in results:
        if result[0] == True:
            counter_for_chosen += 1


    return counter_for_chosen/n

1. Create a random list of the True, False, False.
2. Shuffle the list.
3. The participant chooses a door (now we have 2 doors available for revealing).
4. The host reveal 1 random False door to the participant (now we have 1 door remaining for switching).
5. The monty_hall function returns a tuple of the result (True or False) of the chosen door and the alternative one.
6. Sum for n =1000 tests how many times we have True behind the chosen door. Then divide this number by n. The result is the success ratio when the participant didn't change the door.

If you run this code, you will get approximately 0.33 success ratio.

That could be enough, but I wanted more. I wanted to analyze what the STD of this test is if I run it 1000 times. What is its histogram? What is the analytical solution?

The answers to these questions are in the file analysis.py.

I will not show the code here because it is a little bit long. But I will show here the results :) You can view the code in GitHub.

First, lets explain the difference between n and iter in the code:

1. n - number of independent games the participate plays. Out of it, we can calculate the success ratio - should be about 1/3. When n is bigger - the closer we get to 1/3.
2. iter - number of tests we conduct to get the distribution of the success ratio. Out of it, we have an object with length iter which each value is the success ratio of n independent games. When iter is bigger - the closer we get to success ratio of 1/3.

The analysis has 5 outputs for pairs of n and iter:

1. Histogram of the distribution
2. Cumulative distribution function of the distribution
3. Normal probability plot of the distribution
4. Analytical result - depends only on n (Mean, STD)
5. Estimated results (Mean, STD, Confidence Interval, Median)

Lets review each one of them

Histogram

As you might assume, the histogram looks more and more normal when n and iter increases. Interestingly, iter is more influential. When n=10 and iter=1000 the hist looks more normal than when n=1000 and iter =10.

One more thing to note here is that the x and y axes are different in each plot. Just waned to show the distribution itself. Not necessarily compare each other.

Cumulative distribution function

These plots can be directly compared to each other. As you might guess, the more n and iter are bigger, the more normal the distribution is.

Normal probability plot

This is an interesting one. I didn't know this method until I read the book. You take m random numbers from a normal distribution and m random numbers from an existing distribution. You sort both of the lists. Then you plot on x axis the values from the list generated from the normal distribution and on the y axis the values from your existing distribution. The pairs are drawn according to their position in the list (and both our lists are m in length).

If the resulted plot is a line - it is probable that your existing distribution is a normal distributed. The intercept is the mean of the distribution and the slope is the STD.

Estimated results

You guessed it! The more n and iter are bigger, the more we are close to success ratio 1/3. But you might notice something else, the STD is getting smaller and smaller when n increases also.

The confidence interval of the mean (90%) is derived from the distribution itself and not from an analytical formula.

Analytical result

Here the mean (and median) is 1/3 for all. The STD is derived from:

STD = Sqrt(p * (1 - p) / n)

When p is the success ratio and n is the number of games.

I will not elaborate about it but you can find a link with an explanation below.

Now, compare these numbers to the numbers on iter=1000 from the estimated results. The differences are very small.

Conclusion

I hope you enjoyed this article. I really tried to stay away from complicated concepts and to present to you a different way to look and solve things. If you come across an analytical problem you find it hard to solve, try programming it. You will see things from other perspective.

Also, I encourage you to read the book 'Think Stats'. If you know basic (sometimes it is more than basic but it is worth to learn it on the way) Python and want to understand statistics - you should definitely give it a try.

Links

Book: http://greenteapress.com/thinkstats2/thinkstats2.pdf

GitHub: https://github.com/shteken/monty-hall-problem-analysis

STD explanation: https://en.wikipedia.org/wiki/Binomial_distribution#Normal_approximation