A Monty Simulation

September 10, 2013

(This article was first published on 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.

Monty Hall Simulation

doors = c("A","B","C")

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]
sim.val.stay = as.numeric( (sim.mat[,4]==sim.mat[,1])*1 )
sim.val.change = as.numeric( (sim.mat[,5]==sim.mat[,1])*1 )

## 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)))

plot(head(covergence.prob.stay, 100), pch=16, cex=.5,
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")

To leave a comment for the author, please follow the link and comment on their blog: Statistical Research » R.

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

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.


Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)