Simulate the Monty Hall problem

montyHall.R

# function runs the monty hall problem 'n' times and returns the proportions of 'wins' under sticking and switching strategies

montyHall <- function (n = 1e5) {
  prize <- sample(1:3, n, replace = T)
  firstChoice <- sample(1:3, n, replace = T)
  reveal <- matrix(rep.int(1L, n * 3), ncol = 3)
  reveal[cbind(rep(seq_len(n), 2), c(prize, firstChoice))] <- 0L
  reveal <- max.col(reveal, ties.method = 'random')
  switchTo <- 6L - (firstChoice + reveal)
  c(stick = mean(firstChoice == prize), switch = mean(switchTo == prize))
}

# 100,000 simulations

montyHall(1e5)

# let's try 100 million

t0 <- proc.time()
montyHall(1e8)
proc.time() - t0

# repeat simulations of size 'sampleSize' multiple ('nSamples') times

plotMontyHall <- function (sampleSize = 20, nSamples = 1e4) {
  x <- sapply(rep.int(sampleSize, nSamples), montyHall)
  hist(x[1,], xlim = c(0, 1), xlab = NA, main = "Stick", col = 'seagreen')
  hist(x[2,], xlim = c(0, 1), xlab = NA, main = "Switch", col = 'seagreen')
}

layout(matrix(1:6, ncol = 3, byrow = F))
plotMontyHall(20)
plotMontyHall(100)
plotMontyHall(1000)