Monty Hall Problem and probability, simulated in R, Bayes Nets, bayesian nets.

montyHall.R

library(ggplot2)

montyHall <- function(maxDoors = 3, switch = T, tries = 10) {
  plays <- data.frame()

    for (i in 1:tries) {
    # Create the doors.
    doors <- 1:maxDoors

    # Choose a winning door.
    winningDoor <- sample(doors, 1)

    # Choose the player door.
    playerDoor <- sample(doors, 1)
    o <- playerDoor
    
    if (playerDoor == winningDoor) {
      # The player has the winning door, so Monty reveals one of the other doors, containing a goat.
      montyDoor <- sample(doors[-playerDoor], 1)
    } else {
      # The player has a goat, so Monty reveals one of the other doors, not including the winning door.
      montyDoor <- unlist(sample(as.list(doors[-c(playerDoor, winningDoor)], 1)))
    }

    # Let the player switch doors.
    if (switch == T) {
      playerDoor <- unlist(sample(as.list(doors[-c(playerDoor, montyDoor)], 1)))
    }

    plays <- rbind(plays, data.frame(o, playerDoor, winningDoor, switch, playerDoor == winningDoor))
  }
  
  colnames(plays) <- c('originalDoor', 'playerDoor', 'winningDoor', 'switch', 'win')
  plays$win <- as.logical(plays$win)
  plays
}

switch <- data.frame()
for (count in seq(from=100, to=1000, by=100)) {
  result <- montyHall(switch = T, tries = count)
  switch <- rbind(switch, c(count, sum(result$win) / count))
}
names(switch) <- c('count', 'wins')

stay <- data.frame()
for (count in seq(from=100, to=1000, by=100)) {
  result <- montyHall(switch = F, tries = count)
  stay <- rbind(stay, c(count, sum(result$win) / count))
}
names(stay) <- c('count', 'wins')

p = ggplot() +
  geom_line(data = switch, aes(x = count, y = wins, color='Switch')) +
  geom_line(data = stay, aes(x = count, y = wins, color='Stay')) +
  xlab('# of Games') +
  ylab('Win Percentage') +
  theme(legend.position = 'bottom') +
  labs(colour='Type') +
  scale_color_manual(values = c('red', 'blue'))

print(p)

result1 <- montyHall(switch = F, tries = 500)
result2 <- montyHall(switch = T, tries = 500)

print(sum(result1$win) / nrow(result1))
print(sum(result2$win) / nrow(result2))