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

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 裡是一個串列：

>>> 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]。

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

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

>>> [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)$ 寫成：

>>> 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 的串列不會這樣做。

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

2. 階乘與二項式係數

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

>>> 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}$ 是很好的第一個遞迴練習。基本情況直接來自定義，遞迴情況 則是巴斯卡法則。

>>> choose(5, 2)
10

>>> choose(10, 5)
252

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

>>> choose(3, 5)
0

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)

3. 隨機抽樣

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

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

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

4. 模擬就是一個函式

R 的 replicate(times, expression) 會把一個運算式重複求值許多次。在 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）至少出現一次配對的機率是

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

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

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

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

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

>>> count_fixed_points([])
0

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

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

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     # 至少有一次配對的比例

6. 生日問題

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

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

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

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

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 的標準函式庫沒有對應的函式，那就自己寫。

>>> people_needed(0.5)
23

>>> people_needed(0.9)
41

>>> people_needed(0.99)
57

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

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

6.1. 改用模擬

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

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 不必多寫任何 程式就能處理。

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

7. 註解與 docstring

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

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。
• CS61A 與 John DeNero 的 《Composing Programs》，本文的寫法 來自這裡：高階函式、doctest，以及「Python 會顯示什麼？」
• Python 官方文件：random、 math、 collections.Counter。