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

把 Blitzstein 與 Hwang 第二章的 R 小節，改寫成 CS61A 風格的 Python：把述詞當成值傳遞、用次數估計 P(A|B)，並把蒙提霍爾遊戲拆成純邏輯與輸入輸出兩部分。

## 0. 從上一篇接著走

這篇接續[用 Python 計數與模擬](/zh-hant/notes/simulation-in-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$ 很大時，

$$
P(A\mid B)\approx\frac{n_{AB}}{n_B}.
$$

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

### 1.1. 模擬家庭

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

```python
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 的 `==` 比較的是整個值，而一個串列不會等於一個
數字。

### Python 會顯示什麼？比較一個串列

*單選 · zh-hant*

Python 會顯示什麼？

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

- \[ ] elementwise: `[True, False, True]`
- \[ ] false: `False`
- \[ ] error: `TypeError` 錯誤

你的答案: false

#### 提示

1. Python 的 `==` 會拿左邊的整個值和右邊的整個值比較。

2. 一個串列可能等於一個數字嗎？

#### 參考解答

`False`。Python 把整個串列 `[1, 2, 1]` 和數字 `1` 比較，兩者不相等。它不會逐元素比較；在 R 裡，`child1 == 1` 則會得到 `TRUE FALSE TRUE`。要得到 `[True, False, True]`，請寫 `[c == 1 for c in child1]`。

來源: [Adapted from Blitzstein & Hwang, Introduction to Probability, Chapter 2 and §2.10 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

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

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

### Python 會顯示什麼？用 and 計數

*填空 · zh-hant*

Python 會顯示什麼？

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

你的答案: 2 (允許的絕對誤差： 0)

#### 提示

1. 生成式的結果是 `[True, False, False, True]`。

2. `sum` 把每個 `True` 當作 1。

#### 參考解答

`2`。只有第一組和最後一組的兩個值都是 1，所以串列是 `[True, False, False, True]`，加總為 2。這就是 $A$＝兩個都是女孩時的 $n_{AB}$。

來源: [Adapted from Blitzstein & Hwang, Introduction to Probability, Chapter 2 and §2.10 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

### 1.3. 事件就是述詞

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

### 估計條件機率

*Python · zh-hant*

完成 `conditional_frequency(outcomes, a, b)`。`a` 和 `b` 是述詞，也就是接收一個結果、傳回 `True` 或 `False` 的函式。請傳回 $n_{AB}/n_B$：兩者都成立的結果數，除以 `b` 成立的結果數。如果 `b` 從未成立，傳回 `None`。

```python
def conditional_frequency(outcomes, a, b):
    """Estimate P(A | B) as n_AB / n_B over the list OUTCOMES.

    A and B are predicates: functions from an outcome to True or False.
    Return None when B never occurs.

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

#### 公開測試

Elder child is a girl

```pycon
>>> 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
```

At least one girl

```pycon
>>> 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
```

B never happens

```pycon
>>> 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
```

Other outcomes

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

#### 提示

1. 用 `sum(1 for o in outcomes if b(o))` 數出 $n_B$。

2. 用同樣的方式數出 $n_{AB}$，條件改成 `a(o) and b(o)`。

3. 除法之前先檢查 $n_B = 0$ 的情況。

#### 參考解答

```python
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
```

這個函式對小孩或骰子一無所知，只會呼叫傳進來的述詞；正因如此，一個定義就能回答筆記中的每個問題。

來源: [Adapted from Blitzstein & Hwang, Introduction to Probability, Chapter 2 and §2.10 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

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

```python
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)
```

> **範例 · 1.1 · 一次實際執行的結果**
>
> 用上面的種子與 $10^5$ 個家庭，第一個估計值是 $0.5005$，第二個是 $0.3312$。
> 書中例 2.2.5 的精確答案分別是 $1/2$ 與 $1/3$。知道「**老大**是女孩」和知道
> 「**至少一個**是女孩」是不同的資訊，因此答案也不同。

## 2. 蒙提霍爾問題

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

### 2.1. 永遠不換

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

```python
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. 把蒙提寫成模型

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

### 蒙提會開哪一扇門？

*Python · zh-hant*

針對 1、2、3 號門，完成 `monty_opens(chosen, car)`。蒙提從不打開玩家選的門，也不打開有汽車的門。當有兩扇門可以開（玩家正好選中汽車）時，傳回編號較小的那一扇，讓結果可以預測。

```python
def monty_opens(chosen, car):
    """Return the door Monty opens among doors 1, 2, and 3.

    It is neither CHOSEN nor CAR. When two doors are allowed,
    return the lower-numbered one.

    >>> monty_opens(1, 2)
    3
    """
    pass
```

#### 公開測試

Only one door allowed

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

Player picked the car

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

#### 提示

1. 逐一檢查 `(1, 2, 3)` 中的每扇門，只留下既不是 `chosen` 也不是 `car` 的門。

2. 對允許的門取 `min`，就是編號較小的那一扇。

#### 參考解答

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

當 `chosen != car` 時，只有一扇門可以開；兩者相等時有兩扇，`min` 選編號較小的。實際遊戲中蒙提此時會隨機挑選，筆記中的模擬會把這個隨機性加回去。

來源: [Adapted from Blitzstein & Hwang, Introduction to Probability, Chapter 2 and §2.10 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

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

### 換門獲勝的機率

*填空 · zh-hant*

你選 1 號門，蒙提打開一扇後面是山羊的門，而你永遠選擇換門。你贏得汽車的機率是多少？請以分數作答。

你的答案: 2/3 (允許的絕對誤差： 0)

#### 提示

1. 換門正好在你第一次選錯時會贏。

2. 第一次選錯的機率是 $2/3$。

#### 參考解答

$2/3$。如果第一次就選中汽車（機率 $1/3$），換門就會輸。如果第一次選到山羊（機率 $2/3$），蒙提只能打開另一扇山羊門，所以你換過去的那扇門後面一定是汽車。筆記中的模擬得到 $0.6690$。

來源: [Adapted from Blitzstein & Hwang, Introduction to Probability, Chapter 2 and §2.10 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

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

```python
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.1 · 一次實際執行的結果**
>
> 以種子 110 模擬 $10^5$ 局，換門贏的比例是 $0.6690$，接近 $2/3$。永遠不換只在
> 第一次就選對時才贏，所以換門正好在第一次選錯時贏。

### 2.3. 實際玩一局

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

```python
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)$。

> **定理 · 3.1 · 用補集寫出的貝氏定理**
>
> 若 $0<P(A)<1$ 且 $P(E)>0$，則
>
> $$
> P(A\mid E)=\frac{P(E\mid A)\,P(A)}{P(E\mid A)\,P(A)+P(E\mid A^c)\,P(A^c)}.
> $$
>
> 分母就是 $P(E)$，由全機率公式展開。

### 依證據更新機率

*Python · zh-hant*

完成 `bayes_update(prior, p_evidence_if_a, p_evidence_if_not_a)`。已知 $P(A)$、$P(E\mid A)$ 與 $P(E\mid A^c)$，請用貝氏定理搭配全機率公式傳回 $P(A\mid E)$。測試傳入的是 `Fraction`，所以答案是精確值。假設 $P(E)>0$。

```python
from fractions import Fraction


def bayes_update(prior, p_evidence_if_a, p_evidence_if_not_a):
    """Return P(A | E) from the prior P(A), P(E | A), and P(E | not A).

    >>> bayes_update(Fraction(1, 4), 1, Fraction(2, 3))
    Fraction(1, 3)
    """
    pass
```

#### 公開測試

Monday: at least one girl

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

Tuesday: the elder is a girl

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

Evidence that rules out not-A

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

Evidence that says nothing

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

#### 提示

1. 分子是 $P(E\mid A)\,P(A)$。

2. 分母是在分子之外再加上 $P(E\mid A^c)\,(1-P(A))$。

#### 參考解答

```python
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$ 的機率，所以可以直接當成下一次的先驗再傳進去，這就是依序更新。在 $A$ 與非 $A$ 之下同樣可能出現的證據，不會改變先驗。

來源: [Adapted from Blitzstein & Hwang, Introduction to Probability, Chapter 2 and §2.10 (R), rewritten in Python](http://stat110.net) · Original exercise; the book's R examples are paraphrased, not copied

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

> **範例 · 3.2 · 星期一，然後星期二**
>
> \*\*星期一：\*\*得知 $B$＝至少一個是女孩。若 $A$ 成立，$B$ 一定成立；若不成立，
> 在其餘三種同樣可能的家庭中，$B$ 在兩種裡成立，所以 $P(B\mid A^c)=2/3$。
>
> ```python
> >>> monday = bayes_update(Fraction(1, 4), 1, Fraction(2, 3))
> >>> monday
> Fraction(1, 3)
> ```
>
> \*\*星期二：\*\*又得知 $C$＝老大是女孩。在 $B$ 成立而 $A$ 不成立時，家庭是
> 「女孩、男孩」或「男孩、女孩」，所以 $P(C\mid A^c,B)=1/2$。
>
> ```python
> >>> bayes_update(monday, 1, Fraction(1, 2))
> Fraction(1, 2)
> ```
>
> 兩個答案都和第 1 節的模擬一致。先得知 $B$ 再得知 $C$，等於直接得知「老大是
> 女孩」，這也是檢查計算的好方法。

`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](http://stat110.net)。
- [CS61A](https://cs61a.org/) 與 John DeNero 的
  [《Composing Programs》](https://www.composingprograms.com/)：高階函式、
  doctest，以及把邏輯與輸入輸出分開。
- Python 官方文件：[`random`](https://docs.python.org/3/library/random.html)、
  [`fractions`](https://docs.python.org/3/library/fractions.html)、
  [`input`](https://docs.python.org/3/library/functions.html#input)。