3.2 卜松分配#

本節目標#

理解卜松分配是「單位時間／面積內稀有事件次數」的模型
推導二項 → 卜松的極限（ $n \to \infty, p \to 0, np = \lambda$ ），並用程式驗證收斂
重現統計史上最著名的資料集：普魯士騎兵馬踢致死（Bortkiewicz, 1898）
會用 r.dpois / r.ppois / r.qpois 與 g.rpois，解決容量規劃類問題（網站每分鐘請求數）

問題場景：數「次數」而不是數比例#

3.1 節的二項分配回答「 $n$ 次試驗裡成功幾次」。但很多場景連 $n$ 都說不清：

一小時內客服接到幾通電話？（ $n$ = 一小時內「可能來電的瞬間」？無從定義）
一頁書上打錯幾個字？
一年內某段鐵路出幾次軌？

這類問題的共同特徵：事件稀疏地發生，我們只能數個數，不知道分母。法國數學家 西莫恩・德尼・卜松（Siméon Denis Poisson, 1781–1840） 在 1837 年的《刑事與民事判決中的機率研究》中給出了解答：讓二項分配的試驗數趨於無限、單次機率趨於零、但乘積固定，極限就是一條新的分配。有趣的是，卜松本人並沒有看出它的應用價值——是六十年後的俄國統計學家 Ladislaus Bortkiewicz 用它解釋了一組兵營死亡資料，這個分配才真正「活」了過來（見下文 1898 年的歷史事件）。

理論推導#

二項的極限： $n$ 大、 $p$ 小、 $np = \lambda$ 固定#

設 $Y_n \sim \mathrm{Bin}(n, p_n)$ ，其中 $n p_n \to \lambda$ （ $p_n$ 隨 $n$ 變小）。對固定的 $k$ ：

$$ P(Y_n = k) = \binom{n}{k} p_n^k (1-p_n)^{n-k} \;\xrightarrow[n\to\infty]{}\; \frac{\lambda^k e^{-\lambda}}{k!}

$$ 推導只需三步（代 $p_n = \lambda / n$ ）：

$\binom{n}{k} / n^k = \frac{n(n-1)\cdots(n-k+1)}{k!\, n^k} \to \dfrac{1}{k!}$ （分子有 $k$ 個因子，每個除以 $n$ 後趨近 1）
$p_n^k = (\lambda/n)^k \to \lambda^k / n^k$ ，與第 1 步合併得 $\lambda^k / k!$
關鍵的極限 $(1 - \lambda/n)^n \to e^{-\lambda}$

直覺： $n$ 個人各自以極小機率獨立出事，合起來就是「單位內平均發生 $\lambda$ 次」的稀疏事件流。這就是卜松（Poisson）分配，記作 $Y \sim \mathrm{Pois}(\lambda)$ ：

$$ P(Y = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \qquad k = 0, 1, 2, \ldots

$$ 它只有一個參數 $\lambda$ （單位內的平均次數），期望值與變異數都等於 $\lambda$ （3.5 節推導）。「均值 = 變異數」是卜松分配的指紋，第 14 章檢定離散資料時會用到。

「稀有」要稀有到什麼程度？#

實務經驗法則： $n \ge 20$ 且 $p \le 0.05$ 時， $\mathrm{Bin}(n,p) \approx \mathrm{Pois}(np)$ 已相當準。下面用程式驗證收斂速度。

歷史事件（1898）：普魯士騎兵的馬踢致死資料#

波蘭裔統計學家 Ladislaus Bortkiewicz（1868–1931） 在 1898 年出版《小數法則》（Das Gesetz der kleinen Zahlen），分析了普魯士軍隊 1875–1894 年間 14 個軍團、共 280 個「軍團-年」 中被馬踢死的士兵人數。20 年間共 196 人死於馬踢，平均每軍團每年 $\lambda = 196/280 = 0.7$ 人。他把 280 個軍團-年的死亡人數做成次數表，與 $\mathrm{Pois}(0.7)$ 的理論機率對照——兩者驚人地吻合，這成了卜松分配最經典的驗證案例，之後放射性衰變、電話呼叫量、V2 飛彈轟炸倫敦（1946，Clarke）的建模全都沿用同一條思路。

動手做：用 rstatistics 收斂驗證與經典資料#

收斂：二項 → 卜松#

固定 $\lambda = 10$ ，取三組 $(n, p)$ ，逐點比較 $P(Y=k)$ ：

import rstatistics as r

for n, p in [(100, 0.1), (1000, 0.01), (10000, 0.001)]:
    lam = n * p                      # 固定為 10
    diffs = [abs(r.dbinom(k, n, p) - r.dpois(k, lam)) for k in range(10)]
    print(f"n={n:6d}  p={p:<6}  max|差| = {max(diffs):.2e}")

n=   100  p=0.1     max|差| = 5.31e-03
n=  1000  p=0.01    max|差| = 5.03e-04
n= 10000  p=0.001   max|差| = 5.01e-05

$n$ 每放大 10 倍，最大誤差就縮小約 10 倍——收斂是線性的（誤差 $\approx \lambda^2/(2n)$ 的量級）。以 $n=10000, p=0.001$ 而言，二項與卜松的前幾位機率在第 4 位小數已分不出來。結論： $p$ 只要小於 0.01、 $n$ 上百，就可以放心用卜松。

重現 Bortkiewicz 資料表#

import rstatistics as r

# 普魯士騎兵 1875-1894：死亡 k 匹人數（0 到 4 匹馬踢死）的軍團-年數
deaths = [0, 1, 2, 3, 4]
counts = [144, 91, 32, 11, 2]
tot = sum(counts)                                # 280 個軍團-年
lam = sum(k * c for k, c in zip(deaths, counts)) / tot   # 0.7

print(f"死亡數  觀察次數  觀察頻率   Pois(0.7) 理論")
for k, c in zip(deaths, counts):
    print(f"  {k}      {c:>4}    {c/tot:.4f}     {r.dpois(k, 0.7):.4f}")
print(f"  5+      {0:>4}    0.0000     {1 - r.ppois(4, 0.7):.5f}")

死亡數  觀察次數  觀察頻率   Pois(0.7) 理論
  0       144    0.5143     0.4966
  1        91    0.3250     0.3476
  2        32    0.1143     0.1217
  3        11    0.0393     0.0284
  4         2    0.0071     0.0050
  5+        0    0.0000     0.00079

觀察頻率與理論機率的偏差都在 0.02 上下——對 280 個樣本而言，這正是蒙特卡羅／抽樣誤差的量級（每格標準誤 $\approx \sqrt{p(1-p)/280} \approx 0.03$ ，見 3.1 節）。127 年前的資料與 1837 年的公式就這樣對上了；到了 20 世紀，卡方適合度檢定（第 14 章）可以給這個「吻合」一個正式的 p 值。

模擬對照：`g.rpois`#

g = r.set_seed(42)
n = 200000
xs = g.rpois(n, lambda_=0.7)

cnt = {}
for x in xs:
    cnt[int(x)] = cnt.get(int(x), 0) + 1
for k in range(6):
    print(k, round(r.dpois(k, 0.7), 5), round(cnt.get(k, 0) / n, 5))

fx = list(map(float, xs))
print(r.mean(fx), r.var(fx))   # 理論皆為 0.7

0 0.49659 0.49703
1 0.34761 0.34649
2 0.12166 0.12172
3 0.02839 0.02891
4 0.00497 0.00511
5 0.0007 0.00069
0.700865 0.7033267684087243

模擬均值 0.7009、變異數 0.7033，兩者都貼著理論值 $\lambda = 0.7$ ——「均值 = 變異數」在模擬裡也看得見。

rppois 應用：網站每分鐘請求數#

網站流量是典型的卜松事件流。假設平均每分鐘 $\lambda = 8.3$ 個請求，伺服器每分鐘最多處理 11 個，請問超載（超過 11 個請求）的機率？95% 分位的尖峰流量是多少？

import rstatistics as r

lam = 8.3
print(r.ppois(11, lam))          # P(Y <= 11)：一處理得完的機率
print(1 - r.ppois(11, lam))      # 超載機率
print(r.qpois(0.95, lam))        # 95% 分位數：容量該設多少？

0.8652132378551401
0.13478676214485988
13.0

每分鐘約 13.5% 的機率超過 11 個請求（大約每 7 分鐘超載一次）；要把超載壓到 5% 以下，容量得設到 13 個/分鐘。若想壓得更低，可繼續用 r.ppois(q, lam) 掃描不同的 $q$ 找出滿足 $P(Y \le q) \ge 0.99$ 的最小容量。

與 R 對照#

dpois(0:5, 10)
# [1] 4.539993e-05 4.539993e-04 2.269996e-03 7.566655e-03 1.891664e-02
# [6] 3.783327e-02

set.seed(42)
xs <- rpois(200000, 0.7)
mean(xs); var(xs)
# [1] 0.700865
# [1] 0.7033268

ppois(5, 8.3)     # [1] 0.1652734
qpois(0.95, 8.3)  # [1] 13
1 - ppois(4, 0.7) # [1] 0.0007855354

（本節 R 輸出以 Rscript 實際執行驗證。）函數名 dpois / ppois / qpois / rpois 與參數位置完全一致，Python 版唯一差別是參數要寫 lambda_=8.3（lambda 是 Python 保留字）。

小結#

卜松分配 $\mathrm{Pois}(\lambda)$ 是「單位內稀有事件次數」的模型，pmf $= \lambda^k e^{-\lambda}/k!$ ，只有一個參數。
它是二項分配的極限： $n \to \infty$ 、 $p \to 0$ 、 $np = \lambda$ 固定； $n \ge 100$ 、 $p \le 0.01$ 時誤差已在第 3 位小數之後。
Bortkiewicz（1898）的普魯士騎兵資料（280 軍團-年、 $\lambda = 0.7$ ）與理論頻率吻合，是統計建模的歷史里程碑。
「均值 = 變異數 = $\lambda$ 」是卜松分配的指紋。
容量規劃三問：ppois 算超載機率、qpois 算容量分位數、g.rpois 模擬壓力測試。

練習#

用 $n = 20, p = 0.05$ （ $\lambda = 1$ ）與 $n = 100, p = 0.01$ （ $\lambda = 1$ ）分別計算 $P(Y=0), P(Y=1), P(Y=2)$ ，與 $\mathrm{Pois}(1)$ 對照，誤差各是多少？提示：迴圈比對 r.dbinom 與 r.dpois。
某客服中心平均每 10 分鐘 4 通電話，且來電獨立。求 10 分鐘內恰好 2 通、至少 6 通的機率。提示： $\lambda = 4$ ；「至少 6 通」用 1 - r.ppois(5, 4)。
重做 Bortkiewicz 對照表，但把理論改用你從資料估計的 $\hat\lambda$ （0.7）與固定 $\lambda = 0.65$ 各算一遍，觀察哪個擬合較好。提示： $\hat\lambda =$ 總死亡數 ÷ 280。
網站每分鐘請求 $\lambda = 20$ 。想讓每分鐘超載機率低於 1%，容量至少要幾？用 r.qpois(0.99, 20) 回答。提示：分位數即最小整數 $q$ 使 $P(Y \le q) \ge 0.99$ 。
用 g.rpois 模擬 100000 次 $\mathrm{Pois}(5)$ ，計算模擬變異數與 5 的相對誤差；再把樣本數改成 1 萬次重做，驗證誤差放大 $\sqrt{10}$ 倍。提示：標準誤 $\approx \sqrt{\lambda/n}$ （變異數的估計誤差）。