劉應興的部落格

Poisson 分布，中譯：卜瓦松分布，柏松分布，波松分布。許多不同譯名，還是用英文名吧。

Poisson 分布常用在計數資料，當然計數資料不一定服從 Poisson 分布，什麼資料可以用，最好是透過資料產生程序來看。首先，它是一個點過程，是一個出生過程的產物。也就是說：計數資料是一段時間、一塊區域、或一個立體、或更抽象空間上一個「點」的數量的問題，這些點產生（出生）之後就不會消失（死亡）。其次，它要滿足三個條件：

1. 從 0 開始。
2. 獨立增量︰不相重疊的區段 (時段、地理區塊、立體區塊等) 其數量相互獨立。
3. 區段可無限細分,使每個微小區段幾乎不存在兩個點以上；或每個區段的點數都符合 Poisson 分布。

第三個條件後半段有點循環定義的感覺；不過如果是要用這三個條件描述 Poisson 分布，第三個條件我們可以用前半；只有先有 Poisson 分布再定義 Poisson process 或 Poisson field 時第三個條件會採用後半的描述。

上面三個條件，或再加上它是純出生過程的條件，加上 homogenous 條件，以時間區段描述，數學上的意思是

{N(t), t≧0} 是一個整數值隨機過程 (一組有限或無限, 可數或不可數個, 同性質如本例之「整數值」, 隨機變數, 稱之為隨機過程; 當註標是二維以上時稱隨機域, random field。) 滿足:
0. P[N(t+h)-N(t)≧0] = 1 對任意 t, 任意 h>0 成立；
1. N(0) = 0；
2. 若 [si,ti], i=1,...,k, 不相重疊, 則 N(ti)-N(si), i=1,...,k, 相互獨立；
3. P[N(t+h)-N(t)>1] = o(h), 對任意 t, 任意 h>0 成立；
4. N(t)-N(s), t>s, 的分布只和 t-s 有關。

最後一個條件是 homogenous 條件，一般常假設此條件符合，並稱之為 homogenous Poisson process；若不滿足這一條件，但滿足其他條件，則稱之為 inhomogenous Poisson process。

第三個條件的 o(h) 表示 o(h)/h → 0 當 h → 0。根據條件 0～3 可得一組遞迴微分方程式：

(d/dt) P[N(t) = 0] = -λ(t) P[N(t) = 0]
(d/dt) P[N(t) = k] = -λ(t)(P[N(t) = k] - P[N(t) = k-1]), k≧1

由第一式可解出

P[N(t)= 0] = e^{-Λ(t)}, 其中 Λ(t) = ∫_[0,t] λ(t) dt

而 k ≧ 1 的情形則可依遞迴方法得

P[N(t)= k] = [Λ(t)]^k e^{-Λ(t)}/k! 

也就是說：Poisson 過程結果（計數資料）的分布是一 Poisson 分布。若加上條件 4，λ(t) 是常數 λ，則得 P[N(t)= k] = (λt)^k e^(-λt)/k!。

由於 Poisson 過程具有獨立增量（第２項條件），因此 N(s+t)-N(s) 具有平均數 Λ(s+t)-Λ(s) （在 homogeneous 的情形是 λt）的 Poisson 分布，也就是說：不管時間起點如何，結果都是 Poisson 分布，而平均數是「到達率」 λ(t) 在計數的這時段的積分。而在以時間為標籤的 Poisson 過程，P[N(t)=0] 表示等候第一個點出現所需時間超過 t, 因為等了 t 時間還沒等到。所以在一個等候問題，如果客人到達時間不會相互有關聯，可以假設等候時間的分布是服從 P[T>t] = e^(-Λ(t)), 特例是指數分布 P[T>t] = e^(-λt)。不過等候問題的服務時間，有時候可能就不好假設是指數分布，因為除非是機械性的服務，人的服務時間可能因先前服務時間長短而在後面有所調整，也可能因等候線長度而有所調整，例如銀行服務，雖然不能加窗口，卻可能因等候線太長而有人協助做些動作以加速服務流程。也就是說：在計數過程中，欲了解其計數結果（計數資料）是否可用 Poisson 分布描述，關鍵是能否假設獨立增量。因為條件１只是計數起點的設定（從０開始），基本上與分布形狀無關；條件３則與事件點出現的形式有關，在很多情形假設短時間出現兩個事件點的機率可以忽略是合理的。條件４則是此過程到達率的描述，很多實務應用如等候問題，其到達率明顯不是均勻的，卻仍依均勻到達率 λ(t) = λ 的假設進行分析。

Poisson 過程的「可無限細分」條件也闡述了一個事實：二項分布 p 很小，n 很大，則與 Poisson 分布接近。當 [0,t] 被細分成許多段，各段基本上可認為只可能有 0, 1 事件點，1 的機率是 pi = λ ti , i = 1,...,n，而 Σpi = λt。於是我們等於有 n 個相互獨立的 Bernoulli trials, 各具成功機率 pi。這 n 個 Bernoulli 試作成功次數的分布

P[X = k] = 總共 C(n,k) 項，每一項是 k 個 pi 相乘．(n-k) 個 1-pi 相乘

如果 ti = t/n, pi = λt/n, 則上式近似 (λt)^k e^-(λt)/k!。或者，用動差母函數 (moment generating function, m.g.f., 或譯動差生成函數) 來看，X 的 m.g.f. 是

M(u) = Π[1+(λti+o(ti))(e^u-1)] 
      = e^Σln[1+(λti+o(ti))(e^u-1)] 
      = e^Σ{(λti+o(ti))(e^u-1) - (λti+o(ti))(e^u-1)^2/2 + ...}
      = e^Σ{λti(e^u-1)+o(ti)}
      = e^[λt(e^u-1)+o(1)]

即：若將 Poisson 過程無限細分，即使細分的時段不一樣長，或非一維時間情形各區塊不等大小，只要 max ti/t → 0, 則 以 Bernoulli 試作近似處理的結果，總成功數（總計數）的分布以 Poisson 分布為極限。

如果到達率 λ 的 homogeneous Poisson 過程 N(t) 的事件點出現時獨立地以機率 (p1,...,pm) 分成 m 種型態，則得 m 個同時的點過程

{ N(t) = (N1(t),...,Nm(t)); t ≧ 0 }, 

（黑體的 N 是向量，標準體的 N 是分量總和）結果 m 個分過程是相互獨立的且各 Ni(t) 均為‵ homogeneous Poisson 過程，到達率 λi = λ pi；反之，如果同時看 m 種事件點 N(t) = (N1(t),...,Nm(t)), 其到達率均為常數 λi, 而 N(t) = N1(t)+...+Nm(t) 是有一事件點算一點，顯然 N(t) 本身是一個 Poisson 過程，而 N(t) 相當於 N(t) 在每一事件點出現時以機率 pi = λi/λ 區分屬於哪一類事件點。這個關係直觀上不難理解，數學證明細節上有些煩頊，但也不難。事實上在二項過程（二項獨立試作）也如此，如果每一次試作成功再隨機地分成 m 種結果，就得到 m 個相互獨立的成功數各具二項分布；反之 m 種成功彙總成一個成功，結果就像每次試作成功隨機地歸屬到 m 種特定成功種類之一一樣。事實上在 inhomogeneous 結果也成立，只要 λi(t)/λ(t) = pi 與 t 無關。

Ｐoisson 分布隨機變數無限可分，如常態分布一般。因為 Poisson 分布或其隨機變數就是 Poisson 過程的結果，前面我們等於應用了無限可分性確立了二項分布與 Poisson 分布的關係。而無限可分性同時也暗示了一個性質：當 Poisson 的平均數 λ 足夠大時，可以用常態近似。當然無限可分性不是中央極限定理成立的必要條件，相加性 ( additive ) 就夠充分了，如二項分布、卡方分布等。不過無限可分性讓 λ 不必受限，可以是任何正數。

可能計數資料 X 本質上是平均數 λt 的 Poisson 分布，但 λ 又是隨機的，因此 X 的實際分布不是 Poisson, 這種情形有一個特別的名詞：混合 Poisson 分布 (mixed Poisson distribution 或 Poisson mixture model)。假設給定 Λ = λ ( 此處 Λ 與前面 Λ(t) 意義不同，是代表一個隨機變數，與 t 無關）時 X 的分布是 Poisson, 平均數 λ；但 Λ 本身是隨機的，具有平均數 μ，變異數 υ，則

E[X] = E[E[X|Λ]] = E[Λ] = μ
Var[X] = E[Var[X|Λ]] + Var[E[X|Λ]] = E[Λ] + Var[Λ] = μ + v

本來 Poisson 分布平均數等於變異數，這算是 Poisson 分布的特色，但當 Poisson 的平均數不是確定值而是隨機變數時，變異數超過平均數，這種現象也稱這計數變量有 overdispersion。在常態分布、二項分布等也有 mixture 的問題，事實上 t 分布就是常態分布的一種混合分布，因為 t =z/√(X/ν) 分子是標準常態變量，而分母根號內是一個與 z 獨立的卡方變量除以其自由度。

最常被提及的 Poisson mixture model 是 Poisson-gamma 模型：

X | Λ ～ Poisson(Λ), Λ ～ gamma(α,β)

因此得 計數資料 X 的邊際 p.m.f. 為

P[X = k] = ∫_(0,∞) (λ^k e^(-λ)/k!) λ^(α-1)e(-λ/β)/(Γ(α)β^α) dλ
  = Γ(k+α)β^k/[k!Γ(α)(β+1)^(k+α)]

相當於成功數 α,  成功機率 1/(1+β) 的負二項分布‵ (negative binomial distribution)。其平均數 αβ, 變異數 αβ(1+β)。

若 X|Λ=λ 是 Poi(λt) 而非 Poi(λ), 上面的結果需要調整：一般情形

E[X] = μt ,  Var[X] =  μt + vt^2

而 Poisson-gamma 模型得到的邊際分布 p.m.f. 是

P[X = k] = Γ(k+α)(βt)^k/[k!Γ(α)(βt+1)^(k+α)]

平均數 αβt, 變異數 αβt(1+βt)。

實務上我們在不同地區，不同情況下取得計數資料樣本，可能我們會配適對數線性模型，以幾個我們認為重要的變數解釋到達率的差異：

Yi ～ Poi(μi)
log(μi) = log(ti) + β0 + β1*xi + β2*zi

其中 yi 是觀測到的計數資料，ti 是該計數資料的觀測範圍（地區大小、時段長度之類，直接影響計數數量，在迴歸模型中此例 log(ti) 是已知的 offset）， xi, zi 是影響到到達率或發生率 λ 的解釋變數，也許影響的方式不是像這裡舉例的，對 log(rate) 是線性關係，但不管怎樣，我們通常還是會把它轉成線性形式。但如果有 overdispersion, 你將會發現：不管你怎樣努力找齊可能的解釋變數，模型的殘差就是偏高，描述模型配適是否良好的卡方值居高不下。當然，這其實就是沒找到所有重要的解釋變數，或也許那是難以測量甚至無法測量的，也就是說除去 offset, x, z 之外還有未知的、隨機的因素影響發生率 λ，使得看到的 Y 的變異超出 Poisson 變量應有的 Var(Yi) = E[Yi]。處理 overdispersion 的問題，最直接的當然是找出 mixing 機制，或還有重要解釋變數，或是隨機的 mixing。如果知道混合機制，例如 gamma 分布，我們可改用相應的邊際分布，如負二項來進行對數線性模型或其他廣義線模的分析。另外，一個簡單的調整方法是用最完整的模型，其卡方統計量比自由度大太多，即暗示 overdispersion 的存在。將 χ^2/df 當做應調整變異數倍數，各參數估計的標準誤（差）乘以此倍數，重新計算其 t 統計量，重新計算 p 值。

在實際的計數資料，有時我們會發現 0 異常地多，以致整體並不符合 Poisson 分布，但在非 0 部分又發現其實和 Poisson 很像，例如保險索賠的情形可能就是這樣。這稱為 zero-inflated model, 首先是由 Lambert 在 1992 提出 Poisson 迴歸的情形；1994, Greene 提出負二項分布的零膨脹模型；2000 Daniel 提出了零膨脹二項模型。網路搜尋 "zero-inflated" 可以查到不少資料，如果使用 R 語言做資料分析，

https://fukamilab.github.io/BIO202/04-C-zero-data.html

應可一看，裡面有 R 語言應用示例。其實用 R 語言做統計計算一般不難，主要就是查找所要的分析有什麼函數，在哪個 package，然後下載 package，並查尋函數用法（有哪些鄧引數、什麼作用），目前我對這些還沒實作過，等需要時再查看應該就可以了；前面關於 overdispersion，甚至 GLM 本身也一樣，用 R 語言想來會比其他語言方便——最重要的是它完全免費，不像一些套裝軟體，要不就是大機構才付得起價錢，要不就是盜版、破解版。像我這樣個人，又沒門路盜版沒本事破解的，只能尋求免費軟體了。

 與 Poisson 分布相關的還有一種分布，稱為「複合 Poisson 分布 (compound Poisson distribution)」。複合 Poisson 分布其實是一些獨立同分布隨機變數和的分布：

N ～ Poi(λ),
X1, ..., Xn, ... 獨立同分布, 且與 N 獨立，
Y= Σ{i=1,...,N} Xi

放在 Poisson 過程 N(t) 則形成複合 Poisson 過程 Y(t) = Σ{i≦N(t)} Xi。例如假設一個團體各成員的收入 Xi 是 i.i.d., 而這團體成員人數可用 Poisson 分布描述，則這團體成員總合收入就是一個複合 Poisson 分布。保險公司某種保險每年理賠件數可用 Poisson 描述，而每件理賠申請的理賠人額是 Xi, 獨立同分有，則保險公司每年在這項保險的理賠總金額就是複合 Poisson 變量。

首先我們來看複合 Poisson 分布的期望值和變異數：Xi 同分布故以 X 代。 

E[Ｙ] = E[E[Ｙ|Λ]] = E[N E[X]] = λ E[X]

Var[Y] = E[Var[Y|Λ]] + Var[E[Y|Λ]] 
       = E[Λ Var(X)] + Var[Λ E[X]] 
      = λ Var(X) + λ(E[X])^2 = λ E[X^2]

如果 X 的特性函數 ( ch. f. ) 是 φ(t), 則 Y 的 ch.f. 是

E[e^(itY)] = E[ E[e^(itY)|N] ] = E[ (E[e^(itX)])^N ]
   = E[ (φ(t))^N] = E[ e^{log(φ(t)) N} ]
   = e^{λ[e^log(φ(t)) - 1} = e^{λ(φ(t)-1)}

如果 X 的 m.g.f. 存在，為 m(t), 則 Y 的 m.g.f. 也存在：

M(t)  = E[e^(tY)] = e^[ λ (m(t) - 1) ]

不過通常複合 Poisson 分布的明確形式並不易得，雖然理論上可由 ch.f. 的反算公式得出機率甚至算出 p.d.f.，但實際上可能不容易計算。不過，就統計分析而言，很多時候我們並不需要真正知道群體分布的明確模樣；如果要做模擬，有 X 的分布和 N 的分布就足夠了。

當然關於 Poisson 分布還有很多可以談，不過還是在這裡終止吧。