貝氏統計進階：MCMC、階層貝氏與後驗推論

當後驗算不出來的時候，貝氏統計怎麼辦？

假設你正在分析一份來自三十所學校的學習成效資料，想為每一所學校估計「採用 AI 助教後的平均進步幅度」。你心裡有個合理的先驗信念，手上也有各校的觀測資料，於是你寫下熟悉的貝氏更新式：先驗乘以概似，正比於後驗。看起來一切就緒——直到你真的動手去算那個分母，那個把後驗正規化成機率分布的積分。它是一個橫跨三十多個未知參數的高維積分，沒有封閉形式的解，連數值積分都因為維度太高而徹底失效。

這不是特例，而是現代貝氏分析的常態。真實研究裡的模型幾乎都複雜到無法手算後驗。本文要接續〈貝氏更新〉那篇打下的基礎，補上讓貝氏統計在當代真正「可用」的計算引擎：馬可夫鏈蒙地卡羅（MCMC），以及它如何撐起階層貝氏模型——一種與多層次模型互為表裡的強大框架。

回顧：先驗、概似與後驗

貝氏推論的核心只有一條公式。給定參數 $\theta$ 與觀測資料 $D$：

$$p(\theta \mid D) = \frac{p(D \mid \theta)\, p(\theta)}{p(D)}$$

其中 $p(\theta)$ 是先驗（觀測資料之前對參數的信念），$p(D \mid \theta)$ 是概似（在某組參數下觀測到這份資料的可能性），$p(\theta \mid D)$ 是我們真正想要的後驗。分母 $p(D)$ 稱為邊際概似或證據，是把分子對所有可能的 $\theta$ 積分後得到的正規化常數：

$$p(D) = \int p(D \mid \theta)\, p(\theta)\, d\theta$$

由於 $p(D)$ 不依賴 $\theta$，我們常寫成正比關係：

$$p(\theta \mid D) \propto p(D \mid \theta)\, p(\theta)$$

也就是「後驗正比於概似乘先驗」。這個正比式的形狀我們算得出來，麻煩的永遠是那個積分常數 $p(D)$。

為什麼後驗常常沒有解析解

在入門例子裡，後驗往往有漂亮的封閉形式——這通常是因為我們刻意挑了共軛先驗（稍後在深入段詳談）。但只要模型稍微貼近現實，三個麻煩就會一起浮現。

第一，參數一多，積分維度就爆炸。三十所學校各有自己的平均效應，再加上控制變異程度的超參數，未知量輕易超過三十個。正規化常數是一個三十多維的積分，沒有解析解，數值網格積分也不可行——若每個維度切十格，總格點數是 $10^{30}$，宇宙年齡也算不完。

第二，先驗與概似不共軛。一旦你用了邏輯斯迴歸的概似配上常態先驗，或在模型裡塞進階層結構，分子的代數形式就不再對應任何已知分布，積分自然寫不出來。

第三，我們其實不需要整個後驗的函數式。實務上我們要的是後驗的摘要：平均數、標準差、可信區間、某個事件的機率。如果能直接從後驗分布裡「抽樣」，再對樣本取平均、取分位數，就能繞過那個討厭的積分常數。這正是 MCMC 的切入點。

馬可夫鏈蒙地卡羅（MCMC）的核心構想

蒙地卡羅方法的精神是：與其解析地計算一個期望值，不如從分布中抽出大量樣本，用樣本平均來逼近。若我們能從後驗 $p(\theta \mid D)$ 抽出 $\theta^{(1)}, \theta^{(2)}, \dots, \theta^{(N)}$，那麼任何後驗摘要都可以用樣本估計：

$$\mathbb{E}[g(\theta) \mid D] \approx \frac{1}{N}\sum_{i=1}^{N} g\!\left(\theta^{(i)}\right)$$

問題是：高維後驗根本沒辦法直接抽樣。MCMC 的巧思在於——不直接抽樣，而是建構一條馬可夫鏈，讓這條鏈在參數空間裡遊走，使得它的長期停留分布（穩定分布）恰好等於後驗。只要鏈跑得夠久，沿途記錄下來的位置就近似於從後驗抽出的樣本。

關鍵在於：整個過程只需要算到「後驗的正比式」就夠了。因為 MCMC 比較的永遠是兩個位置的後驗「比值」，那個未知的正規化常數 $p(D)$ 在比值裡上下相消，徹底消失。這就是 MCMC 能繞開高維積分的根本原因。

Metropolis–Hastings 的直覺

最基礎的 MCMC 演算法是 Metropolis–Hastings。想像參數空間是一片地形，後驗密度的高低就是地形的海拔——後驗機率大的地方是高地，我們希望鏈多停留在那裡。流程如下：

1. 從目前位置 $\theta_{\text{current}}$ 出發，依某個提議分布隨機提出一個候選位置 $\theta_{\text{prop}}$（例如在附近隨機跳一步）。
2. 計算接受比值（以對稱提議分布為例）：

$$r = \frac{p(\theta_{\text{prop}} \mid D)}{p(\theta_{\text{current}} \mid D)} = \frac{p(D \mid \theta_{\text{prop}})\, p(\theta_{\text{prop}})}{p(D \mid \theta_{\text{current}})\, p(\theta_{\text{current}})}$$

注意正規化常數已經消掉了。

3. 若 $r \geq 1$（候選地海拔更高），就接受，移動過去；若 $r < 1$，則以機率 $r$ 接受、以機率 $1-r$ 留在原地。

「往高處走，但偶爾也接受往低處走」這條規則，讓鏈既會集中在高機率區，又不會困死在局部高點。長期下來，鏈停留在各處的頻率，正比於該處的後驗密度。

Gibbs sampling 的直覺

當參數有很多個時，Gibbs sampling 提供另一條捷徑。它不一次更新所有參數，而是輪流逐一更新每個參數，每次都從「固定其他參數時，這個參數的條件後驗」中直接抽樣：

$$\theta_1^{(t+1)} \sim p(\theta_1 \mid \theta_2^{(t)}, \theta_3^{(t)}, \dots, D)$$ $$\theta_2^{(t+1)} \sim p(\theta_2 \mid \theta_1^{(t+1)}, \theta_3^{(t)}, \dots, D)$$

依此類推，跑完一輪算一次迭代。Gibbs 的好處是不需要調提議分布、也沒有「拒絕」這回事——只要每個條件後驗都能直接抽樣（這在共軛結構下常成立），它就非常有效率。階層模型的條件後驗往往恰好是共軛的，這也是 Gibbs sampling 與階層貝氏天生契合的原因。Gibbs 其實可視為 Metropolis–Hastings 的一個特例，其接受機率恆為 1。

收斂診斷：怎麼知道鏈跑夠了

MCMC 有一個無法迴避的問題：理論保證鏈「最終」會收斂到後驗，但沒告訴你要跑多久。剛開始那段鏈還停留在起點附近、尚未進入後驗的高機率區，這段稱為 burn-in（暖機期），通常要丟掉。我們需要診斷工具判斷鏈是否已經收斂、混合是否良好。

Trace plot（軌跡圖）：把每個參數的取樣值對迭代次數畫出來。理想的軌跡圖看起來像一條穩定的、毛毛蟲般上下抖動的橫帶——表示鏈在後驗周圍均勻遊走、沒有趨勢。若軌跡圖呈現緩慢漂移、或長時間卡在某個區域，就是還沒收斂的警訊。

R-hat（$\hat{R}$，又稱 Gelman–Rubin 統計量）：這是最常用的量化指標。做法是同時跑多條起點不同的鏈，比較「鏈與鏈之間的變異」對上「每條鏈內部的變異」。若所有鏈最終都探索同一個後驗分布，這兩種變異應該趨於一致，$\hat{R}$ 會接近 1。實務上常用的門檻是 $\hat{R} < 1.01$；若明顯大於 1（例如 1.1 以上），代表不同的鏈還在探索不同區域，尚未收斂，結果不可信。

另一個常看的是有效樣本數（ESS）。MCMC 樣本彼此相關（這一步的位置高度依賴上一步），所以 $N$ 個樣本所含的獨立資訊量遠少於 $N$。ESS 估計「相當於多少獨立樣本」，太低就代表鏈混合不佳，需要跑更久或改進取樣器。

階層貝氏與部分匯集

回到開頭那三十所學校。對「估計各校效應」這件事，有兩種極端做法。

完全不匯集（no pooling）：每所學校各自估計，互不參考。問題是，只收了二十份樣本的小學校，估計值會非常不穩，極端值多半是雜訊。

完全匯集（complete pooling）：把所有學校混成一鍋，估一個共同的平均效應。問題是，這抹平了學校之間真實存在的差異。

階層貝氏提供了優雅的第三條路——部分匯集（partial pooling）。它假設各校的效應 $\theta_j$ 並非彼此獨立，而是來自一個共同的母體分布：

$$\theta_j \sim \mathcal{N}(\mu, \tau^2), \qquad j = 1, \dots, 30$$

這裡 $\mu$ 是全體學校的整體平均效應，$\tau$ 控制學校之間的離散程度，兩者都是待估的超參數，並各自有自己的先驗（這就是「階層」的由來：參數的先驗本身又有參數）。

這個結構帶來一個漂亮的後果：每所學校的後驗估計，會自動朝整體平均 $\mu$ 收縮（shrinkage），而收縮的力道取決於該校資料量。資料多、證據強的學校，後驗幾乎不被拉動；資料少、雜訊大的學校，則被往整體平均拉得多一些。換句話說，小學校「借用」了其他學校的資訊來穩定自己的估計——這正是部分匯集的核心價值，也是它在預測上常勝過兩種極端做法的原因。

帶數字的小範例：收縮如何運作

考慮一個簡化情境：某校觀測到的平均進步分數為 $\bar{y}_j = 12$ 分，而全體三十校的整體平均 $\mu = 7$ 分。假設該校觀測的抽樣變異數為 $\sigma^2/n_j = 9$（標準誤 3 分），學校間變異 $\tau^2 = 9$。

在常態—常態階層模型下，該校的後驗平均是觀測值與整體平均的加權平均，權重為收縮因子 $B$：

$$B = \frac{\sigma^2/n_j}{\sigma^2/n_j + \tau^2} = \frac{9}{9 + 9} = 0.5$$

後驗估計為：

$$\hat{\theta}_j = (1 - B)\,\bar{y}_j + B\,\mu = 0.5 \times 12 + 0.5 \times 7 = 9.5$$

原始觀測值 12 分被收縮到 9.5 分——往整體平均拉近了一半。

現在比較兩所資料量不同的學校。若某大校樣本很多，$\sigma^2/n_j = 1$（標準誤 1 分），則：

$$B = \frac{1}{1 + 9} = 0.1, \qquad \hat{\theta}_j = 0.9 \times 12 + 0.1 \times 7 = 11.5$$

只被拉動 0.5 分。資料越充足，收縮越少、越信任自己的觀測——這正符合直覺：證據強的個體應該保留自己的特色，證據弱的個體則該謙虛地向群體靠攏。

後驗預測檢查：模型配得好不好

估出後驗只是一半工作；我們還得問：這個模型產生的資料，看起來像不像真實資料？後驗預測檢查（posterior predictive check, PPC）就是做這件事的工具。

其核心是後驗預測分布——把後驗的不確定性整合進去，模擬出「未來可能觀測到的資料」：

$$p(\tilde{y} \mid D) = \int p(\tilde{y} \mid \theta)\, p(\theta \mid D)\, d\theta$$

實作上，我們從後驗抽出許多組 $\theta^{(i)}$，每一組都用來模擬一份假想資料集 $\tilde{y}^{(i)}$。接著挑一些檢定統計量（例如資料的最大值、變異數、零的比例），比較「模型模擬出的分布」與「真實資料的觀測值」。若真實值落在模擬分布的尾端、顯得格格不入，就說明模型在某方面系統性地與資料不符——可能漏掉了某種變異來源、或分布假設不對。PPC 是貝氏工作流程裡檢驗模型適配度的標準步驟，比單看點估計誠實得多。

可信區間 vs 信賴區間

貝氏與頻率學派最常被混淆的，是可信區間（credible interval）與信賴區間（confidence interval）。兩者數值上有時接近，詮釋卻天差地別。

可信區間是後驗分布的直接產物。一個 95% 可信區間 $[a, b]$ 的意思就是字面上那樣：「給定觀測資料與模型，參數落在 $[a, b]$ 之間的後驗機率是 95%。」參數在這裡被當成有機率分布的隨機量，區間是固定的，機率敘述直接落在參數上。常見做法是取後驗的 2.5% 與 97.5% 分位數（等尾區間），或取後驗密度最高的區域（HPD 區間）。

信賴區間則是頻率學派的概念。它的 95% 不是某次區間的機率，而是「程序的長期表現」：若無限次重複抽樣、每次都照同樣方法造一個區間，其中約 95% 會涵蓋那個固定但未知的真參數。對任何單一一次算出來的區間，你不能說「真值有 95% 機率在裡面」——真值要嘛在、要嘛不在。

這個差異不只是哲學品味。當你想對決策者說「這個效應有 95% 機率為正」時，只有貝氏可信區間能合法地這樣講；信賴區間在語法上不允許這種敘述。多數人對信賴區間的直覺理解，其實是可信區間的定義。

重點回顧

• 後驗常無解析解，根源是正規化常數 $p(D)$ 那個高維、非共軛的積分；MCMC 透過抽樣繞過它，且全程只需後驗的正比式。
• MCMC 建構一條穩定分布等於後驗的馬可夫鏈：Metropolis–Hastings 用「接受／拒絕」往高機率區走，Gibbs sampling 則逐一從條件後驗抽樣，特別適合階層結構。
• 收斂診斷不可省略：trace plot 看是否像穩定的毛毛蟲、$\hat{R}$ 應接近 1（門檻常取 $<1.01$）、ESS 反映自相關後的有效資訊量。
• 階層貝氏透過部分匯集讓各組估計朝整體平均收縮，收縮力道由資料量決定，小樣本群體借力於整體而更穩健。
• 可信區間對參數做直接的機率敘述，與信賴區間的「程序長期涵蓋率」詮釋截然不同，兩者不可混用。

深入探討（研究所視角）

共軛先驗的計算捷徑。 在 MCMC 普及之前，貝氏分析高度依賴共軛性。若先驗與後驗屬於同一參數族，就稱該先驗對該概似為共軛先驗。經典例子是 Beta–二項：以 $\text{Beta}(\alpha, \beta)$ 為成功率 $p$ 的先驗，觀測到 $s$ 次成功、$f$ 次失敗後，後驗恰為 $\text{Beta}(\alpha + s, \beta + f)$。更新只是把資料計數加進超參數，完全不必積分。常見的共軛對還有常態—常態（已知變異數時）、Gamma–Poisson、Dirichlet–多項等，背後的統一結構是指數族分布：對指數族概似，其充分統計量的線性更新就構成共軛先驗的封閉式後驗。共軛性在今天仍有兩重價值：其一，它讓 Gibbs sampling 的每個條件後驗都能直接抽樣，是高效取樣器的基石；其二，它提供了快速的原型與檢驗基準。但代價是先驗形狀被概似綁死——你不能因為實質知識而自由指定先驗，只能在共軛族內挑超參數。MCMC 的歷史意義，正是把貝氏從「先驗必須共軛」這道枷鎖中解放出來。

與多層次模型的連結。 階層貝氏與頻率學派的多層次模型（混合效應模型），在數學上是同一個結構的兩種詮釋。常態—常態階層模型 $\theta_j \sim \mathcal{N}(\mu, \tau^2)$，正是多層次模型裡的隨機截距：$\mu$ 是固定效應，$\theta_j - \mu$ 是各組的隨機效應，$\tau^2$ 是組間變異成分。兩者的部分匯集機制本質相同——頻率學派用最佳線性無偏預測（BLUP）得到收縮估計，貝氏則由後驗自然產生收縮。差別在於：頻率學派的混合模型通常把變異成分 $\tau^2$ 當成插入式（plug-in）點估計，忽略其不確定性；而完整的貝氏處理會給 $\tau^2$ 一個先驗，讓它的不確定性一路傳遞到各組效應的後驗裡。當組數少、$\tau^2$ 本身估得不準時，這個差異會放大——貝氏的可信區間會誠實地變寬，反映「我們連組間變異有多大都還沒把握」。也因此，當資料層級多、變異成分估計困難時，全貝氏配 MCMC 往往比頻率學派的近似更穩健。

估計理論的視角。 收縮估計並非貝氏專利，它與 James–Stein 估計量的深刻結果相通：在估計三個以上的常態平均時，把各別估計朝共同中心收縮，其總均方誤差可以一致地低於各自獨立的最大概似估計——這個反直覺的事實，正是部分匯集為何在預測上有效的理論根基。階層貝氏可視為把這種收縮以資料驅動、且帶不確定性量化的方式自動實現。

限制與假設。 上述一切建立在幾個前提上：模型設定正確（概似形式、階層結構、分布假設）、MCMC 確實收斂（診斷通過不等於萬無一失，複雜後驗仍可能藏有未被探索的模態）、先驗合理且其影響在資料充足時可被淹沒。收縮會引入偏誤以換取變異的降低——當各組真實效應差異極大、母體常態假設不成立時，過度收縮反而有害。後驗預測檢查能揪出部分模型不適配，但無法證明模型為真。所有貝氏結論都是「在此模型與此先驗下」的條件性陳述，換一組合理假設可能得到不同結果——這既是貝氏方法的誠實之處，也是使用者必須清楚溝通的責任。