Monty Hall Problem#

FIZ371 - Scientific & Technical Calculations | 08/03/2021

Emre S. Tasci emre.tasci@hacettepe.edu.tr

import numpy as np
import random

N = 100000

#np.random.seed(371)

3 doors, 1 door is opened#

tries = np.random.randint(1,4,(N,5))
tries[:,2:5] = 0

# 0th col: Door with the prize
# 1st col: Door that the contestant has chosen
# 2nd col: Door that Monty opens
# 3rd col: 1 if the user wins by sticking, 0 otherwise
# 4th col: 1 if the user wins by switching, 0 otherwise
# tries
filt = tries[:,0]==tries[:,1]
#tries[filt]
for i in np.array(range(N))[filt]:
    possible = np.array([1,2,3])
    possible = possible[possible != tries[i,0]]
    tries[i,2] = random.choice(possible)

for i in np.array(range(N))[filt]:
    possible = np.array([1,2,3])
    possible = possible[possible != tries[i,0]]
    tries[i,3] = 1
nfilt = np.invert(filt)
#tries[nfilt]
for i in np.array(range(N))[nfilt]:
    possible = np.array([1,2,3])
    possible = possible[possible != tries[i,0]]
    possible = possible[possible != tries[i,1]]
    tries[i,2] = possible[0]
    tries[i,4] = 1
tries
array([[1, 1, 3, 1, 0],
       [1, 2, 3, 0, 1],
       [1, 1, 3, 1, 0],
       ...,
       [3, 1, 2, 0, 1],
       [1, 1, 3, 1, 0],
       [1, 1, 2, 1, 0]])
print("# of total wins if sticks  : {:} [prob: {:.5f}]"\
      .format(np.sum(tries[:,3]),np.sum(tries[:,3])/N))
print("# of total wins if switches: {:} [prob: {:.5f}]"\
      .format(np.sum(tries[:,4]),np.sum(tries[:,4])/N))
# of total wins if sticks  : 33075 [prob: 0.33075]
# of total wins if switches: 66925 [prob: 0.66925]

4 doors, 2 are opened simultaneously#

tries4 = np.random.randint(1,5,(N,6))
tries4[:,2:] = 0

# 0th col: Door with the prize
# 1st col: Door that the contestant has chosen
# 2nd col: Door that Monty opens
# 3rd col: Door that Monty opens
# 4th col: 1 if the user wins by sticking, 0 otherwise
# 5th col: 1 if the user wins by switching, 0 otherwise
# tries4
filt4 = tries4[:,0] == tries4[:,1]

tries4[filt4,4] = 1
tries4[np.invert(filt4),5] = 1
a = np.array([1,2,3,4])
for i in range(N):
    b = a[a != tries4[i,0]]
    b = b[b != tries4[i,1]]
    b = np.sort(np.random.choice(b,size=2,replace = False))
    #print(b);
    tries4[i,2:4] = b
tries4
array([[1, 4, 2, 3, 0, 1],
       [4, 3, 1, 2, 0, 1],
       [3, 4, 1, 2, 0, 1],
       ...,
       [3, 4, 1, 2, 0, 1],
       [3, 1, 2, 4, 0, 1],
       [1, 3, 2, 4, 0, 1]])
print("# of total wins if sticks  : {:} [prob: {:.5f}]"\
      .format(np.sum(tries4[:,4]),np.sum(tries4[:,4])/N))
print("# of total wins if switches: {:} [prob: {:.5f}]"\
      .format(np.sum(tries4[:,5]),np.sum(tries4[:,5])/N))
# of total wins if sticks  : 25147 [prob: 0.25147]
# of total wins if switches: 74853 [prob: 0.74853]

4 doors, a door is opened at each stage#

N=10000
tries41 = np.random.randint(1,5,(N,12))
tries41[:,2:] = 0

# 0th  col: Door with the prize
# 1st  col: Door that the contestant has chosen first (stick - *)
# 2nd  col: Door that Monty opens first
# 3rd  col: Door that the contestant switches on the first ask (switch - *)
# 4th  col: Door that Monty opens second (given the contestant has stuck on the 1st)
# 5th  col: Door that Monty opens second (given the contestant has switched on the 1st)
# 6th  col: Door that the contestant switches on the second ask (stick - switch)
# 7th  col: Door that the contestant switches on the second ask (switch - switch)
# 8th  col: 1 if the user has stick-stick win
# 9th  col: 1 if the user has stick-switch win
# 10th col: 1 if the user has switch-stick win
# 11th col: 1 if the user has switch-switch win

tries41[:10]
array([[2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
a = np.array([1,2,3,4])

for i in range(N):
    
    # The 1st door Monty opens
    b = a[a != tries41[i,0]] # Can't open the one with the prize
    b = b[b != tries41[i,1]] # Can't open the contestant's pick
    tries41[i,2] = np.random.choice(b)
    
    # The door the contestant switches (1st)
    b = a[a != tries41[i,1]] # Can't switch to the already picked
    b = b[b != tries41[i,2]] # Can't switch to the opened door
    tries41[i,3] = np.random.choice(b)

    # The 2nd door Monty opens (given the contestans has stuck)
    b = a[a != tries41[i,0]] # Can't open the one with the prize
    b = b[b != tries41[i,1]] # Can't open the contestant's pick
    b = b[b != tries41[i,2]] # Can't open the already opened
    tries41[i,4] = np.random.choice(b)    
    
    # The 2nd door Monty opens (given the contestans has switched)
    b = a[a != tries41[i,0]] # Can't open the one with the prize
    b = b[b != tries41[i,3]] # Can't open the contestant's pick
    b = b[b != tries41[i,2]] # Can't open the already opened
    tries41[i,5] = np.random.choice(b)    
    
    # The door the contestant switches (after having stuck)
    b = a[a != tries41[i,2]] # Can't switch to the opened door
    b = b[b != tries41[i,4]] # Can't switch to the opened door
    b = b[b != tries41[i,1]] # Can't switch to the already picked
    tries41[i,6] = np.random.choice(b)
    
    # The door the contestant switches (after having switched once)
    b = a[a != tries41[i,2]] # Can't switch to the opened door
    b = b[b != tries41[i,5]] # Can't switch to the opened door
    b = b[b != tries41[i,3]] # Can't switch to the already picked
    tries41[i,7] = np.random.choice(b)

tries41[:10]
array([[2, 1, 3, 4, 4, 1, 2, 2, 0, 0, 0, 0],
       [3, 1, 2, 3, 4, 1, 3, 4, 0, 0, 0, 0],
       [4, 1, 3, 2, 2, 1, 4, 4, 0, 0, 0, 0],
       [4, 3, 2, 1, 1, 3, 4, 4, 0, 0, 0, 0],
       [4, 4, 3, 2, 2, 1, 1, 4, 0, 0, 0, 0],
       [4, 4, 3, 1, 1, 2, 2, 4, 0, 0, 0, 0],
       [1, 4, 2, 3, 3, 4, 1, 1, 0, 0, 0, 0],
       [1, 1, 4, 3, 3, 2, 2, 1, 0, 0, 0, 0],
       [2, 4, 3, 1, 1, 4, 2, 2, 0, 0, 0, 0],
       [2, 3, 1, 2, 4, 3, 2, 4, 0, 0, 0, 0]])
# stick - stick wins
tries41[tries41[:,0] == tries41[:,1],8] = 1

# stick - switch wins
tries41[tries41[:,0] == tries41[:,6],9] = 1

# switch - stick wins
# If they first switch, then stick, 
# treat the 3rd column as the final choice
tries41[tries41[:,0] == tries41[:,3],10] = 1

# switch - switch wins
tries41[tries41[:,0] == tries41[:,7],11] = 1

tries41[:10]
array([[2, 1, 3, 4, 4, 1, 2, 2, 0, 1, 0, 1],
       [3, 1, 2, 3, 4, 1, 3, 4, 0, 1, 1, 0],
       [4, 1, 3, 2, 2, 1, 4, 4, 0, 1, 0, 1],
       [4, 3, 2, 1, 1, 3, 4, 4, 0, 1, 0, 1],
       [4, 4, 3, 2, 2, 1, 1, 4, 1, 0, 0, 1],
       [4, 4, 3, 1, 1, 2, 2, 4, 1, 0, 0, 1],
       [1, 4, 2, 3, 3, 4, 1, 1, 0, 1, 0, 1],
       [1, 1, 4, 3, 3, 2, 2, 1, 1, 0, 0, 1],
       [2, 4, 3, 1, 1, 4, 2, 2, 0, 1, 0, 1],
       [2, 3, 1, 2, 4, 3, 2, 4, 0, 1, 1, 0]])
print("# of total wins for stick-stick   : {:} [prob: {:.5f}]"\
      .format(np.sum(tries41[:,8]),np.sum(tries41[:,8])/N))
print("# of total wins for stick-switch  : {:} [prob: {:.5f}]"\
      .format(np.sum(tries41[:,9]),np.sum(tries41[:,9])/N))
print("# of total wins for switch-stick  : {:} [prob: {:.5f}]"\
      .format(np.sum(tries41[:,10]),np.sum(tries41[:,10])/N))
print("# of total wins for switch-switch : {:} [prob: {:.5f}]"\
      .format(np.sum(tries41[:,11]),np.sum(tries41[:,11])/N))
# of total wins for stick-stick   : 2468 [prob: 0.24680]
# of total wins for stick-switch  : 7532 [prob: 0.75320]
# of total wins for switch-stick  : 3813 [prob: 0.38130]
# of total wins for switch-switch : 6187 [prob: 0.61870]
tries41[:10]
array([[2, 1, 3, 4, 4, 1, 2, 2, 0, 1, 0, 1],
       [3, 1, 2, 3, 4, 1, 3, 4, 0, 1, 1, 0],
       [4, 1, 3, 2, 2, 1, 4, 4, 0, 1, 0, 1],
       [4, 3, 2, 1, 1, 3, 4, 4, 0, 1, 0, 1],
       [4, 4, 3, 2, 2, 1, 1, 4, 1, 0, 0, 1],
       [4, 4, 3, 1, 1, 2, 2, 4, 1, 0, 0, 1],
       [1, 4, 2, 3, 3, 4, 1, 1, 0, 1, 0, 1],
       [1, 1, 4, 3, 3, 2, 2, 1, 1, 0, 0, 1],
       [2, 4, 3, 1, 1, 4, 2, 2, 0, 1, 0, 1],
       [2, 3, 1, 2, 4, 3, 2, 4, 0, 1, 1, 0]])
ss = tries41[:,[0,1,6,3,7]] 
ss[:10] # pr, st-st, st-sw, sw-st, sw-sw
array([[2, 1, 2, 4, 2],
       [3, 1, 3, 3, 4],
       [4, 1, 4, 2, 4],
       [4, 3, 4, 1, 4],
       [4, 4, 1, 2, 4],
       [4, 4, 2, 1, 4],
       [1, 4, 1, 3, 1],
       [1, 1, 2, 3, 1],
       [2, 4, 2, 1, 2],
       [2, 3, 2, 2, 4]])