Solving Monty Hall Problem with Monte Carlo Simulation in R
Monty Hall Problem 是許多第一次接觸機率論或初統的同學會很容易卡住的問題。同學們常常會思索著 information inflow 跟 update 條件機率必不太直觀，因此也不容易接受 Monty Hall Problem 的結論。以下我們嘗試在 R 當中寫一個 function 來模擬 Monty Hall Problem，然後再透過 Monte Carlo 模擬分別在兩種情況下反覆玩這個遊戲後的結果。
解析解
我們可以令 A 表示贏得車的事件；令 R 表示第一次選擇猜中車的事件
在選擇「換」的情況下：
𝑃(𝐴)=𝑃(𝐴∩𝑅)+𝑃(𝐴∩𝑅′)=𝑃(𝐴|𝑅)𝑃(𝑅)+𝑃(𝐴|𝑅′)𝑃(𝑅′)=0×13+1×23=23
至於「不換」的情況我們也可以如法炮製，答案就會是 1/3
定義函數
在這裡，我們先以最經典的三門版本來寫。
monty_hall = function(switch = logical()){
#randomly arrange the door
doors = c(1,2,3)
names(doors) = sample(c("goat", "car", "goat"))
prize_door = which(names(doors)=='car') #find the index of car
#Now Guess
guess = sample(doors, 1)
if(guess == prize_door){
revealed_door = sample(doors[doors != prize_door], 1) #randomly open one of the goat door
}else{
revealed_door = doors[doors != prize_door & doors != guess] #open the other goat door
}
#Show the return
if(switch){
switched_door = doors[doors != guess & doors != revealed_door]
return(prize_door == switched_door)
} else {
return(prize_door == guess)
}
}首先，這個函數的輸入只有一個，就是是否要「換門」；而輸出則是「是否得到大獎」。如此一來我們只要反覆執行這個函數，並且統計「平均而言我們得到大獎的相對頻率」，就可以估算換與不換情況下贏得大獎的機率。
解釋
那這個函數的內容是什麼意思呢？又是如何運作的呢？我們大致解釋如下：
doors代表遊戲初始的三道門
透過names()函式來將羊（槓龜）以及車（大獎）分別隨機地賦予給三道門prize_door則代表那道背後是車的門guess代表最起初我們選擇哪一道門只能靠猜，所以隨機地選。這裡要改成都固定選某一道門也是可以的，並不會影響結果。revealed_door表示被主持人打開的門
接下來就會進入比較麻煩的條件判斷，因為會打開哪一扇門取決於我們初始的猜測猜到哪一道門。假如我們選中的門背後就是車，那主持人會在剩下兩道被後是羊的門隨機選擇一道開啟。若起初的猜測並不是車的話，代表我們選到羊，那主持人就只能開啟唯一一道剩下羊的門，表現在程式裡就是對我們選到的以及車對應的那些門取否。
然後就輕鬆了，只要根據輸入來決定打開哪一扇門就行了。如果不想換，那就直接打開初始猜測，然後比較看看跟車對應的門是不是同一扇就行了。
如果要換的話，那就是在剩下的唯二道門（主持人一定會打開一定數量的門直到剩下兩道，其中一道會是初始猜測的門），扣除掉初始猜測的那一道門。注意到，這代表 Monty Hall Problem 當然有多門版本的拓展，我們會再最下面提供程式範例。
模擬結果
接著我們就可以來實驗看看了。透過 replicate() 函式來重複執行我們自定義的函式。
mean(replicate(monty_hall(T), n = 10000))
mean(replicate(monty_hall(F), n = 10000))
result = replicate(monty_hall(T), n = 10000)
barplot(c(sum(result), length(result)-sum(result)),
main = 'The frequency of winning the car in
Monty Hall Problem in 10000 rounds',names.arg = c('Change', 'Do not change'), ylim = c(0,10000))可以看到實驗的結果如下
可以發現，如果我們總是選擇換門的話，在一萬次遊戲中，平均而言可以贏得大獎6700次，大約是 2/3，意味著選擇換門能贏得大獎的機率很接近 2/3。而不換的話則是有接近 1/3 的機率贏得大獎。
多門問題
如同上文提到的，我們可以拓展原先的問題，變成多道門的版本。例如原先有 100 道門，其中只有一道門背後是車，其餘 99 道門後都是山羊。主持人一樣會在你做出選擇後，打開 98 道門，留下你選的初始選擇跟有大獎的那道門，讓你二選一。遊戲的邏輯類似，所以我們只需要簡單修改我們的函式就可以再次運用蒙地卡羅法來模擬條件機率。
# n = 4
monty_hall_n = function(switch = logical(), n){
#randomly arrange the door
doors = 1:n
names(doors) = sample(rep(c("goat", "car"), c(n-1,1)))
prize_door = which(names(doors)=='car') #find the index of car
#Now Guess
guess = sample(doors, 1)
if(guess == prize_door){
revealed_door = sample(doors[doors != prize_door], n-2) #randomly open one of the goat door
}else{
revealed_door = doors[!doors %in% c(prize_door, guess)] #open the other goat door
}
#Show the return
if(switch){
switched_door = doors[!doors %in% c(guess, revealed_door)]#'%in%' is value matching
return(prize_door == switched_door)
} else {
return(prize_door == guess)
}
}
monty_hall_n(T, 3)
mean(replicate(monty_hall_n(T, 4), n = 10000))
mean(replicate(monty_hall_n(F, 4), n = 10000))
mean(replicate(monty_hall_n(T, 100), n = 10000))
mean(replicate(monty_hall_n(F, 100), n = 10000))
result = replicate(monty_hall_n(T, 100), n = 10000)
barplot(c(sum(result), length(result)-sum(result)),
main = 'The frequency of winning the car in
Monty Hall Problem with 100 doors in 10000 rounds',
names.arg = c('Succeed', 'Fail'), ylim = c(0,10000))這個 n 可以是任意正整數。當然，當 n = 3 時會跟上面的三門問題得到一樣的結果。至於這個條件機率該怎麼計算，跟三門問題並無二致，我們就留待讀者自行思考了。
延伸討論 — 如果我們不總是決定「換」與「不換」？
上面提到的都是如果我們始終堅定地採取換或者不換的策略之下，我們可以贏得大獎的相對頻率。我們可以思考一下，如果我們回到三門問題，然後每次遇到換與不換的抉擇時，都採用拋硬幣的方式來決定，那我們贏得大獎的機率會是多少呢？
我們一樣可以模擬一下，透過隨機選擇函式的輸入，我們模擬每次抉擇時都用拋硬幣的方式來決定換與不換。
monty_hall(sample(c(T, F), 1))
mean(replicate(expr = monty_hall(sample(c(T, F), 1)), 10000))執行後可以發現，贏得大獎的相對頻率很接近 1/2 。直覺來說，剩下的兩道門中我隨機地選擇一道，當然獲獎的機率就是 1/2 囉！這其實也就代表著，我們完全沒有用到 information inflow ，而能否得到大獎的機率則會按照我們分配給「換」與「不換」的比重來乘以對應的機率得到。
