# 用 Python 計數與模擬：配對問題與生日問題

把 Blitzstein 與 Hwang 第一章的 R 小節，改寫成 CS61A 風格的 Python：小而有名字的函式、doctest，以及當成值傳遞的模擬。

## 0. 這篇筆記要做什麼

Blitzstein 與 Hwang 的《Introduction to Probability》第一章結尾有一小節 R
程式：向量、階乘、隨機抽樣，以及兩個用來驗證前文結果的模擬。這篇筆記用
Python 走一遍同樣的內容。

這不是逐行翻譯。R 以向量為核心；Python 以名稱、函式和串列（list）為核心。
所以下面的程式採用 CS61A 的習慣：

- 每個想法都寫成一個**有名字的小函式**，並附上 docstring；
- 每個 docstring 都帶有 **doctest**，也就是 Python 可以自動檢查的範例；
- 模擬是一個**被當成值傳遞的函式**，而不是一個特殊指令。

全部只用標準函式庫（`math`、`random`、`collections`），因此在任何 Python 3
環境，以及本頁的練習裡都能執行。

## 1. 從向量到串列

R 的 `c(3,1,4,1,5,9)` 在 Python 裡是一個串列：

```python
>>> v = [3, 1, 4, 1, 5, 9]
>>> sum(v), max(v), min(v), len(v)
(23, 9, 1, 6)
```

數列 $1,2,\dots,n$ 在 R 寫成 `1:n`，在 Python 寫成 `range(1, n + 1)`。
range 會在第二個參數**之前**停下，所以要加 1。從 $m$ 倒數到 $n$ 則需要負的
步長：`range(m, n - 1, -1)`。

### 1.1. 位置從零開始

最大的差別是索引。R 的位置從 1 開始編號，Python 從 0 開始。所以 `v[0]` 是第
一個元素，書中所說的第 $i$ 個元素是 `v[i - 1]`。

### Python 會顯示什麼？索引

*填空 · zh-hant*

Python 會顯示什麼？

```python
>>> v = [3, 1, 4, 1, 5, 9]
>>> v[1]
```

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

#### 提示

1. Python 的位置從 0 開始編號，不是從 1。

2. `v[0]` 是 `3`，下一個是什麼？

#### 參考解答

`1`。Python 的位置從 0 開始，所以 `v[1]` 是**第二個**元素。在 R 裡，`v[1]` 則是第一個元素 `3`。

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

取出第 1、3、5 個元素（R 的 `v[c(1,3,5)]`）可以用串列生成式走過這些位置。
去掉第 2 到第 4 個元素（R 的 `v[-(2:4)]`）則是把兩段切片接起來：

```python
>>> [v[i - 1] for i in (1, 3, 5)]
[3, 4, 5]
>>> v[:1] + v[4:]
[3, 5, 9]
```

負的索引在 Python 裡意思不同：`v[-1]` 是**最後**一個元素，而不是「除了第一
個以外的全部」。

### 1.2. 逐元素運算要寫出來

R 會自動把算術套用到向量的每個元素上；Python 則用串列生成式明確寫出來。向量
$\left(1,\frac{1}{2^2},\frac{1}{3^2},\dots,\frac{1}{100^2}\right)$ 寫成：

```python
>>> squares = [1 / k ** 2 for k in range(1, 101)]
>>> squares[:3]
[1.0, 0.25, 0.1111111111111111]
```

次方用 `**`。在 Python 裡 `^` 是位元互斥或（XOR），所以 `2 ^ 3` 是 `1`，
不是 `8`。這是把 R 程式搬到 Python 時最容易犯的錯之一。

R 還會把較短的向量**循環使用**（recycling）：`v + 3` 會把每個元素都加 3。
Python 的串列不會這樣做。

### Python 會顯示什麼？沒有循環使用

*單選 · zh-hant*

這裡會發生什麼事？

```python
>>> v = [3, 1, 4]
>>> v + 3
```

- \[ ] added: `[6, 4, 7]`
- \[ ] appended: `[3, 1, 4, 3]`
- \[ ] error: `TypeError` 錯誤

你的答案: error

#### 提示

1. 兩個串列之間的 `+` 會把它們接起來。`3` 是串列嗎？

2. Python 不會把較短的值重複到和較長的一樣長。

#### 參考解答

Python 會引發 `TypeError: can only concatenate list (not "int") to list`。對串列而言，`+` 是串接，兩邊都必須是串列。要讓每個元素加 3，寫成 `[x + 3 for x in v]`，得到 `[6, 4, 7]`。

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

## 2. 階乘與二項式係數

標準函式庫兩者都有：`math.factorial(n)` 與 `math.comb(n, k)`。

```python
>>> import math
>>> math.factorial(20)
2432902008176640000
>>> math.comb(10, 3)
120
```

Python 的整數沒有固定大小，所以 `math.factorial(1000)` 是一個精確的 2568 位
整數。上限只在轉成浮點數時出現：`float(math.factorial(171))` 會引發
`OverflowError`，正好是 R 開始回傳 `Inf` 的地方。$n$ 很大時，改用對數計算。
因為 $n!=\Gamma(n+1)$，自然對數 $\log(n!)$ 就是 `math.lgamma(n + 1)`；
$\log\binom{n}{k}$ 則是 `lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1)`。

自己寫 $\binom{n}{k}$ 是很好的第一個遞迴練習。基本情況直接來自定義，遞迴情況
則是巴斯卡法則。

> **定理 · 2.1 · 巴斯卡法則**
>
> 對整數 $n\ge 1$ 與 $1\le k\le n-1$，
>
> $$
> \binom{n}{k}=\binom{n-1}{k-1}+\binom{n-1}{k}.
> $$
>
> 在 $n$ 個人中指定某一個人。大小為 $k$ 的委員會要嘛包含這個人，剩下 $k-1$
> 人從其餘 $n-1$ 人中選；要嘛不包含，$k$ 人全都從其餘 $n-1$ 人中選。

### 用遞迴計算二項式係數

*Python · zh-hant*

用**遞迴**完成 `choose(n, k)`，依據巴斯卡法則 $\binom{n}{k}=\binom{n-1}{k-1}+\binom{n-1}{k}$。不要使用 `math.comb` 或階乘。$k<0$ 或 $k>n$ 時傳回 0。

```python
def choose(n, k):
    """Return the number of ways to choose K items from N, using Pascal's rule.

    >>> choose(5, 2)
    10
    """
    pass
```

#### 公開測試

Small case

```pycon
>>> choose(5, 2)
10
```

Middle of a row

```pycon
>>> choose(10, 5)
252
```

Choose none or all

```pycon
>>> choose(4, 0)
1
>>> choose(4, 4)
1
```

Too many

```pycon
>>> choose(3, 5)
0
```

#### 提示

1. 有兩種基本情況：不可能的選法，以及一個都不選或全選。

2. $k$ 為 0 或等於 $n$ 時，恰好只有一種選法。

3. 其他情況，傳回兩個較小 `choose` 呼叫的和。

#### 參考解答

```python
def choose(n, k):
    if k < 0 or k > n:
        return 0
    if k == 0 or k == n:
        return 1
    return choose(n - 1, k - 1) + choose(n - 1, k)
```

第一個基本情況處理不可能的要求，第二個是巴斯卡三角形的邊緣。每次遞迴呼叫都會讓 $n$ 變小，所以遞迴一定會結束。它會重複計算很多次，$n$ 小時沒有問題，這也正是之後學習記憶化（memoization）的好理由。

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

## 3. 隨機抽樣

`random` 模組扮演 R 中 `sample` 的角色。書中的位置是 $1,\dots,n$，所以從
`range(1, n + 1)` 抽樣。

| R | Python | 意思 |
| --- | --- | --- |
| `sample(n, k)` | `random.sample(range(1, n + 1), k)` | 不放回抽 $k$ 個 |
| `sample(n, k, replace=TRUE)` | `random.choices(range(1, n + 1), k=k)` | 放回抽 $k$ 個 |
| `sample(n)` | `random.sample(range(1, n + 1), n)` | 隨機排列 |
| `sample(letters, 7)` | `random.sample(string.ascii_lowercase, 7)` | 七個不重複的字母 |
| `sample(4, 3, replace=TRUE, prob=…)` | `random.choices(range(1, 5), weights=[0.1, 0.2, 0.3, 0.4], k=3)` | 加權、放回 |

`random.choices` 一定是放回抽樣；標準函式庫沒有加權的不放回抽樣。

這些數字是偽隨機的，由一個確定性的演算法產生。這其實很有用：先呼叫
`random.seed(110)`，模擬就能重現，別人執行你的程式也會得到同樣的數字。

## 4. 模擬就是一個函式

R 的 `replicate(times, expression)` 會把一個運算式重複求值許多次。在 Python
裡，我們直接把實驗本身，也就是一個不帶參數的函式傳進去，要呼叫幾次都可以：

```python
def replicate(times, experiment):
    """Call EXPERIMENT, a function of no arguments, TIMES times.
    Return the results as a list.

    >>> replicate(3, lambda: 'hi')
    ['hi', 'hi', 'hi']
    """
    return [experiment() for _ in range(times)]
```

這就是 CS61A 所說的**高階函式**。`replicate` 對紙牌或生日一無所知，它只知道
怎麼重複。`lambda` 把一個運算式包成函式，讓它延後、而且每次呼叫都重新求值，
而不是只算一次。

## 5. 配對問題

把標上 $1,\dots,n$ 的一副牌洗亂，再一張一張翻開。如果第 $i$ 張翻開的正好是
$i$ 號牌，就叫做一次**配對**。書中證明（例 1.6.4）至少出現一次配對的機率是

$$
1-\frac{1}{2!}+\frac{1}{3!}-\cdots+(-1)^{n+1}\frac{1}{n!}\;\approx\;1-\frac{1}{e}\approx 0.632.
$$

要模擬，先得數出一副牌裡有幾次配對。在 Python 的算術裡，`True` 和 `False`
的行為就像 `1` 和 `0`，和 R 的 `TRUE`、`FALSE` 一樣。

### Python 會顯示什麼？把真假值相加

*填空 · zh-hant*

Python 會顯示什麼？

```python
>>> deck = [2, 1, 3, 4]
>>> sum([card == position for position, card in enumerate(deck, start=1)])
```

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

#### 提示

1. 裡面的串列是 `[False, False, True, True]`。

2. 在算術中，`True` 當作 `1`，`False` 當作 `0`。

#### 參考解答

`2`。3 號和 4 號牌在自己的位置上，所以串列是 `[False, False, True, True]`，加總就是在數 `True` 的個數。

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

### 數出配對數

*Python · zh-hant*

完成 `count_fixed_points(deck)`。`deck` 依翻開順序列出牌；第 $i$ 張翻開的若是 $i$ 號牌，就算一次配對。請傳回配對的次數。

```python
def count_fixed_points(deck):
    """Return how many cards in DECK sit in their own position.

    Positions are numbered from 1, like the cards.

    >>> count_fixed_points([1, 2, 3])
    3
    >>> count_fixed_points([2, 1, 3])
    1
    """
    pass
```

#### 公開測試

Every card matches

```pycon
>>> count_fixed_points([1, 2, 3])
3
```

One match

```pycon
>>> count_fixed_points([2, 1, 3])
1
```

No matches

```pycon
>>> count_fixed_points([3, 1, 2])
0
>>> count_fixed_points([2, 1])
0
```

Empty deck

```pycon
>>> count_fixed_points([])
0
```

#### 提示

1. `enumerate(deck, start=1)` 會同時給出每個位置和那張牌。

2. 把每張牌和它的位置比較，再把結果加總。

#### 參考解答

```python
def count_fixed_points(deck):
    return sum(card == position for position, card in enumerate(deck, start=1))
```

每次比較的結果是 `True` 或 `False`，`sum` 會把它們當作 1 或 0。空的牌堆沒有配對，加總為 0。

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

寫好 `count_fixed_points` 之後，洗一次牌和完整的模擬都只要一兩行：

```python
import random

def shuffled_matches(n):
    """Shuffle cards 1..N and count the matches."""
    deck = random.sample(range(1, n + 1), n)
    return count_fixed_points(deck)

random.seed(110)
results = replicate(10 ** 4, lambda: shuffled_matches(100))
sum(r >= 1 for r in results) / 10 ** 4     # 至少有一次配對的比例
```

> **範例 · 5.1 · 一次實際執行的結果**
>
> 以種子 110 對 100 張牌洗 $10^4$ 次，至少有一次配對的比例是 $0.6267$，
> 理論值是 $1-1/e\approx0.6321$。換一個種子或 Python 版本，比例可能稍有不同；
> 兩者的差距大約以 $1/\sqrt{10^4}=0.01$ 的量級縮小。

## 6. 生日問題

$k$ 個人、365 個同樣可能的生日（忽略 2 月 29 日）時，至少兩人同一天生日的
機率是

$$
1-\frac{365\cdot364\cdots(365-k+1)}{365^k}.
$$

`math.prod` 可以直接算出乘積，精確的整數運算讓結果一直準確到最後的除法：

```python
>>> k = 23
>>> 1 - math.prod(range(365 - k + 1, 366)) / 365 ** k
0.5072972343239854
```

寫成帶有 doctest 的函式，一次乘一個因子：

```python
def birthday_match(k, days=365):
    """Probability that at least two of K people share a birthday.

    >>> round(birthday_match(23), 3)
    0.507
    >>> birthday_match(1)
    0.0
    >>> birthday_match(366)
    1.0
    """
    no_match = 1.0
    for i in range(k):
        no_match *= (days - i) / days
    return 1 - no_match
```

R 的 `qbirthday(p)` 回答反過來的問題：要多少人，機率才會達到 $p$？Python
的標準函式庫沒有對應的函式，那就自己寫。

### 需要多少人？

*Python · zh-hant*

題目已提供 `birthday_match(k, days)`。請完成 `people_needed(p, days)`：傳回讓 `birthday_match(k, days)` 至少為 `p` 的最小人數 $k\ge1$。這正是 R 的 `qbirthday` 所計算的。假設 $0<p\le1$。

```python
def birthday_match(k, days=365):
    """Probability that at least two of K people share a birthday."""
    no_match = 1.0
    for i in range(k):
        no_match *= (days - i) / days
    return 1 - no_match


def people_needed(p, days=365):
    """Return the smallest number of people for which
    birthday_match reaches at least P.

    >>> people_needed(0.5)
    23
    """
    pass
```

#### 公開測試

Even odds

```pycon
>>> people_needed(0.5)
23
```

Ninety percent

```pycon
>>> people_needed(0.9)
41
```

Ninety-nine percent

```pycon
>>> people_needed(0.99)
57
```

A shorter year

```pycon
>>> people_needed(0.5, days=2)
2
```

#### 提示

1. 從一個人開始，逐次加人，直到機率夠大。

2. 適合用 `while` 迴圈：當 `birthday_match(k, days) < p` 時就繼續。

#### 參考解答

```python
def people_needed(p, days=365):
    k = 1
    while birthday_match(k, days) < p:
        k += 1
    return k
```

迴圈一定會停：有 `days + 1` 個人時必定有人同天生日，機率會達到 1。

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

### 6.1. 改用模擬

模擬時，放回抽出 23 個生日，再看最常見的生日出現了幾次。
`collections.Counter` 做的就是 R 中 `tabulate` 的工作：

```python
from collections import Counter

def most_shared(k, days=365):
    """Draw K birthdays; return the size of the largest group sharing one."""
    birthdays = random.choices(range(1, days + 1), k=k)
    return max(Counter(birthdays).values())

random.seed(110)
r = replicate(10 ** 4, lambda: most_shared(23))
sum(x >= 2 for x in r) / 10 ** 4           # 至少一組同天生日
sum(x >= 3 for x in r) / 10 ** 4           # 至少一組三人同天
```

模擬能回答公式回答不了的問題。R 的 `pbirthday(23, coincident=3)` 給出 23
人中出現**三人**同天生日的機率約為 $0.014$，`qbirthday(0.5, coincident=3)`
則給出 88。Python 沒有處理三人同天的內建函式，但 `most_shared` 不必多寫任何
程式就能處理。

> **範例 · 6.1 · 一次實際執行的結果**
>
> 以種子 110 模擬 $10^4$ 組，23 人中有人同天生日的比例是 $0.5039$（精確值
> $0.5073$），出現三人同天的比例是 $0.013$。88 人時，三人同天的比例是
> $0.5073$，與 R 給出的 88 人一致。

生日機率不相等時，公式會變得很難，模擬卻依然簡單：把 `weights=` 傳給
`random.choices` 就好。

## 7. 註解與 docstring

兩種語言都用 `#` 開始一行註解。Python 另外有 docstring，也就是函式最上方說明
它做什麼的字串。裡面放了 doctest，說明文件同時也是測試。執行

```bash
python3 -m doctest simulation.py -v
```

Python 就會檢查檔案裡每一個 `>>>` 範例。docstring 寫函式承諾了什麼，
註解寫某一行為什麼這樣寫。

## 參考資料

- Joseph K. Blitzstein、Jessica Hwang，《Introduction to Probability》，
  第 1 章 §1.6（例 1.6.4，配對問題）與 §1.8（R）。
  課程資源：[stat110.net](http://stat110.net)。
- [CS61A](https://cs61a.org/) 與 John DeNero 的
  [《Composing Programs》](https://www.composingprograms.com/)，本文的寫法
  來自這裡：高階函式、doctest，以及「Python 會顯示什麼？」
- Python 官方文件：[`random`](https://docs.python.org/3/library/random.html)、
  [`math`](https://docs.python.org/3/library/math.html)、
  [`collections.Counter`](https://docs.python.org/3/library/collections.html#collections.Counter)。