Skip to content

Instantly share code, notes, and snippets.

@alekrutkowski
Created April 24, 2026 15:16
Show Gist options
  • Select an option

  • Save alekrutkowski/082cde7ca6a2b70c0a3f23f248d7dfa4 to your computer and use it in GitHub Desktop.

Select an option

Save alekrutkowski/082cde7ca6a2b70c0a3f23f248d7dfa4 to your computer and use it in GitHub Desktop.
Monty Hall problem simulation in R

See: https://en.wikipedia.org/wiki/Monty_Hall_problem

You can run it in the browser

# "Pedagogical" version

drawOne <- function(x)
  x[sample.int(length(x), 1)]

remove <- setdiff

game <- function(unused) {
  doors <- 1:3
  winning_door <- drawOne(doors)
  # Round 1:
  choice1 <- drawOne(doors)
  # Round 2:
  revealed_by_monty <- 
    doors |>
    remove(winning_door) |>
    remove(choice1) |>
    drawOne()
  choice2 <- data.frame(
    # Column names = strategies:
    `Stick to the first choice` = 
      choice1,
    `Re-draw randomly from the still closed` =
      doors |>
      remove(revealed_by_monty) |>
      drawOne(),
    `Switch from the first choice` =
      doors |>
      remove(revealed_by_monty) |>
      remove(choice1) |>
      drawOne(),
    check.names=FALSE)
  # Report which strategy won (TRUE) or lost (FALSE)
  choice2 |>
    lapply(\(choice) choice==winning_door) |>
    as.data.frame(check.names=FALSE)
}

# Which strategy is the best one?
set.seed(54321)
message('Calculating...')
percent_of_wins <-
  1:10000 |>
  lapply(game) |>
  do.call(rbind, args = _) |> 
  colMeans() |>
  (\(x) x*100)()
message('Per cent of wins by strategy:')
print(percent_of_wins)
# More concise and faster version

pick <- function(x) x[sample.int(length(x), 1)]

play <- function() {
  doors <- 1:3
  win <- pick(doors)
  first <- pick(doors)
  monty <- pick(setdiff(doors, c(win, first)))
  closed <- setdiff(doors, monty)

  c(
    `Stick to the first choice` = first,
    `Re-draw randomly from the still closed` = pick(closed),
    `Switch from the first choice` = setdiff(closed, first)
  ) == win
}

set.seed(54321)

percent_of_wins <- rowMeans(vapply(seq_len(1e4), \(i) play(), logical(3))) * 100

print(percent_of_wins)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment