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 show the cumulative results of the first 100 iterations comparing the two different 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")

A Month Hall Problem Simulation. http://t.co/LXNkN8fbOn

#simulation #r #rstats