public
Created

Monty Hall Monte Carlo Code

  • Download Gist
MontyMonteAll.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
#### Monty Hall Monte Carlo
#### Rob Mealey
library(ggplot2)
library(RColorBrewer)
library(reshape2)
 
### Function: Run simulation n times and plot results in stacked bar histograms
montyMonte <- function(n,titleSize=7,legendTitle=5,ytextSize=5,xtextSize=5){
Picks <- c()
Opens <- c()
WinningDoor <- c()
WinsIfSwitch <- c()
PlacedDf <- matrix(nrow=n, ncol=3)
OpensDf <- matrix(nrow=n, ncol=3)
PicksDf <- matrix(nrow=n, ncol=3)
Doors <- c('Door 1', 'Door 2', 'Door 3')
colnames(PlacedDf) <- Doors
colnames(PicksDf) <- Doors
colnames(OpensDf) <- Doors
########## Simulation Loop #########
for (i in seq(n)){
## 1. Randomly place prize behind one of three doors
PlacePrize <- c(1,0,0)[sample(1:3,3)]
## 2. Randomly pick one of three doors
YouPick <- Doors[sample(1:3,1)]
## 3. Monty either randomly opens one of the two doors left over if
## you happen to pick the correct door or picks the only door left
## if you pick one of two incorrect doors
MontyOpens <- ifelse(PlacePrize[Doors==YouPick]==1,
Doors[!Doors%in%YouPick][sample(1:2,1)],
Doors[(!Doors%in%c(YouPick,Doors[PlacePrize==1]))])
PrizeIsBehind <- Doors[PlacePrize==1]
## 4. If the prize is behind the leftover door, you win if you switch.
## Else you win if you stick on your original choice.
WinIfSwitch <- ifelse(PlacePrize[!Doors%in%c(YouPick,MontyOpens)]==1,1,0)
Picks <- c(Picks, YouPick)
Opens <- c(Opens, MontyOpens)
WinningDoor <- c(WinningDoor, PrizeIsBehind)
WinsIfSwitch <- c(WinsIfSwitch, WinIfSwitch)
### Write results to data frames
PlacedDf[i,] <- PlacePrize
PicksDf[i,YouPick] <- 2
OpensDf[i,MontyOpens] <- 3}
########## End Simulation Loop #########
WinsIfSwitches <- ifelse(WinsIfSwitch==1,
'Switch Door = Win','Switch Door = Lose')
Games <- data.frame(Picks, Opens, WinningDoor, WinsIfSwitches)
Wins <- sum(WinsIfSwitch)/n
Games <- melt(Games,measure.vars=c('Picks', 'Opens',
'WinningDoor', 'WinsIfSwitches'))
Games$variable <- ordered(Games$variable, levels=c('WinsIfSwitches',
'WinningDoor',
'Opens','Picks'))
PicksDf[is.na(PicksDf)] <- 0
OpensDf[is.na(OpensDf)] <- 0
ResultsDf <- rbind(
data.frame('Type'=rep('Placed',n*3),melt(PlacedDf,measure.vars=Doors)),
data.frame('Type'=rep('Picked',n*3),melt(PicksDf,measure.vars=Doors)),
data.frame('Type'=rep('Opens',n*3),melt(OpensDf,measure.vars=Doors)))
colnames(ResultsDf) <- c('Type','Trial','Door','value')
# Plot stacked bar histograms of your picks, monty's opens, winning doors
# and win if switch
ggplot(Games, aes(x=variable, fill=factor(value)))
last_plot() + geom_histogram()
last_plot() + scale_x_discrete(labels=rev(c('Your Picks',"Monty's Opens",
'Winning Door','Switch=Win/Lose')))
last_plot() + scale_fill_brewer(type='qual',palette=6) + xlab('') + ylab('')
last_plot() + theme_bw() + coord_flip()
last_plot() + opts(title =
paste('Monty Hall Monte Carlo Total Simulation Results, N = ',n,', Pct Switches Win = ',Wins,sep=''),
legend.position='bottom',legend.title=theme_blank())
last_plot() + opts(plot.title = theme_text(size=titleSize),
legend.text = theme_text(size=legendTitle),
axis.text.y = theme_text(size=ytextSize),
axis.text.x = theme_text(size=xtextSize))
ggsave(paste('MontyMonteHistograms',n,'.png',sep=''),width=5, height=3)
WinsIfSwitches <- factor(WinsIfSwitches)
ResultsDf$lineTypes <- ordered(rep(WinsIfSwitches,3*3),
levels=rev(levels(WinsIfSwitches)))
ResultsDf$Trial <- ordered(ResultsDf$Trial,levels=rev(seq(nrow(ResultsDf))))
return(ResultsDf)}
 
trialLengths <- c(3,10,100,1000,10000)
resultsList <- list()
for (i in seq(length(trialLengths))){resultsList[[i]] <- montyMonte(trialLengths[i])}

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.