Simulation of Monty Hall game1

Single run of the game

classic_monty <- function() {
  
  prize <- sample(1:3,1)     # Assign the prize
  choice <- sample(1:3,1)    # Pick a door

  # Monty picks a door (not your choice, and not where the prize is)
  monty <- sample(c(1:3)[-c(choice,prize)], 1)
  
  # if the prize is not under your door, you win if you switch
  # if the prize IS under your door, you win if you stick
  return(win = ifelse(prize!=choice, "Switch", "Stick"))
}

Replicate runs

#n <- 2^(1:16)  # 2, 4, 8, 16, ... 65536
n = c(1,seq(10,1000,by=10))
runs <- data.frame(n=numeric(), switch=numeric(), stay=numeric())
for (trials in n) {
  run <- table(replicate(trials, classic_monty()))
  runs <- runs %>%  add_row(n=trials, switch=(sum(run["Switch"]))/trials, stay=(sum(run["Stick"]))/trials)
}
# Handle zero-occurrence trials
runs[is.na(runs)]<-0

Look at the results

head(runs)
##    n    switch      stay
## 1  1 1.0000000 0.0000000
## 2 10 0.8000000 0.2000000
## 3 20 0.7500000 0.2500000
## 4 30 0.7333333 0.2666667
## 5 40 0.7750000 0.2250000
## 6 50 0.6200000 0.3800000

Plot the results

#ggplot(data = runs, aes(x=log(n, base=2), y=switch)) +
ggplot(data = runs, aes(x=n), y=switch) +
  geom_line(aes(y=switch), color="blue") +
  geom_line(aes(y=stay), color="red") +
  ylim(c(0,1)) +
  labs(x="number of trials", y="proportion of wins")

  1. 1↩︎