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)
R has vectorized operators and using those is much more efficient in this case:

Monty_Hall2 &lt;- function(repetitions){
    doors &lt;- c('A','B','C')
    winning_door &lt;- sample(doors, repetitions, replace=TRUE)
    contestant_choose &lt;- sample(doors, repetitions, replace=TRUE)
    stay_wins &lt;- winning_door == contestant_choose
    switch_wins &lt;- winning_door != contestant_choose
    series &lt;- 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.

&gt; system.time(tmp &lt;- Monty_Hall(10000))
   user  system elapsed 
   6.04    0.00    6.04 
&gt; system.time(tmp &lt;- 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.