**Statistical Research » R**, and kindly contributed to R-bloggers)

I was listening to Science Friday from Sep 6th. and one of the discussions by Ira Flatow was on the well known *Monty Hall Problem*. This problem has been hashed out many times and in fact I was first introduced to the probability aspects of the problem while taking my introductory to statistics class when I was an undergraduate. Though the solution is not initially intuitive it is fairly straight forward and a discussion is available on Wikipedia. I figured I would take a few minutes to play the game a couple million times. This problem is so ubiquitous that Mythbusters did an episode on this back in 2011 where they too kept on playing the game over and over. So there is no need for me to detail it out but here’s a brief simulation showing the probability of winning, if you change doors, is 0.66. The following graphs shows the cumulative results of the first 100 iterations comparing the two difference strategies.

library(e1071) doors = c("A","B","C") nsim=1000000 sim.mat = matrix(NA, nrow=nsim, ncol=5) winning.loc = rdiscrete(nsim,doors, prob=rep(1,length(doors))) first.choice = rdiscrete(nsim,doors, prob=rep(1,length(doors))) sim.mat[,1] = winning.loc sim.mat[,2] = first.choice for(i in 1:length(winning.loc)){ if(winning.loc[i] != first.choice[i]){ ## Contestant has selected the goat. Therefore the host must show the other goat and not the winner. sim.mat[i,3] = doors[ !(doors %in% c(winning.loc[i], first.choice[i])) ] } else { rand.goat.options = doors[ !(doors %in% c(winning.loc[i], first.choice[i])) ] ## Contestant has selected the winner and the host then randomly select one of the other two. rand.goat = rdiscrete(1, rand.goat.options , prob=rep(1,length(rand.goat.options)) ) sim.mat[i,3] = rand.goat } ## Second choice changes sim.mat[i,5] = doors[ !(doors %in% c(sim.mat[i,3], first.choice[i])) ] } ## Second choice is the same as the first choice sim.mat[,4] = sim.mat[,2] ## names = TRUE.VAL, FIRST.CHOICE, GOAT.LOC.SHOWN, SECOND.CHOICE.STAY sim.val.stay = as.numeric( (sim.mat[,4]==sim.mat[,1])*1 ) sim.val.change = as.numeric( (sim.mat[,5]==sim.mat[,1])*1 ) mean(sim.val.stay) mean(sim.val.change) ## Probability by iteration covergence.prob.change = cumsum(sim.val.change)/cumsum(rep(1,length(sim.val.change))) covergence.prob.stay = cumsum(sim.val.stay)/cumsum(rep(1,length(sim.val.stay))) par(mfrow=c(1,2)) plot(head(covergence.prob.stay, 100), pch=16, cex=.5, ylab="Probability", main="Probability of Winning\nKeeping the Same Door") plot(head(covergence.prob.change, 100), pch=16, cex=.5, ylab="Probability", main="Probability of winning\nChanging Doors")

**leave a comment**for the author, please follow the link and comment on his blog:

**Statistical Research » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...