用 Python 看條件機率：兩個小孩、蒙提霍爾問題，以及依證據更新

0. 從上一篇接著走

這篇接續用 Python 計數與模擬， 上一篇改寫了第一章的 R 小節。Blitzstein 與 Hwang《Introduction to Probability》的第二章談條件機率，它的 R 小節（§2.10）用模擬驗證兩個著名的 結果：兩個小孩問題與蒙提霍爾問題（Monty Hall problem）。寫法和上一篇一樣： 有名字的小函式、doctest，以及把函式當成值傳遞。只用標準函式庫就夠了，這次 也不需要上一篇的 replicate。

1. 把條件機率看成次數

條件機率的頻率學派解讀是：把實驗重複 $n$ 次，數出 $B$ 發生的次數 $n_B$， 以及 $A$ 與 $B$ 同時發生的次數 $n_{AB}$。當 $n$ 很大時，

我們只看 $B$ 發生的那些回合，再問其中 $A$ 發生了多少次。

1.1. 模擬家庭

假設每個家庭有兩個小孩，每個小孩各自獨立、以相同機率是女孩（1）或男孩（2）。 R 會建立兩個向量，每個小孩一個；Python 也可以這樣做，再用 zip 把它們配成 一個家庭一個 tuple：

import random

random.seed(110)
n = 10 ** 5
child1 = random.choices([1, 2], k=n)     # 每個家庭的老大
child2 = random.choices([1, 2], k=n)     # 老二
families = list(zip(child1, child2))

families[0] 可能是 (1, 2)：老大是女孩、老二是男孩。把同一個家庭的兩個 小孩放在同一個 tuple 裡，就不可能不小心把一個家庭的老大和另一個家庭的老二 配在一起。

1.2. 沒有逐元素比較

R 的 child1 == 1 會把向量的每個元素都拿來和 1 比較，傳回一個由 TRUE、 FALSE 組成的向量。Python 的 == 比較的是整個值，而一個串列不會等於一個 數字。

>>> child1 = [1, 2, 1]
>>> child1 == 1

所以 Python 用串列生成式把迴圈寫出來。R 的 & 和 | 變成 Python 的 and 和 or，寫在生成式裡面：

>>> pairs = [(1, 1), (1, 2), (1, 1)]
>>> [c1 == 1 and c2 == 1 for c1, c2 in pairs]
[True, False, True]

>>> pairs = [(1, 1), (1, 2), (2, 1), (1, 1)]
>>> sum([c1 == 1 and c2 == 1 for c1, c2 in pairs])

1.3. 事件就是述詞

與其每個問題都重寫一次計數，不如把每個事件描述成一個述詞（predicate）： 一個接收結果、傳回 True 或 False 的函式。這樣一個函式就能估計任何條件 機率。這正是 CS61A 所說的「把函式當成參數傳遞」。

>>> families = [(1, 1), (1, 2), (2, 1), (2, 2)]
>>> both_girls = lambda f: f == (1, 1)
>>> conditional_frequency(families, both_girls, lambda f: f[0] == 1)
0.5

>>> families = [(1, 1), (1, 2), (2, 1), (2, 2)]
>>> both_girls = lambda f: f == (1, 1)
>>> round(conditional_frequency(families, both_girls, lambda f: 1 in f), 3)
0.333

>>> families = [(1, 1), (1, 2), (2, 1), (2, 2)]
>>> both_girls = lambda f: f == (1, 1)
>>> conditional_frequency(families, both_girls, lambda f: False) is None
True

>>> rolls = [1, 2, 3, 4, 5, 6]
>>> conditional_frequency(rolls, lambda r: r == 6, lambda r: r % 2 == 0)
0.3333333333333333

def conditional_frequency(outcomes, a, b):
    n_b = sum(1 for o in outcomes if b(o))
    n_ab = sum(1 for o in outcomes if a(o) and b(o))
    return None if n_b == 0 else n_ab / n_b

寫好 conditional_frequency 之後，書中的兩個問題各只要一行：

both_girls = lambda f: f == (1, 1)
elder_girl = lambda f: f[0] == 1
at_least_one_girl = lambda f: 1 in f

conditional_frequency(families, both_girls, elder_girl)
conditional_frequency(families, both_girls, at_least_one_girl)

2. 蒙提霍爾問題

三扇門後面，一扇是汽車，另外兩扇是山羊。你先選一扇門。知道汽車在哪裡的主持人 蒙提，會打開另外兩扇中的一扇，讓你看到山羊，然後問你要不要換。該換嗎？

2.1. 永遠不換

假設你總是選 1 號門。那麼「永遠不換」這個策略只有在汽車就在 1 號門後面時 才會贏，連蒙提都不必模擬：

random.seed(110)
cardoor = random.choices([1, 2, 3], k=n)
sum(car == 1 for car in cardoor) / n

以種子 110 執行一次得到 $0.3313$，接近 $1/3$。

2.2. 把蒙提寫成模型

要模擬換門，就需要蒙提的選擇。如果你選的門後面是山羊，蒙提只有一扇門可以開； 如果你選的正是汽車，他有兩扇可以選，實際遊戲中他會隨機挑一扇。為了讓程式容易 測試，下一題改成固定的選擇：永遠開允許的門中編號最小的那一扇。隨機性之後可以 再從外面加回來。

>>> monty_opens(1, 2)
3
>>> monty_opens(2, 1)
3
>>> monty_opens(3, 1)
2

>>> monty_opens(1, 1)
2
>>> monty_opens(2, 2)
1
>>> monty_opens(3, 3)
1

def monty_opens(chosen, car):
    return min(door for door in (1, 2, 3) if door != chosen and door != car)

蒙提開門之後，換門就是改選那扇既不是你的、也不是他開的門。1、2、3 號門加起來 是 6，所以那扇門是 6 - chosen - opened。

要驗證這個答案，就模擬整個遊戲，並讓蒙提在有選擇時隨機挑門：

def play_switching():
    """Pick door 1, watch Monty open a goat door, switch; return True on a win."""
    car = random.randint(1, 3)
    chosen = 1
    opened = random.choice([d for d in (1, 2, 3) if d != chosen and d != car])
    return 6 - chosen - opened == car

random.seed(110)
sum(play_switching() for _ in range(n)) / n

2.3. 實際玩一局

書中最後寫了一個互動式的 monty() 函式。在 Python 裡，input 扮演 R 中 scan 的角色：它會印出提示，並傳回玩家輸入的內容，而且一定是字串，所以要 用 int(...) 把門號轉成整數。f-string 取代 R 的 paste，startswith 則取代 substr(reply, 1, 1)。

import random

def switch_door(chosen, opened):
    """The one door that is neither CHOSEN nor OPENED."""
    return 6 - chosen - opened

def monty():
    """Play one round of the Monty Hall game at the terminal."""
    doors = [1, 2, 3]
    car = random.choice(doors)                       # 汽車所在的門
    print("Monty Hall says 'Pick a door, any door!'")
    chosen = int(input('Your door (1, 2, or 3): '))
    # 蒙提打開一扇既不是你的、也不是汽車的門。
    opened = random.choice([d for d in doors if d != chosen and d != car])
    print(f'Monty opens door {opened}!')
    reply = input('Would you like to switch (y/n)? ')
    if reply.strip().lower().startswith('y'):        # "yes"、"yep"、"Yeah!"
        chosen = switch_door(chosen, opened)
    print('You won!' if chosen == car else 'You lost!')

把它存成 monty.py，執行 python3 -i monty.py，之後想玩幾次就呼叫幾次 monty()。本頁的練習無法執行它，因為練習環境沒有鍵盤輸入。

注意這個程式的結構。所有不需要有人坐在鍵盤前就能測試的決定（switch_door， 以及練習中的 monty_opens）都是純函式；只有最外層的 monty 會和玩家互動。 CS61A 把這稱為「把計算與輸入輸出分開」，這也是前面的練習能用 doctest 檢查 邏輯的原因。

3. 後驗變成先驗

第二章結尾（圖 2.7）描述了證據如何一步一步到來。還沒有任何證據時，$P(A)$ 是 先驗機率。得知 $B$ 之後，用貝氏定理算出後驗機率 $P(A\mid B)$。這個後驗又成為 下一步的先驗，再得知 $C$ 就得到 $P(A\mid B,C)$。

>>> bayes_update(Fraction(1, 4), 1, Fraction(2, 3))
Fraction(1, 3)

>>> bayes_update(Fraction(1, 3), 1, Fraction(1, 2))
Fraction(1, 2)

>>> bayes_update(Fraction(1, 100), 1, 0)
Fraction(1, 1)

>>> bayes_update(Fraction(1, 5), Fraction(1, 2), Fraction(1, 2))
Fraction(1, 5)

def bayes_update(prior, p_evidence_if_a, p_evidence_if_not_a):
    numerator = p_evidence_if_a * prior
    return numerator / (numerator + p_evidence_if_not_a * (1 - prior))

一次更新一項證據，就是一連串的函式呼叫，每次的輸出成為下一次的輸入。再看一次 兩個小孩的家庭，令 $A$＝兩個都是女孩，所以 $P(A)=1/4$。

>>> monday = bayes_update(Fraction(1, 4), 1, Fraction(2, 3))
>>> monday
Fraction(1, 3)

>>> bayes_update(monday, 1, Fraction(1, 2))
Fraction(1, 2)

fractions.Fraction 讓計算保持精確，$1/3$ 就一直是 $1/3$，不會變成帶有 捨入誤差的 $0.333\ldots$。

參考資料

• Joseph K. Blitzstein、Jessica Hwang，《Introduction to Probability》， 第 2 章：例 2.2.5（兩個小孩）、蒙提霍爾問題、圖 2.7（依序更新）與 §2.10（R）。課程資源：stat110.net。
• CS61A 與 John DeNero 的 《Composing Programs》：高階函式、 doctest，以及把邏輯與輸入輸出分開。
• Python 官方文件：random、 fractions、 input。