I thought I had familiarized myself enough with the basic workings of linear algebra after Week 1 that I thought Rachel Thomas’ (fast.ai)’s Linear algebra for coders course would be a good way to deepen my knowledge by applying these newly learned concepts in Python.

But embarrassingly, even the first warmup problem was beyond me. It required the use of a Markov matrix, which is the application of linear algebra to probability. Then I realized I had no background in probability, so I tried Brilliant’s Intro to Probability, whereupon I roundly botched a bunch of basic probability questions.

I found it really odd that I had just learned these mathematical techniques (basic linear algebra and multivariable calculus) to solve problems in a new way, yet while these probability problems could be solved computationally, they are meant to be solved “intuitively.” Of course, that seems like a cruel joke, since their counterintuitive solutions make such problems seem more like Zen koans.

The Monty Hall problem

Take, for example, the Monty Hall problem:

Suppose you’re on a game show, and you’re given the choice of three doors: Behind one door is a car; behind the others, goats.

You pick a door, say No. 1, and the host, who knows what’s behind the doors, opens another door, say No. 3, which has a goat. He then says to you, “Do you want to pick door No. 2?”

Is it to your advantage to switch your choice?

Spoiler alert: The answer is that you do switch, because you would then have a \(\frac{2}{3}\) probability of choosing the car.

Wait, what? How does a single door go from having a \(\frac{1}{3}\) chance that the car is behind it to a \(\frac{2}{3}\) chance? It seemed to me that you would now have a 50/50 chance, since there were two doors left.

People tend to explain the rationale for the solution by increasing the number of doors to 100, and the number of doors revealed to 98. But this didn’t change anything for me; you still ended up with 2 doors, and thus a 50/50 chance.

No number of doors, visual explanation, or animation seemed to help me understand it. Finally, I managed to find a comment in an old ELI5 thread that made it click for me:

As simply as possible: Don’t think of it as three doors. Think of it as your door, and Monty’s doors. The odds that you picked the right door are 1 in 3, and the odds that you didn’t are 2 in 3, right?

When Monty gets rid of one bad choice, he doesn’t change the odds that your door is right - it’s still 1 in 3. That means he’s also not changing the odds that you aren’t right - it’s still 2 in 3.

Therefore you’re not picking one door - you’re picking two doors at the same time and getting the best possible outcome. If either of Monty’s doors was right, you win; If both of Monty’s doors were bad, you lose.

Explained this way, it then makes sense that there would only be a 50/50 chance if you had two doors to begin with. Or, if you had fallen into “the probability trap of treating non-random information as if it were random,” as it has been described elsewhere.

Statistical approach

Let’s first try to simulate the game in Python. I’ll skip over explaining what’s going on in the code here, as it’s pretty elementary.

I just added perpetual loops to avoid the problem of bad input (if you were going to get other people to play, for example).

import random

def monty_hall():
    doors = ["a","b","c"]
    car = random.choice(doors)

    while True:
        first_choice = input("Choose a door (a/b/c): ")
        if first_choice in doors:
            break

    for door in doors:
        if door is not car and door is not first_choice:
            monty_opens = door
    while True:
        action = input(f"Monty opens door {monty_opens}; do you switch? (y/n) ")
        if action == "y": # switch doors
            final_choice = [door for door in doors if door is not first_choice and door is not monty_opens][0]
            break
        elif action == "n":
            final_choice = first_choice
            break
        else:
            continue

    print(f"Your final choice is door {final_choice}.")
    if final_choice == car:
        print(f"You win! The car is behind door {car}.")
    else:
        print(f"You lose. The car was behind door {car}.")

Next, let’s automate the process so that the computer decides randomly whether or not to switch, and returns us the information we need to analyze the outcomes (stay/switch and win/lose).

def monty_hall():
    doors = ["a","b","c"]
    car = random.choice(doors)
    first_choice = random.choice(doors)

    for door in doors:
        if door is not car and door is not first_choice:
            monty_opens = door

    switch_doors = random.choice([True, False])
    final_choice = [door for door in doors if door is not first_choice and door is not monty_opens][0] if switch_doors else first_choice
    result = True if final_choice == car else False
    return {"switch_doors": switch_doors, "result": result}

In [18]: monty_hall()
Out[18]: {'switch_doors': False, 'result': False}

Just as we can get the computer to flip a million coins for us, we can also get it to confront the Monty Hall problem a million times (if we want). This is called a Monte Carlo simulation. I was already familiar with this concept, but to offer a quick summary here, it basically leverages the law of large numbers: If you perform the same experiment enough times, the real result eventually converges on the theoretical result.

So let’s code it with a simple loop, tallying the results in two lists, and taking advantage of the truthiness of booleans in Python to quickly analyze the results:

def monte_monty(times):
    switch, stay = [], []
    for _ in range(times):
        x = monty_hall()
        if x["switch_doors"]:
            switch.append(x["result"])
        else:
            stay.append(x["result"])
    if switch: # avoid division by zero errors if you do 1 or 2 trials
        print(f"Switching won {sum(switch)} out of {len(switch)} times; probability of {round(sum(switch)/len(switch) * 100, 2)}%.")
    if stay:
        print(f"Staying won {sum(stay)} out of {len(stay)} times; probability of {round(sum(stay)/len(stay) * 100, 2)}%.")

To reiterate the explanation above, if you perform the same experiment enough times, the real result eventually converges on the theoretical result. Watch what happens with just 25 simulations:

In [20]: monte_monty(25)
Switching won 7 out of 12 times; probability of 58.33%.
Staying won 4 out of 13 times; probability of 30.77%.

It is already clear after just 25 simulations that the probability of switching is significantly greater than that of staying. Let’s keep going.

In [21]: monte_monty(50)
Switching won 18 out of 29 times; probability of 62.07%.
Staying won 7 out of 21 times; probability of 33.33%.

In [22]: monte_monty(100)
Switching won 32 out of 51 times; probability of 62.75%.
Staying won 16 out of 49 times; probability of 32.65%.

In [23]: monte_monty(1000)
Switching won 330 out of 496 times; probability of 66.53%.
Staying won 169 out of 504 times; probability of 33.53%.

In [24]: monte_monty(10000)
Switching won 3376 out of 5053 times; probability of 66.81%.
Staying won 1675 out of 4947 times; probability of 33.86%.

In [25]: monte_monty(100000)
Switching won 33289 out of 49759 times; probability of 66.9%.
Staying won 16789 out of 50241 times; probability of 33.42%.

In [26]: monte_monty(1000000)
Switching won 333204 out of 500463 times; probability of 66.58%.
Staying won 166298 out of 499537 times; probability of 33.29%.

In [27]: monte_monty(10000000)
Switching won 3331217 out of 4998552 times; probability of 66.64%.
Staying won 1666441 out of 5001448 times; probability of 33.32%.

After 1000 runs, we pretty much obtain the theoretical probability of \(\frac{2}{3}\) chance of a winning switch and \(\frac{1}{3}\) chance of a winning stay, and after 10 million runs, we come within a hair’s breadth of it.

In this way, the Monte Carlo method is not that different from stochastic gradient descent.

Probabilistic approach

There is another, less brute-force, more theoretical way of working out this answer using Bayesian inference.

Let’s rephrase the problem. Instead of focusing on the probability of a certain action (stay/switch) winning, we’ll examine the probability of a car being behind door \(x\) after Monty opens door \(y\), or in the language of conditional probability, given that Monty opens door \(y\). That is,

Let’s start with the unconditional probability (i.e., before the game starts) of the car being behind each of the doors, \(P(\textrm{car}_x)\). There are three doors and one car, so

Next, let’s determine the unconditional probabilities (i.e., before the game starts) of Monty choosing each of the doors. Since he cannot choose the door you will have chosen, that leaves him with two doors:

Now that we have our building blocks, let’s plug in what we can into the theorem:

where \(H\) is the hypothesis we’re interested in and \(E\) is the plot-twist, game-changing event.

Let’s assume our first choice is door A. If we want to know the probability of the car being behind door A after Monty opens door C (i.e., we win by staying), we can fill in what we know:

which leaves one last probability to calculate, \(P(\textrm{Monty}_c|\textrm{car}_a)\). That is the probability that Monty will open door C if the car is behind door A.

Monty can’t choose your door (A) or the door with the car (also A), leaving him two choices. He is therefore equally likely to choose door B or C, which means that \(P(\textrm{Monty}_c|\textrm{car}_a) = \tfrac{1}{2}\).

Though the theorem is quite convoluted, this answer agrees with our intuition and the Monte Carlo results.

What about the probability that we win by switching (i.e., the probability that the car is behind door B after Monty opens door C)?

In the case of \(P(\textrm{Monty}_c|\textrm{car}_b)\), we have chosen door A and the car is behind door B. Monty knows this, and he also can’t choose the door you’ve chosen, door A.

His only choice, then, is to choose door C, and therefore \(P(\textrm{Monty}_c|\textrm{car}_b) = 1\).

Again, the theorem agrees with the Monte Carlo results, although whether or not it also agrees with your intuition is a different question!

Where do we go from here?

I know that the Monte Carlo method can be used to make predictions about the movements of the stock market. However, I’m totally new to probability theory, so I’m not yet sure where to go after Bayes’ theorem. I’ll update this section once I find out and post about more generalized and sophisticated ways to apply Monte Carlo and Bayes’ theorem. Stay tuned!

References

• The Monty Hall problem and 3 ways to solve it, Anthony O’Farrell
• Bayes’ Theorem with Lego, Will Kurt
• How to solve the Monty Hall problem using Bayesian inference, Max Pagels

// filed under Monte Carlo Bayesian inference in Probability    // share on Twitter