18
$\begingroup$

Problem

For purely recreational purposes I would like to solve the Monty Hall problem with Mathematica using the function Probability (dedicated to the calculation of probabilities).

About the Monty Hall problem and its solution

Here is a possible formulation of the famous Monty Hall problem:

Suppose you’re given the choice of three doors: behind one door is a car, each door having the same probability of hiding it; behind the others, goats. You pick a door and the game organizer, who knows what’s behind the doors, opens another door which has a goat. They then says to you: “Do you want to pick the other door?”.

Is it to your advantage to switch your choice? Or more precisely: what is the probability that the car is behind the other door?

This is a well-known probability problem, and its solution may sometimes appear counterintuitive. The answer being: yes it is advantageous to switch your choice, the probability of finding the car behind the other door is $\frac{2}{3}$.

One way to arrive at this result is to use Bayes’ theorem. Let $C_i$ denote the event “the car is behind the door $i$”. We consider the case where door 3 has just been chosen. At this point: $P(C_1) = P(C_2) = P(C_3) = \frac{1}{3}$.

With disjunction of cases, one can notice that if the car is behind door 1, the organizer shall open door 2; if the car is behind door 2, the organizer shall open door 1; and finally if the car is behind door 3, the organizer may open either door 1 or 2 (each outcome being equiprobable).

We can then consider that door 1 has been opened by the organizer (thus discovering a goat behind it), while denoting this event $O_1$. To determine the probability that the car is behind the other door (door 2), we can calculate the conditional probability using the information we’ve just obtained:

$$ P(C_2 | O_1) = \frac{P( O_1 | C_2) P(C_2)}{P(O_1)} = \frac{P( O_1 | C_2) P(C_2)}{\sum_{i=1}^3 P(O_1 | C_i) P(C_i)} = \frac{\frac{1}{3}}{\frac{1}{2}} = \frac{2}{3}. $$

One can notice that the same reasoning applies regardless of the door chosen initially and the door opened subsequently. We can then conclude that the probability of finding the car behind the other door is always $\frac{2}{3}$.

My attempt to solve the problem with Mathematica

Obviously, it is very simple here to simulate the situation with Mathematica a large number of times in order to obtain the probability numerically. But I’m trying to solve the problem analytically using the function Probability to get an exact result.

I therefore took up the situation described above: the door 3 has been chosen, and the door 1 has been subsequently opened by the organizer, and we want to determine the probability that the winning door is the other door (door 2). So I tried:

In[1]:= Probability[
 (c == 2) \[Conditioned] (o == 1 && (c == 1 \[Implies] (o == 2)) && (c == 2 \[Implies] (o == 1))),
 {
    c \[Distributed] DiscreteUniformDistribution[{1, 3}],
    o \[Distributed] DiscreteUniformDistribution[{1, 2}]
 }
]

I considered two random variables in Mathematica: c, the number of the winning door, following a discrete uniform distribution between 1 and 3; and o, the number of the opened door, following a discrete uniform distribution between 1 and 2 (since door 3 has been chosen, it can no longer be opened). The Probability function considers a priori that these variables are independent. So I used the expression after \[Conditioned] to express the door opened by the organizer, and the link between that door and the winning door.

Unfortunately, I don’t get the expected result:

Out[1]= 1/2

I think I understand why Mathematica comes up with this output: it simplifies the expression after \[Conditioned] to o == 1 && c != 1 and eliminates information about o (since it considers the variables as independent) thus leading to the aforementioned result.

Thenceforth, I am not sure how to model the problem with the Probability function in such a way as to correctly express the link between the winning door and the opened door.

Improve this question
$\endgroup$
3
  • $\begingroup$ +1 for a correct description of the problem, so often not the case. $\endgroup$
    – mckenzm
    Jul 5, 2020 at 2:47
  • $\begingroup$ [short description]: There're 3 doors, each which may have the prize behind it. You can pick one, then the host will eliminate one of the wrong doors that you didn't pick. Then you can switch your selection from your current door to the other, but should you do it? Then, yes, because you had a 1-in-3 of your current door being correct, while there was a 2-in-3 of one of the other doors being correct, so you can get that 2-in-3 by selecting the remaining door. $\endgroup$
    – Nat
    Jul 5, 2020 at 5:10
  • $\begingroup$ Unfortunately, I do get the expected result with the given formula. That supports my consideration the question is solved with this. One might consult WolframAlpha too, with supports my answer: Switch. $\endgroup$
    – Steffen Jaeschke
    Jul 7, 2020 at 14:01

3 Answers 3

Reset to default
25
$\begingroup$

I've looked into this myself and I don't think the problem is with Mathematica. The problem is how to represent the choice of the host. Here's an attempt I tried:

So the basic idea here is: I pick a number between from 1 to 3 and so does the car. The host picks randomly between the numbers 1 and 2 and adds that number (mod 3) to mine to pick a different door than I did. Then you condition on the host's number not being the car.

So what does this give?

unif[n_] := DiscreteUniformDistribution[{1, n}];
Probability[
 Conditioned[
  myChoice == car,
  Mod[myChoice + hostChoice, 3, 1] != car
  ],
 {
  myChoice \[Distributed] unif[3],
  car \[Distributed] unif[3],
  hostChoice \[Distributed] unif[2]
  }
]

1/2

Ugh... that doesn't look right, does it? Surely something went wrong here. Let's just simulate this thing, because numbers don't lie:

simulation = AssociationThread[{"MyChoice", "Car", "HostChoice"}, #] & /@ 
  RandomVariate[
   ProductDistribution[unif[3], unif[3], unif[2]],
   10000
];
Dataset[simulation, MaxItems -> 10]

enter image description here

I'm turning the numbers into Assocations for making the code more readable. So let's do some counting:

CountsBy[
 Select[simulation, Mod[#MyChoice + #HostChoice, 3, 1] =!= #Car &],
 #MyChoice === #Car &
]
N[%/Total[%]]

<|True -> 3392, False -> 3310|>

<|True -> 0.506118, False -> 0.493882|>

Ok, so maybe Probability wasn't wrong after all. What we're seeing here is the real reason why the Monty Hall problem is difficult: the outcome depends crucially on how you model the behaviour of the host. In this description it is -in principle- possible for the host to pick the door with the car. We just condition that possibility away.

But this is different from the actual behaviour of the host: If you pick the door with the car, the host selects randomly between the two remaining doors. If you don't pick the car, the host doesn't pick randomly at all! This is where our calculation breaks down: we always assume the host picks between two doors, but that's not how it works and that's why the Monty Hall problem is trickier than it appears, even when you think you understand it.

To put it succinctly: the line hostChoice \[Distributed] unif[2] is plainly wrong. The host's choice is a combination between a deterministic choice and unif[2] that depends on myChoice.

As for the question how to reproduce the correct answer with Probability and Conditioned: I don't think that it's possible to represent this type of conditionality (i.e., the distribution of one random variable depending on another random variable) with the tools currently given. The only thing that comes close is ParameterMixtureDistribution, but I don't think that will help here.

Edit

I'm happy to let you know that I actually did manage to squeeze Monty Hall into ParameterMixtureDistribution with some torture. First of all, we will need to be able to define probability distributions such as "a random choice from the numbers in a list by weight". I defined such a distribution as follows:

Clear[discreteNumberDistribution]
discreteNumberDistribution[lst_List -> weights_List, {min_, max_}] := 
  With[{nWeights = weights/Total[weights]},
   ProbabilityDistribution[
    Sum[nWeights[[i]]*KroneckerDelta[\[FormalX], lst[[i]]], {i, Length[lst]}],
    {\[FormalX], min, max, 1}
   ]
];

So now we can do things like:

RandomVariate @ discreteNumberDistribution[{2, 3} -> {2, 10}, {1, 3}]

3 (* most likely *)

Now we can define the mixture distribution of my choice, the car and the host choice as follows:

mixture = ParameterMixtureDistribution[
  ProductDistribution[
   discreteNumberDistribution[{\[FormalM]} -> {1}, {1, 3}], (* my choice *)
   discreteNumberDistribution[{\[FormalC]} -> {1}, {1, 3}], (* car *)
   discreteNumberDistribution[ (* host choice *)
    Range[3] -> (Boole[! (\[FormalM] == # || \[FormalC] == #)] & /@ Range[3]),
    {1, 3}
   ]
  ],
  {
   \[FormalM] \[Distributed] DiscreteUniformDistribution[{1, 3}],
   \[FormalC] \[Distributed] DiscreteUniformDistribution[{1, 3}]
   }
];

So let's ask Mathematica again:

Probability[myChoice == car, {myChoice, car, host} \[Distributed] mixture]

1/3

and

Probability[
  otherChoice == car \[Conditioned] otherChoice != myChoice && otherChoice != host, 
  {
    {myChoice, car, host} \[Distributed] mixture, 
    otherChoice \[Distributed] DiscreteUniformDistribution[{1, 3}]
  }
]

2/3

Victory!

Improve this answer
$\endgroup$
6
  • $\begingroup$ Nice answer! Yes, you’re right about how difficult it is to model the situation correctly. In my question, I am convinced that Mathematica’s answer is correct, and that my Monty Hall model is not. Kudos on modeling the situation with ParameterMixtureDistribution! $\endgroup$
    – mlpo
    Jul 5, 2020 at 17:18
  • $\begingroup$ Just a note though, with this way of getting 1/3 I feel like we’re asking Mathematica for the probability of having two equal independent uniformly distributed variables between 1 and 3—here myChoice and car (which is obviously 1/3). To take into account the choice of host (and to get the answer to the problem), I’d rather use: Probability[otherChoice == car [Conditioned] otherChoice != myChoice && otherChoice != host, {{myChoice, car, host} [Distributed] mixture, otherChoice \[Distributed] DiscreteUniformDistribution[{1, 3}]}] (which correctly outputs 2/3!). What do you think about it? $\endgroup$
    – mlpo
    Jul 5, 2020 at 17:19
  • $\begingroup$ @mlpo Yeah, once you have the full multivariate distribution mixture, you can derive other properties from it. Good thinking :) $\endgroup$
    – Sjoerd Smit
    Jul 5, 2020 at 19:18
  • $\begingroup$ +1 I wonder if there's an easy way to generate a tree diagram with this or to otherwise visualize it. $\endgroup$
    – Sterling
    Oct 6, 2021 at 2:42
  • 1
    $\begingroup$ @Sterling You modify the line \[FormalC] \[Distributed] DiscreteUniformDistribution[{1, 3}] accordingly. No idea if Mathematica will actually be able to work it out, though. $\endgroup$
    – Sjoerd Smit
    Oct 7, 2021 at 7:09
5
$\begingroup$

I know that this is not what you want, but for completeness I'll add a Monte-Carlo version:

findthecar[numberofdoors_Integer /; numberofdoors >= 3] :=
  Module[{car, goats, myfirstchoice, notmyfirstchoice, organizerschoice, mysecondchoice},
    (* the car is behind a random door *)
    car = RandomInteger[{1, numberofdoors}];
    (* there are goats behind the other doors *)
    goats = Complement[Range[numberofdoors], {car}];
    (* at first I choose a random door *)
    myfirstchoice = RandomInteger[{1, numberofdoors}];
    (* these are the doors I did not choose yet *)
    notmyfirstchoice = Complement[Range[numberofdoors], {myfirstchoice}];
    (* the organizer opens a door that is not my choice and that has a goat *)
    organizerschoice = RandomChoice@Intersection[notmyfirstchoice, goats];
    (* my second choice is not my first and not the organizer's *)
    mysecondchoice = RandomChoice@Complement[Range[numberofdoors],
                                             {myfirstchoice, organizerschoice}];
    (* is the car behind my second chosen door? *)
    mysecondchoice == car]

Try it out a million times for three doors, and see that I find the car in about 2/3 of cases:

Table[findthecar[3], {10^6}] // Counts
(*    <|True -> 666122, False -> 333878|>    *)

More generally, in a game of $n$ doors where the organizer opens $k$ goat-revealing doors and I am given the option to switch, my probability for finding the car before and after switching are

$$ P_{\text{no switch}}(n,k) = \frac{1}{n},\\ P_{\text{switch}}(n,k) = \frac{n-1}{n(n-k-1)}, $$

respectively. As $P_{\text{switch}}(n,k)>P_{\text{no switch}}(n,k)$ in all cases, we should always switch after the organizer's goat-reveal.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for the answer! I was indeed looking for a symbolic solution, but the title was generic. I also created a numerical solution, and maybe I should have added it in the first place. +1 for completness :) $\endgroup$
    – mlpo
    Jul 5, 2020 at 9:58
0
$\begingroup$

I would like to present my version of code for simulation:

Clear[pick]
pick[n_Integer, switch : (True | False)] := 
    Module[{bonuses = {"Goat", "Goat", "Car"}, samples},
    samples = Which[
        switch == False, RandomChoice[bonuses, n],
        switch == True, Table[MapAt[DeleteCases[#, "Goat", 1, 1] &, TakeDrop[RandomSample[bonuses], 1], 2][[2, 1]], n]
    ];
    Counts[samples]/n // N // KeySort
];

pick[5000, False]

<|"Car" -> 0.3352, "Goat" -> 0.6648|>

pick[5000, True]

<|"Car" -> 0.6694, "Goat" -> 0.3306|>

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.