Single Trial
monty_hall <- function() {
key <- sample(1:3, size = 1)
goat <- setdiff(1:3, key)
contestant <- sample(1:3, size = 1)
monty <- ifelse(contestant == key, sample(goat, size = 1), setdiff(goat, contestant))
switch <- setdiff(1:3, c(contestant, monty))
result <- ifelse(switch == key, "Switching wins", "Staying wins")
# result
c("Key" = key, "Contestant" = contestant, "Monty" = monty, "Switch" = switch, "Result" = result)
}
monty_hall()
## Key Contestant Monty Switch Result
## "2" "3" "1" "2" "Switching wins"
N trials
N <- 30
monty_result <- replicate(N, monty_hall())
t(monty_result)
## Key Contestant Monty Switch Result
## [1,] "2" "3" "1" "2" "Switching wins"
## [2,] "1" "1" "2" "3" "Staying wins"
## [3,] "2" "2" "3" "1" "Staying wins"
## [4,] "1" "2" "3" "1" "Switching wins"
## [5,] "1" "1" "2" "3" "Staying wins"
## [6,] "3" "2" "1" "3" "Switching wins"
## [7,] "1" "1" "3" "2" "Staying wins"
## [8,] "3" "2" "1" "3" "Switching wins"
## [9,] "1" "3" "2" "1" "Switching wins"
## [10,] "2" "2" "3" "1" "Staying wins"
## [11,] "2" "3" "1" "2" "Switching wins"
## [12,] "3" "1" "2" "3" "Switching wins"
## [13,] "2" "3" "1" "2" "Switching wins"
## [14,] "2" "2" "1" "3" "Staying wins"
## [15,] "1" "2" "3" "1" "Switching wins"
## [16,] "3" "3" "1" "2" "Staying wins"
## [17,] "1" "3" "2" "1" "Switching wins"
## [18,] "3" "2" "1" "3" "Switching wins"
## [19,] "1" "3" "2" "1" "Switching wins"
## [20,] "2" "1" "3" "2" "Switching wins"
## [21,] "1" "1" "2" "3" "Staying wins"
## [22,] "1" "2" "3" "1" "Switching wins"
## [23,] "1" "3" "2" "1" "Switching wins"
## [24,] "3" "2" "1" "3" "Switching wins"
## [25,] "1" "1" "2" "3" "Staying wins"
## [26,] "2" "1" "3" "2" "Switching wins"
## [27,] "3" "2" "1" "3" "Switching wins"
## [28,] "2" "2" "1" "3" "Staying wins"
## [29,] "2" "1" "3" "2" "Switching wins"
## [30,] "1" "3" "2" "1" "Switching wins"
table(unlist(monty_result[5, ]))
##
## Staying wins Switching wins
## 10 20
sum(monty_result[5, ] == "Switching wins")/N
## [1] 0.6666667
cumsum(monty_result[5, ] == "Switching wins")
## [1] 1 1 1 2 2 3 3 4 5 5 6 7 8 8 9 9 10 11 12 13 13 14 15 16 16 17 18 18 19 20
cumsum(monty_result[5, ] == "Staying wins")
## [1] 0 1 2 2 3 3 4 4 4 5 5 5 5 6 6 7 7 7 7 7 8 8 8 8 9 9 9 10 10 10
y_switch <- cumsum(monty_result[5, ] == "Switching wins")
y_stay <- cumsum(monty_result[5, ] == "Staying wins")
Plot
monty_plot <- function(N) {
monty_result <- replicate(N, monty_hall())
y_switch <- cumsum(monty_result[5, ] == "Switching wins")
y_stay <- cumsum(monty_result[5, ] == "Staying wins")
p_wins <- sum(monty_result[5, ] == "Switching wins")/N
plot(x = 1:N,
y = y_switch / N,
pch = 23,
col = "red",
bg = "red",
ylim = c(0, 4/5),
xlab = "Number of Trials",
ylab = "Proportion of Wins",
yaxt = "n",
cex = 0.7)
axis(side = 2,
at = c(0, 1/3, 2/3, 3/4),
labels = c("0", "1/3", "2/3", "3/4"), las = 2)
points(x = 1:N,
y = y_stay / N,
pch = 22,
col = "blue",
bg = "blue",
cex = 0.7)
abline(h = c(1/3, 2/3), lty = 3)
title(main = "Monty Hall Simulation")
legend("topleft",
inset = 0.05,
legend = c("Switching wins", "Staying wins"),
pch = c(23, 22),
col = c("red", "blue"),
pt.bg = c("red", "blue"))
text(x = N / 5, y = 1 / 2,
labels = paste0("P(Switching wins) = ", format(p_wins, digits = 2, nsmall = 2)))
p_wins
}
monty_plot(30)
## [1] 0.5666667
