R Notebook

The so-called Monty Hall or Let’s Make a Deal problem has caused much consternation over the years. It is named for an old television program. A contestant is presented with three doors. Behind one door is a fabulous prize; behind the other two doors are virtually worthless prizes. The contestant chooses a door. The host of the show, Monty Hall, then opens one of the remaining two doors, revealing one of the worthless prizes. Because Monty is the host, he knows which doors conceal the worthless prizes and always chooses one of them to reveal, but never the door chosen by the contestant. Then the contestant is offered the choice of keeping what is behind her original door or trading for what is behind the remaining unopened door. What should she do?

There are two popular answers. • There are two unopened doors, they are equally likely to conceal the fabulous prize, so it doesn’t matter which one she chooses. • She had a 1/3 probability of choosing the right door initially, a 2/3 chance of getting the prize if she trades, so she should trade.

Solving the problem by simulation: Let’s calculate the probability of winning if you don’t change your inital choice after Monte reveals one door without a prize.

nSims<-10000
prizeMat<-rmultinom(n=nSims, size=1, prob=c(rep(1,3))) # 1 indicates prize, 0 no prize.
doorChoseMat<-rmultinom(n=nSims, size=1, prob=c(rep(1,3)))
winnerMat<-prizeMat*doorChoseMat
sum(apply(winnerMat,2,sum))/nSims 

[1] 0.3394

~1/3 as expected

Next, Let’s next calculate the probability of winning if you do change your inital choice after Monte reveals one door without prize. Note that we know this has to equal 2/3 so that the total probability sums to 1.

nSims<- 10000
prizeMat<-rmultinom(n=nSims, size=1, prob=c(rep(1,3))) # 1 indicates prize, 0 no prize.
doorChoseMat<-rmultinom(n=nSims, size=1, prob=c(rep(1,3)))
winCount<-0
for(i in 1:nSims){ 
  indexChoosen<-which(doorChoseMat[,i]==1, arr.ind = TRUE) # index of the door choosen by contestant
  index0s<-which(prizeMat[,i]==0, arr.ind = TRUE) # positions that contain 0's
  setMonteCanReveal<-setdiff(index0s,indexChoosen)
  doorMonteReveals<-ifelse(length(setMonteCanReveal)>1,sample(x=setMonteCanReveal,size=1),setMonteCanReveal)
  doorChoosen<-setdiff(c(1,2,3),c(doorMonteReveals,indexChoosen))
  if(prizeMat[doorChoosen,i]==1){winCount<-winCount+1}
  # browser()
}
winCount/nSims

[1] 0.6734

Probability is ~2/3 if you change your door