平均每小時 5 個瑕疵品：
poisson.pmf 和 cdf 各在算什麼？

先用「生產線上每小時數瑕疵品」把機率分佈、卜瓦松公式、使用條件講清楚，中間兩支動畫：一支在生產線上一小時一小時數，一支打開記憶體看 4 行程式各做了什麼。

把觀念合起來：4 行程式在記憶體裡做了什麼

SciPy 是建立在 NumPy 之上的科學計算套件，scipy.stats 是它的統計模組，裡面有一百多種機率分佈。每一種分佈都是一個物件，用法都一樣：換一種分佈，方法名稱不變。

參數 mu 就是 λ。方法的簽名是 pmf(k, *args, **kwds)，poisson.shapes 實測是 'mu'，所以 poisson.pmf(5, 5) 和 poisson.pmf(5, mu=5) 一樣。也可以先「凍結」參數：dist = poisson(5)，之後寫 dist.pmf(5)、dist.cdf(10)，不用每次都傳 5。

實驗 1 用 JavaScript 即時計算卜瓦松機率。本頁用 Node.js 跑同一段計算程式，和 SciPy 1.16.2 的 pmf、cdf、sf 比對了 λ＝0.5～15（每 0.5）× k＝0～40，共 3,690 個數字，相對誤差都在 10⁻⁹ 以內（畫面顯示到小數 4 位）。實驗 2 的每一根長條，是事先用 NumPy 2.4.3 實際抽出來的次數（和 §6 程式 3 同一份亂數）。

常見的錯誤想法

三支程式都用到 SciPy（程式 3 還用到 NumPy），要先安裝：pip install scipy（會一起裝好 NumPy；Google Colab 已經內建，直接貼上就能跑）。

程式 1：考卷的程式，再用公式對一次，列出 k＝0～10

import math
from scipy.stats import poisson

lambda_poisson = 5
print(poisson.pmf(5, lambda_poisson))                       # 考卷的程式：剛好 5 個的機率
print(5**5 * math.exp(-5) / math.factorial(5))              # 用公式自己算一次
for k in range(11):
    print(f'k = {k:>2}：pmf {poisson.pmf(k, 5):.4f}　cdf {poisson.cdf(k, 5):.4f}')

0.17546736976785068
0.1754673697678507
k =  0：pmf 0.0067　cdf 0.0067
k =  1：pmf 0.0337　cdf 0.0404
k =  2：pmf 0.0842　cdf 0.1247
k =  3：pmf 0.1404　cdf 0.2650
k =  4：pmf 0.1755　cdf 0.4405
k =  5：pmf 0.1755　cdf 0.6160
k =  6：pmf 0.1462　cdf 0.7622
k =  7：pmf 0.1044　cdf 0.8666
k =  8：pmf 0.0653　cdf 0.9319
k =  9：pmf 0.0363　cdf 0.9682
k = 10：pmf 0.0181　cdf 0.9863

程式 2：四個選項各算一次

from scipy.stats import poisson

mu = 5
print('A｜超過 5 個 P(X > 5)：', poisson.sf(5, mu))
print('B｜剛好 5 個 P(X = 5)：', poisson.pmf(5, mu), '｜少於 5 個 P(X < 5)：', poisson.cdf(4, mu))
print('D｜10 個以下 P(X ≤ 10)：', poisson.cdf(10, mu), '｜10 個以上 P(X ≥ 10)：', poisson.sf(9, mu))
print('平均、變異數：', poisson.mean(mu), poisson.var(mu))

A｜超過 5 個 P(X > 5)： 0.384039345166937
B｜剛好 5 個 P(X = 5)： 0.17546736976785068 ｜少於 5 個 P(X < 5)： 0.44049328506521257
D｜10 個以下 P(X ≤ 10)： 0.9863047314016171 ｜10 個以上 P(X ≥ 10)： 0.03182805730620481
平均、變異數： 5.0 5.0

程式 3：模擬 10 萬小時，平均發生率固定和不固定

import numpy as np
from scipy.stats import poisson

rng = np.random.default_rng(0)
hours = 100_000
steady = rng.poisson(5, size=hours)                     # 每小時的平均發生率都是 5
rate = np.where(np.arange(hours) % 2 == 0, 2, 8)        # 一小時 2、一小時 8：整體平均還是 5
unsteady = rng.poisson(rate)
for name, x in (('平均固定  ', steady), ('平均不固定', unsteady)):
    print(f'{name}：平均 {x.mean():.3f}，變異數 {x.var():.3f}，剛好 5 個的比例 {(x == 5).mean():.4f}')
print(f'卜瓦松(5) 理論值：平均 5，變異數 5，剛好 5 個 {poisson.pmf(5, 5):.4f}')

平均固定  ：平均 4.998，變異數 5.036，剛好 5 個的比例 0.1751
平均不固定：平均 5.004，變異數 13.975，剛好 5 個的比例 0.0648
卜瓦松(5) 理論值：平均 5，變異數 5，剛好 5 個 0.1755