How to simulate Monty Hall problem in computer?

Question

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 [but the door is not opened], 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?

This is Monty Hall problem, and we know that a switching strategy really does win two out of three times on the average. But is it possible to simulate it in Matlab?

Answer 1

Some pseudocode:

1. Set winswap = 0, winnoswap = 0 and simcount = 0.
2. Increment simcount by 1.
3. Set car = random value from a choice of $\{ 1, 2, 3\}$.
4. Set guess = random value from a choice of $\{1, 2, 3\}$.
5. Set goat = random value from a choice of $\{1, 2, 3\}$ that is not equal to either car or guess.
6. Set swap = random value from a choice of $\{1, 2\}$.
7. If swap = 2, set guess = the first value of $\{1, 2, 3\}$ that is not equal to either goat or guess.
8. If guess = car and swap = 1, increment winnoswap by 1.
   If guess = car and swap = 2, increment winswap by 1.
9. If simcount is below some predetermined value, return to step 2.
10. Compare the values of winswap and winnoswap.

Answer 2

After the interesting Mythbusters episode on the problem, I wrote a bit of Python code to simulate the results for fun. The code can be simplified considerably by breaking it into two trials, one where you always stay with your choice, and the other where you always swap. A further simplification can be had with the following two observations:

1. If you always keep your choice, you only win when you select the car (guess == car)
2. If you always swap, you only win when you do not select the car (guess != car)

This allows one to avoid computing the positions of the goats and the door Monty decides to remove. I am not really familiar with MATLAB, so I will just include the Python code below, but I hope it gives the basic idea for porting to your language of choice.

import random
import sys

def main():
  if len(sys.argv) < 2:
    print "usage: monty.py trials"
    sys.exit(1)

  t = int(sys.argv[1])

  # Stick with door
  stick_winners = 0
  for i in xrange(t):
    car = random.randint(1, 3)
    choice = random.randint(1, 3)
    if car == choice:
      stick_winners += 1
  print "(stick) won:", stick_winners, "out of", t, float(stick_winners) / float(t)

  # Switch door
  switch_winners = 0
  for i in xrange(t):
    car = random.randint(1, 3)
    choice = random.randint(1, 3)
    # If car != choice, open a goat door and switch, win the car
    if car != choice:
      switch_winners += 1
  print "(switch) won:", switch_winners, "out of", t, float(switch_winners) / float(t)

if __name__ == "__main__":
  main()

Answer 3

Python code n doors k revealed

Here another way in python to compare simulations with the general result in

Monty hall problem extended.

If extended to $n$ doors, $k$ doors are revealed by the host after the first choice, the probability of winning with switching and without switching changes slightly.

Here is the general formula for the Monty Hall problem with $n$ doors and $k$ revealed doors:

• If the player does not switch :

$$P(win) = 1/n$$

• If the player switches and the number of revealed doors is $k=n-2$ :

$$P(win) =\frac{n-1}{n}$$

• And in general for any number k of revealed doors:

$$0\leq k\leq n-2$$

$$P(win)=\frac{n-1}{n}\cdot \frac{1}{n-k-1} = \frac{1}{n} \cdot \frac{n-1}{n-k-1} \geq \frac{1}{n}$$

n=15
k=13# To simulate the original case k should be k=n-2 where p(win)=(n-1)/n
p_win_sw=((n-1)/n)*(1/(n-k-1))
p_win_no_sw=1/n


print('Total number of doors: ',n,'Number of revealed doors: ',k)
print("**********************************************************")
print('Probability of winning with switching: ',p_win_sw," = ", Fraction.from_float(p_win_sw).limit_denominator())
print('Probability of winning with NO switching:',p_win_no_sw," = ", Fraction.from_float(p_win_no_sw).limit_denominator())





import random
from fractions import Fraction

def simulate_game(switch, num_doors, k_revealed):
    # Define the doors as a list of goats and a car
    doors = ['goat'] * (num_doors - 1) + ['car']

    # Shuffle the doors randomly
    random.shuffle(doors)

    # User chooses a door
    first_choice = random.randint(0, num_doors - 1)

    # Host reveals k_revealed doors with goats
    revealed_doors = []
    for i in range(num_doors):
        if i != first_choice and doors[i] == 'goat' and len(revealed_doors) < k_revealed:
            revealed_doors.append(i)

    # User decides whether to switch or keep their original choice
    if switch:
        # If switching, the user must choose a door that was not already revealed
        remaining_doors = [i for i in range(num_doors) if i not in revealed_doors and i != first_choice]
        second_choice = random.choice(remaining_doors)
        return doors[second_choice] == 'car'
    else:
        # If not switching, the user keeps their original choice
        return doors[first_choice] == 'car'
    

# Simulate the game for a specified number of times
num_simulations = 100000
num_doors = n
k_revealed = k

num_wins_with_switching = 0
num_wins_without_switching = 0
for i in range(num_simulations):
    if simulate_game(True, num_doors, k_revealed):
        num_wins_with_switching += 1
    if simulate_game(False, num_doors, k_revealed):
        num_wins_without_switching += 1

# Print the results
print("**********************************************************")
print("Simulated probability of winning with switching: ", num_wins_with_switching / num_simulations)
print("Simulated probability of winning without switching: ", num_wins_without_switching / num_simulations)

Answer 4

Code I created based off of Peter Phipps pseudocode for school assignment (very helpful by the way):

n = int(input("Enter the number of iterations: "))

winswap = 0.0
winnoswap = 0.0
simcount = 0.0

for i in range(n):
    simcount = simcount + 1
    picked_door = random.randint(1,3)
    car = random.randint(1,3)
    g = random.randint(1,3)

    completed = False 
    while not completed:
        if g != car and g != picked_door:
            completed = True
        else:
            g = random.randint(1,3)      

    swap = random.randint(1,2)


    if swap == 2:
        if g != 1 and picked_door != 1:
            picked_door = 1
        elif g != 2 and picked_door != 2:
            picked_door = 2
        else:
            picked_door = 3 



    if picked_door == car and swap == 1:
        winnoswap = winnoswap + 1
    if picked_door == car and swap == 2:
        winswap = winswap + 1

print "The amount of no swap wins is", winnoswap
print "The amount of swap wins is", winswap

total_wins = winnoswap + winswap

percentswap = (winswap/total_wins)*100
percentnoswap = (winnoswap/total_wins)*100

print "Percentage wins of swapping is", percentswap,"%"
print "Percentage wins of not swapping is", percentnoswap,"%"

Answer 5

I have created a general Monty Hall simulation here.

Answer 6

Here is my MATLAB solution:

    clear all
    clc

trial=input('Enter the number of trials: ');
win_swap = 0; win_noswap = 0; counter = 0;
for i = 1:trial 
  car = randi(3);
  guess = randi(3);
  goat = randi(3);

  if(car~=goat)
  counter = counter+1;

    if (car == guess)
    win_noswap = win_noswap+1;

    elseif (car ~= guess)
    win_swap = win_swap +1;

    end % car = guess end

  end % car ~= goat end

end %for end
  R = sprintf('IN %d CONVENIENT TRIALS:\n',counter);
  disp(R)
  NS =sprintf('The amount of no swap wins is %d,',win_noswap);
  disp(NS)
  S = sprintf('The amount of swap wins is %d,',win_swap);
  disp(S)
%Percentage  
total_wins = win_noswap + win_swap;
percent_swap = (win_swap/counter)*100;
percent_noswap = (win_noswap/counter)*100;

  PNS=sprintf('Percentage wins of not swapping is %% %d',percent_noswap);
  disp(PNS)
  PS=sprintf('Percentage wins of swapping is %% %d',percent_swap);
  disp(PS)
  disp('ALWAYS SWAP THE DOOR!');
  %Pie Chart
  percent_pie = [percent_swap percent_noswap];
  pie(percent_pie)

and a sample output can be found here.