Increase efficiency for an R simulator of the Monty Hall Puzzle
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.