Increase efficiency for an R simulator of the Monty Hall Puzzle

The Monty Hall Problem is a simple puzzle involving probability that even stumps professionals in careers dealing with some heavy-duty math. Here's the basic problem:

Suppose you're on a game show, and you're given the choice of three doors: Behind one door is a car; behind the others, goats. You pick a door, say No. 1, and the host, who knows what's behind the doors, opens another door, say No. 3, which has a goat. He then says to you, "Do you want to pick door No. 2?" Is it to your advantage to switch your choice?

You can find numerous explanations of the solution here: http://en.wikipedia.org/wiki/Monty_Hall_problem Goal of my simulation: Prove that a switching strategy will win you the car 2/3 of the time.

I got curious and wanted to write a little function that simulates the problem many times and returns the proportion of wins if you switched and the proportion of wins if you stayed with your first choice. The function then plots the cumulative wins.

This function takes a long time to run once I get to about 10,000 simulations. I know I don't need this many simulations to prove this but I'd love to hear some ideas on how to make it more efficient. Thanks for your feedback!

Monty_Hall=function(repetitions){
doors=c('A','B','C')
stay_wins=0
switch_wins=0
series=data.frame(sim_num=seq(repetitions),cum_sum_stay=replicate(repetitions,0),cum_sum_switch=replicate(repetitions,0))
for(i in seq(repetitions)){
winning_door=sample(doors,1)
contestant_chooses=sample(doors,1)
if(contestant_chooses==winning_door)
stay_wins=stay_wins+1
else
switch_wins=switch_wins+1
series[i,'cum_sum_stay']=stay_wins
series[i,'cum_sum_switch']=switch_wins
}
plot(series$sim_num,series$cum_sum_switch,col=2,ylab='Cumulative # of wins',
xlab='Simulation #',main=sprintf('%d Simulations of the Monty Hall Paradox',repetitions),type='l')
lines(series$sim_num,series$cum_sum_stay,col=4)
legend('topleft',legend=c('Cumulative wins from switching',
'Cumulative wins from staying'),col=c(2,4),lty=1)
result=list(series=series,stay_wins=stay_wins,switch_wins=switch_wins,
proportion_stay_wins=stay_wins/repetitions,
proportion_switch_wins=switch_wins/repetitions)
return(result)
}
#Theory predicts that it is to the contestant's advantage if he
#switches his choice to the other door.  This function simulates the game
#many times, and shows you the proportion of games in which staying or
#switching would win the car.  It also plots the cumulative wins for each strategy.
Monty_Hall(100)
• parallel computing is a solution and it's been done in R with this very problem. See this post: librestats.com/2012/03/15/… – Tyler Rinker Oct 15 '12 at 21:53
• For extra credit, enhance your code to allow N doors, of which M<N have prizes. Amaze and entertain your friends at parties! – Carl Witthoft Oct 16 '12 at 11:25

R has vectorized operators and using those is much more efficient in this case:

Monty_Hall2 <- function(repetitions){
doors <- c('A','B','C')
winning_door <- sample(doors, repetitions, replace=TRUE)
contestant_choose <- sample(doors, repetitions, replace=TRUE)
stay_wins <- winning_door == contestant_choose
switch_wins <- winning_door != contestant_choose
series <- data.frame(sim_num = seq_len(repetitions),
cum_sum_stay = cumsum(stay_wins),
cum_sum_switch = cumsum(switch_wins))
plot(series$sim_num,series$cum_sum_switch,col=2,ylab='Cumulative # of wins',
xlab='Simulation #',main=sprintf('%d Simulations of the Monty Hall Paradox',repetitions),type='l')
lines(series$sim_num,series$cum_sum_stay,col=4)
legend('topleft',legend=c('Cumulative wins from switching',
'Cumulative wins from staying'),col=c(2,4),lty=1)
list(series = series,
stay_wins = series[repetitions,"cum_stay_wins"],
switch_wins = series[repetitions,"cum_switch_wins"],
proportion_stay_wins = series[repetitions,"cum_stay_wins"]/repetitions,
proportion_switch_wins = series[repetitions,"cum_switch_wins"]/repetitions)
}

That is, generate all the answers and guesses, check for each one (vectorized) which (staying or switching) would be a win, and then accumulate those results. The rest is just to put the data in the same format that your function already had.

> system.time(tmp <- Monty_Hall(10000))
user  system elapsed
6.04    0.00    6.04
> system.time(tmp <- Monty_Hall2(10000))
user  system elapsed
0.01    0.02    0.03

So approximately 1/600 of the time. Or, since I didn't do proper benchmarking, say between 2 and 3 orders of magnitude quicker.