八爪 Paul 與未來哥番外篇：Bin(n,p) 的 95% C.I.

圖一：Clopper-Pearson（紅）及 The Suffocated（藍）的 95% C.I. 的覆蓋機率；此處 n=8

上一篇蒙 匿名讀者《東南西北》網誌的宋以朗先生賜教，指出我構作的 Bin(n,p) 95% C.I. 與統計學上的所謂 Clopper-Pearson interval 不同。Clopper-Pearson interval 並非構作二項分佈信賴區間的唯一方法，所以不同於它並非稀奇事，不過前文中我構作 C.I. 時犯了一個錯誤，所以本文將一併解釋。

為二項式分佈找一個信賴區間，貌似簡單，其是統計學中一個著名難題。文獻中最早的構作方法，出現於 C.J. Clopper 及統計學家 E.S. Pearson 合著的一篇 1934 年文章^[1]之中。其作法十分簡單：若觀察得到 X = x，那麼所需的 (1-α)×100% C.I.，就是一個區間 (a,b) ，其中端點 a 及 b 是等式

Pr(X ≤ x | p = a) = α/2 及 Pr(X ≥ x | p = b) = α/2

的解。這種構作方式，與我於前篇提及，用 hypothesis testing 的角度去構作的方式，基本精神相同，只是我先前的做法在 X=n 的時候，考慮的是 one-sided test，而 Clopper-Pearson 考慮的，永遠都是 two-sided test。

讀者可能會想，是不是因為我從 one-sided test 出發，所以出事？確實有點關係，不過原因沒那麼單純。詳細情形容後再說，暫時且繼續談 Clopper-Pearson interval。

Clopper-Pearson interval 又通稱為 exact interval，原因是它保證有不少於 95% 的覆蓋機率 (coverage probability)。若你未聽過覆蓋機率這個詞語，就表示你念大學時很可能並非專攻統計學。且讓我們複習一下甚麼叫 95% C.I.。如果 X 是一個隨機變數，其分佈取決於參數 p，那麼一個 X 的 95% C.I. 是一個隨機區間 (a(X), b(X))，其端點 a(X), b(X) 為只取決於 X 而非 p 的函數，而且：

對任何固定的 p 而言，都有 Pr{ (a(X), b(X)) ∋ p } = 0.95。^[2]

舉個例，當我們說常態分佈 X~N(μ,σ²) 的 two-sided 95% C.I. 是 (X - 1.96σ, X + 1.96σ)，意思就是 a(X) = X - 1.96σ，而 b(X) = X + 1.96σ。 要留意 X 是隨機變數，所以 (a(X), b(X)) 也是隨機的。當我們有了 (X - 1.96σ, X + 1.96σ) 這個 95% C.I.，我們可以觀察 X 的一個 outcome x，從而計算該 C.I. 的一個實現 (realization)，而這個實現的區間 (x - 1.96σ, x + 1.96σ)，也就是一般中學教科書所稱的 95% C.I.。實際上，這種懶惰的說法混淆了「隨機區間」跟「隨機區間的實現」兩者，不過下面我會依然從俗，不作區分。

Pr{ (a(X), b(X)) ∋ p } 的數值，就是隨機區間 (a(X), b(X)) 的所謂覆蓋機率，而所謂 95% C.I.，實際上就是一個覆蓋機率等於 0.95 的隨機區間。 

在上面常態分佈的例子，不難驗證，不管 μ 的數值為何，Pr{ (X - 1.96σ, X + 1.96σ) ∋ μ } = Pr(μ - 1.96σ < X < μ + 1.96σ) = 0.95 這條等式都成立。然而，對許多隨機分佈，特別是離散分佈來說，覆蓋機率恰好等於 0.95 的隨機區間未必存在。不巧地，二項式分佈正是這樣的情形。這其實不難解釋。如果 X ~ Bin(n,p)，而我們記 X = 0, 1, 2, ... n 的每一種情形的 likilihood function 為

f_x(p) = C(n,x) p^x (1-p)^n-x,  x=0,1,...,n,

而且將 (a(x), b(x)) 覆蓋 p 的真假值寫成 p 的 membership function:

δ_x(p) = 1 if a(x)< p < b(x), or 0 otherwise.

那麼 (a(X), b(X)) 的 coverage probability 就是

Pr{ (a(X), b(X)) ∋ p } = f₀(p) δ₀(p) + f₁(p) δ₁(p) + ... + f_n(p) δ_n(p).

很明顯，每個 likelihood function f_x(p) 都是連續函數：

圖二

然而，除非 (a(x), b(x)) = [0,1] 或者它是空集，否則 δ_n(p) 對 p 並不連續，因此覆蓋機率 Pr{ (a(X), b(X)) ∋ p } 亦非連續函數，更不可能等於 0.95 這個 constant (hence continuous) function。

舉例說，在八爪魚 Paul 的情形 (n=8)，可以計算 Clopper-Pearson interval 如下：

[a(0), b(0)) = [0, 0.3694),
(a(1), b(1)) = (0.0032, 0.5265),
(a(2), b(2)) = (0.0319, 0.6509),
(a(3), b(3)) = (0.0852, 0.7551),
(a(4), b(4)) = (0.1570, 0.8430),
(a(5), b(5)) = (0.2449, 0.9148),
(a(6), b(6)) = (0.3491, 0.9681),
(a(7), b(7)) = (0.4735, 0.9968),
(a(8), b(8)] = (0.6306, 1.0000].

因此如圖三左所示，δ₀(p) 是一個於 (0, 0.3694) 以內為 1，以外為 0 的階梯函數 (step function)，而 f₀(p) δ₀(p) 則在 p=0.3694 處不連續（圖三右）：

圖三：（左）f₀(p), δ₀(p) 及（右） f₀(p) δ₀(p)

對其他 f_x(p) δ_x(p) 也有同樣情形：

圖四：f_x(p) δ_x(p), x = 0,1,...,8.

所以覆蓋機率 f₀(p) δ₀(p) + f₁(p) δ₁(p) + ... + f₈(p) δ₈(p) 並不連續，更非恆等於 0.95：

圖五：Bin(8,p) 的 95% Clopper-Pearson interval 的 coverage probability

換句話說，雖然 Clopper-Pearson interval 俗稱為 "exact interval"，但它一點也不 exact。事實上，在 n=8 的情形，它名義上為 95% 的 C.I. 實際上的覆蓋機率最少也差不多等於 0.97，超出 0.95 兩個百分點。

如果是連續分佈，它的覆蓋機率不是採取 f₀(p) δ₀(p) + f₁(p) δ₁(p) + ... + f_n(p) δ_n(p) 這種有限和的形式，而是採取  ∫ f_x(p) δ_x(p) dx 這種積分形式。因此，儘管每個 δ_x(p) 並不連續，但積分本身卻有機會連續，亦有機會可以設計成恆等於 0.95。

對離散分佈來說，往往不可能構作 exact C.I.，因此要退而求其次，但是放棄了覆蓋機率 = 0.95 這個目標之後，要換上一個甚麼目標，就要視乎情況。有些人希望保證覆蓋機率不少於0.95，Clopper-Pearson interval 則正正為此而設。

而我於前文所作嘗試，也是希望保證覆蓋機率不少於 0.95，可惜構作時卻犯了錯誤。這個做法將 x=n 的 95% C.I. 定為 (a,1]，其中 a 為 Pr(X ≤ x | p = a) = 0.05 的解。類似地，x=0 的 95% C.I. 被定為 [0,b)，其中 b 為 Pr(X ≥ x | p = b) = 0.05 的解。

原則上，這種取單邊走向的構作方式並沒有錯，對某些 n 來說，這種做法有時也可以獲得覆蓋機率不少於 0.95，同時區間闊度又比 Clopper-Pearson interval 窄的 95% C.I.。例如當 n=1，用前文所述方式得到的 95% C.I. 為
X=0: [0, 0.95),
X=1: [0.05, 1],
而相應的 Clopper-Pearson 95% C.I. 為
X=0: [0, 0.975),
X=1: [0.025, 1],
所以 Clopper-Pearson interval 較差。然而我這種構作方式卻不能保證對任何 n 都有不少於 95% 的覆蓋機率。當 n=8，就正好出事：我對 X=8 所作的 C.I. 為 (0.6877,1]，所以，依照同樣原理，對 X=0 的 C.I. 就是 [0, 1-0.6877) = [0, 0.3123)。因此，f₀(p) δ₀(p) 及 f₈(p) δ₈(p) 有如下圖所示：

圖六

當 p=0.6877，剛好有 f₀(p) δ₀(p) + f₈(p) δ₈(p) = 0。因此，若要覆蓋機率在 p=0.6877 的時候不小於 0.95，我們必須有 f₁(p) δ₁(p) + ... + f₇(p) δ₇(p) ≥ 0.95。然而，這是絕不可能的，因為即使 δ₁(p), ..., δ₇(p) 全部等於 1，f₁(0.6877) + ... + f₇(0.6877) 仍不及 0.95 這麼大！

雖然我於前文的構作錯誤，但這並不是說 Clopper-Pearson interval 是唯一保證覆蓋機率達到 95% 的構作方式。事實上，只要小心平衡各個 (a(x), b(x)) 的端點的先後次序及整體分佈，我們仍可獲得覆蓋機率不少於 0.95，同時區間闊度又比 Clopper-Pearson interval 窄的 95% C.I.。例如以下就是 the suffocated 用某種貪婪算法 (greedy algorithm)，對 n=8 所構作的 95% C.I.：

X=0: [0, 0.3156),
X=1: (0.0063, 0.5001),
X=2: (0.0464, 0.6845),
X=3: (0.1111, 0.7108),
X=4: (0.1929, 0.8071),
X=5: (0.2892, 0.8889),
X=6: (0.3155, 0.9536),
X=7: (0.4999, 0.9937),
X=8: (0.6844, 1.0000].

圖一是這個 C.I. 及 Clopper-Pearson intervals 的覆蓋機率，讀者可見我作的 C.I. 的覆蓋機率，比 Clopper-Pearson intervals 的覆蓋機率更接近 0.95。我的 C.I. 的平均覆蓋機率 0.9741，也比Clopper-Pearson interval 的 0.9860 低。此外，我的 C.I. 除了在 X=2 及 X=6 時比 Clopper-Pearson intervals 略長，其餘的都比 Clopper-Pearson intervals 短。

所以，要構作比 Clopper-Pearson "exact" C.I. 更 exact 的 C.I.，是可能的。然而這需要非常小心的電腦計算，而 Clopper-Pearson interval 可貴的地方，在於它可以用簡單而且直觀的方法，就能保證獲得覆蓋機率不小於名義信賴度的 C.I.。

鳴謝：多謝 匿名讀者 宋先生於前文留言，令這系列本來是吹水的文章被迫變成認真的學術討論。

註：
[1] C.J. Clopper and E.S. Pearson (1934), The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika, 26(4): 404-413.
[2] 寫成 (a(X), b(X)) ∋ p 而不是 p ∈ (a(X), b(X))，是為了強調有關的事件為「隨機區間 (a,b) 覆蓋固定的點 p」，而非「隨機的 p 落在固定區間 (a,b) 當中」。）