貝氏統計進階（進階篇）：HMC 機制、模型比較與後驗預測檢驗

當後驗已經算出來：你真的「驗證」過你的貝氏模型嗎？

入門篇帶你走過了從先驗到後驗、用 MCMC（Markov Chain Monte Carlo）採樣、用階層結構借力跨群體的整套流程。假設你照做了：跑完 4 條鏈、看了 trace plot、得到一條漂亮的後驗分布。現在問題來了——這個模型「對」嗎？

很多人在這裡停下來，把後驗均值報出去就交差。但一位謹慎的貝氏分析者會繼續追問三件事：

1. 採樣到底可不可信？ 鏈收斂了不代表它探索完了整個後驗。
2. 這個模型比另一個好嗎？ 兩個競爭模型，怎麼公平比較？
3. 它對「沒看過的資料」會怎麼預測？ 後驗本身不會說話，能預測的是後驗預測分布（posterior predictive distribution）。

這篇進階篇就鎖定這三個入門少談、但決定研究成敗的環節：現代梯度式 MCMC 的內部機制、貝氏模型比較與資訊準則，以及後驗預測檢驗。我們會把數學攤開，也會帶數字算給你看。

從隨機漫步到漢米爾頓動力學：HMC 為什麼快

入門篇大概率用的是 Metropolis-Hastings 或 Gibbs 採樣。這些方法在低維度時堪用，但維度一高就災難性失效——它們本質上是「隨機漫步（random walk）」，在高維後驗的狹長山脊上像醉漢一樣慢慢挪，自相關高得嚇人。

現代工具（Stan、PyMC、NumPyro）的預設引擎是 Hamiltonian Monte Carlo（HMC，漢米爾頓蒙地卡羅），以及它的自動調參版本 NUTS（No-U-Turn Sampler）。核心想法借自物理：把要採樣的參數 $\theta$ 想成一個有位能的粒子，再給它一個輔助的「動量（momentum）」變數 $r$，讓粒子依照物理定律滑行。

我們定義位能（potential energy）為負對數後驗：

$$ U(\theta) = -\log p(\theta \mid \mathcal{D}) = -\log p(\mathcal{D} \mid \theta) - \log p(\theta) + \text{const} $$

動能（kinetic energy）取一個高斯動量，$K(r) = \tfrac{1}{2} r^\top M^{-1} r$。系統的總能量是漢米爾頓量（Hamiltonian）：

$$ H(\theta, r) = U(\theta) + K(r) $$

接著依照漢米爾頓方程式讓粒子演化：

$$ \frac{d\theta}{dt} = \frac{\partial H}{\partial r} = M^{-1} r, \qquad \frac{dr}{dt} = -\frac{\partial H}{\partial \theta} = -\nabla_\theta U(\theta) $$

關鍵在第二式：演化方向用到了 後驗的梯度 $\nabla_\theta U$。這就是 HMC 之所以快的原因——它不像隨機漫步那樣盲目試探，而是利用後驗的局部斜率資訊，沿著高機率區域有方向地滑行，一步就能走很遠卻仍維持高接受率。

實作上用 leapfrog（蛙跳）積分器離散化這條軌跡，半步更新動量、整步更新位置、再半步更新動量：

$$ r_{t+\frac{\epsilon}{2}} = r_t - \frac{\epsilon}{2}\nabla_\theta U(\theta_t), \quad \theta_{t+\epsilon} = \theta_t + \epsilon\, M^{-1} r_{t+\frac{\epsilon}{2}}, \quad r_{t+\epsilon} = r_{t+\frac{\epsilon}{2}} - \frac{\epsilon}{2}\nabla_\theta U(\theta_{t+\epsilon}) $$

走完 $L$ 步後，因為離散化有誤差，再用一個 Metropolis 步驟校正，接受機率為：

$$ \alpha = \min\!\left(1,\ \exp\big(H(\theta_t, r_t) - H(\theta^*, r^*)\big)\right) $$

理想軌跡上 $H$ 守恆，所以 $\alpha \approx 1$，幾乎每步都接受。NUTS 進一步自動決定要走幾步（$L$）——它持續延伸軌跡，直到軌跡開始「掉頭（U-turn）」折返為止，省去手動調 $L$ 的痛苦。

這對你的意義：當你看到 Stan 警告「divergent transitions（發散轉移）」，那不是小事。它代表 leapfrog 在後驗某些高曲率區域（典型如階層模型的變異數參數接近 0 處）積分失準，粒子「飛出去」了。這意味著該區域的後驗被系統性地漏採，你的推論可能有偏。標準解法是把模型改寫成 non-centered parameterization（非中心化參數化），放鬆那個導致漏斗形（funnel）幾何的相依結構。

比較模型：別再用後驗機率硬比

假設你有兩個模型 $M_1$、$M_2$。最「純」的貝氏做法是算貝氏因子（Bayes factor）：

$$ \text{BF}_{12} = \frac{p(\mathcal{D} \mid M_1)}{p(\mathcal{D} \mid M_2)}, \quad p(\mathcal{D} \mid M_k) = \int p(\mathcal{D} \mid \theta, M_k)\, p(\theta \mid M_k)\, d\theta $$

這個邊際似然（marginal likelihood，又稱 evidence）對先驗極度敏感——把一個無關緊要參數的先驗從 $N(0,1)$ 改成 $N(0,100)$，貝氏因子可以變幾個數量級。而且高維積分本身難算。所以實務上，預測導向的比較更受青睞。

現代主流是 WAIC（Widely Applicable Information Criterion） 和 PSIS-LOO（Pareto-Smoothed Importance Sampling Leave-One-Out cross-validation）。兩者都在估計同一個東西：期望對數逐點預測密度（expected log pointwise predictive density, ELPD），也就是模型對新觀測值的平均對數預測能力。

LOO 的精神是「留一交叉驗證」：對每個觀測 $y_i$，用「去掉 $y_i$ 的其餘資料」訓練，再看模型對 $y_i$ 預測得多好：

$$ \text{elpd}_{\text{loo}} = \sum_{i=1}^{n} \log p(y_i \mid y_{-i}) = \sum_{i=1}^{n} \log \int p(y_i \mid \theta)\, p(\theta \mid y_{-i})\, d\theta $$

天真做法要重跑 $n$ 次 MCMC，貴到不可行。PSIS-LOO 的巧妙在於：用同一份完整後驗樣本，透過重要性權重（importance weights）把「全資料後驗」重新加權成「去掉第 $i$ 點的後驗」：

$$ w_i^{(s)} \propto \frac{1}{p(y_i \mid \theta^{(s)})} $$

直覺是——一個對 $y_i$ 解釋力很強的參數值，在去掉 $y_i$ 之後就該被降權。原始重要性權重的尾巴常常很重、方差爆炸，於是 PSIS 用一個廣義 Pareto 分布去平滑那條尾巴，並回報形狀參數 $\hat{k}$ 作為診斷：$\hat{k} < 0.7$ 才可信，超過就代表那個資料點是模型難以涵蓋的高影響力點（influential point），需要對它真的重跑一次。

動手算一下：用 WAIC 比較兩個模型

WAIC 的公式把「擬合度」和「複雜度懲罰」拆開：

$$ \text{WAIC} = -2\Big(\underbrace{\text{lppd}}_{\text{擬合}} - \underbrace{p_{\text{WAIC}}}_{\text{有效參數懲罰}}\Big) $$

其中逐點對數預測密度（log pointwise predictive density）用 $S$ 個後驗樣本估計：

$$ \text{lppd} = \sum_{i=1}^{n} \log\left(\frac{1}{S}\sum_{s=1}^{S} p(y_i \mid \theta^{(s)})\right) $$

而有效參數數用後驗對數似然的方差估計（這很妙——模型越靈活，不同後驗樣本對同一點的似然分歧越大）：

$$ p_{\text{WAIC}} = \sum_{i=1}^{n} \operatorname{Var}_{s}\!\big[\log p(y_i \mid \theta^{(s)})\big] $$

來個迷你數值例子。假設只有 $n=3$ 個資料點，我們各自用 $S=4$ 個後驗樣本算出每點的似然 $p(y_i\mid\theta^{(s)})$：

先算 lppd（每點取似然平均、再取 log、再加總）：

• $y_1$：平均 $= (0.30+0.40+0.20+0.50)/4 = 0.35$，$\log 0.35 = -1.0498$
• $y_2$：平均 $= (0.10+0.15+0.05+0.20)/4 = 0.125$，$\log 0.125 = -2.0794$
• $y_3$：平均 $= (0.60+0.55+0.65+0.50)/4 = 0.575$，$\log 0.575 = -0.5534$

$$ \text{lppd} = -1.0498 - 2.0794 - 0.5534 = -3.6826 $$

再算 $p_{\text{WAIC}}$（每點先取 $\log$ 似然的樣本方差，再加總）。以 $y_1$ 為例，四個 log 似然是 $\log 0.30, \log 0.40, \log 0.20, \log 0.50 = -1.204, -0.916, -1.609, -0.693$，其樣本方差約 $0.118$。同理 $y_2 \approx 0.260$、$y_3 \approx 0.0103$，加總得 $p_{\text{WAIC}} \approx 0.388$。

$$ \text{WAIC} = -2(\,-3.6826 - 0.388\,) = -2(-4.0706) = 8.141 $$

WAIC 越小越好。若另一個模型 $M_2$ 算出 $\text{WAIC}=11.5$，則 $M_1$ 勝出。但別只看點估計：實務上要連同 WAIC 差值的標準誤一起看，若 $\Delta\text{WAIC}$ 小於約 2 倍標準誤，兩模型其實難分軒輊，貿然宣稱誰贏是過度詮釋。

後驗預測檢驗：讓模型自己生資料給你看

判斷模型好壞最直觀的方式之一，是問：「如果這個模型是真的，它會生出什麼樣的資料？跟我手上的真實資料像嗎？」這就是後驗預測檢驗（posterior predictive check, PPC）。

後驗預測分布是把參數的不確定性「積分掉」後，對新資料 $\tilde{y}$ 的預測：

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

操作上不必真的積分：你已有 $S$ 個後驗樣本 $\theta^{(s)}$，對每個樣本從似然模型抽一筆模擬資料集 $\tilde{y}^{(s)}$，就得到 $S$ 個「模型生成的平行世界」。接著選一個檢驗統計量（test statistic） $T(\cdot)$——例如資料的最大值、變異數、或零值比例——比較真實資料的 $T(y)$ 落在這 $S$ 個 $T(\tilde{y}^{(s)})$ 分布的哪裡。

定量指標是後驗預測 p 值（Bayesian p-value）：

$$ p_B = \Pr\big(T(\tilde{y}) \ge T(y) \mid \mathcal{D}\big) \approx \frac{1}{S}\sum_{s=1}^{S} \mathbb{1}\big[T(\tilde{y}^{(s)}) \ge T(y)\big] $$

$p_B$ 接近 $0.5$ 表示真實資料的這個面向落在模型生成的中央，模型抓得住；接近 $0$ 或 $1$ 則代表模型在這個面向系統性地偏離真實。

看一個例子：過度離散（overdispersion）的計數資料

你用 Poisson 模型擬合一份計數資料（如「每位學生每週發問次數」）。Poisson 有個鐵律：均值等於變異數。但真實的學習行為往往有少數同學狂發問、多數人沉默，變異數遠大於均值——這叫過度離散。

做 PPC：取 $T = \text{樣本變異數}$。你算出真實資料的變異數 $T(y) = 8.7$。然後從後驗預測抽 1000 組模擬資料，發現它們的變異數大多落在 $2.1$ 到 $4.5$ 之間，沒有任何一組超過 $8.7$。於是 $p_B \approx 0$——這是一記響亮的警告：Poisson 撐不住這份資料的離散程度。合理的修正是改用負二項分布（Negative Binomial）或加一層隨機效應。注意：傳統 AIC 或單看後驗均值，都不會這麼直白地把這個結構性錯配指給你看。

重點回顧

• HMC/NUTS 之所以高效，是因為它用後驗梯度 $\nabla_\theta U$ 導航，避免隨機漫步的高自相關；divergent transitions 代表後驗某區域被漏採，常需 non-centered 重參數化。
• 模型比較優先用預測導向的 PSIS-LOO 或 WAIC，而非對先驗極度敏感的貝氏因子；比較時要看差值的標準誤，差距小於約 2 SE 別硬分勝負。
• WAIC 把「擬合（lppd）」與「有效參數懲罰（$p_{\text{WAIC}}$）」分離，且後者用後驗對數似然的方差自動估計模型複雜度。
• 後驗預測檢驗讓模型自己生資料，用檢驗統計量 + 貝氏 p 值揪出結構性錯配（如 Poisson 對過度離散資料失效）。
• 「後驗算出來」只是起點：採樣診斷、模型比較、預測檢驗，三者缺一就稱不上完整的貝氏工作流程（Bayesian workflow）。

深入探討（研究所視角）

當資料量大到 MCMC 也跑不動：變分推論（Variational Inference, VI）。 MCMC 給的是漸近精確的樣本，但對大規模資料或複雜模型可能太慢。VI 改變策略——不採樣，而是把後驗推論轉成最佳化問題。它假設一族容易處理的分布 $q_\phi(\theta)$（如平均場高斯，mean-field Gaussian），透過調整 $\phi$ 去逼近真後驗，目標是最小化 KL 散度 $\mathrm{KL}(q_\phi \,\|\, p(\theta\mid\mathcal{D}))$。由於後者含難算的 evidence，等價地改成最大化證據下界（Evidence Lower Bound, ELBO）：

$$ \mathcal{L}(\phi) = \mathbb{E}_{q_\phi}\big[\log p(\mathcal{D}, \theta)\big] - \mathbb{E}_{q_\phi}\big[\log q_\phi(\theta)\big] = \log p(\mathcal{D}) - \mathrm{KL}\big(q_\phi \,\|\, p(\theta\mid\mathcal{D})\big) $$

因為 $\log p(\mathcal{D})$ 對 $\phi$ 是常數，最大化 ELBO 就等於最小化那個 KL。配合 reparameterization trick（把 $\theta = \mu + \sigma\odot\varepsilon$, $\varepsilon\sim N(0,I)$，讓梯度能穿過隨機節點）與隨機梯度上升，就得到 ADVI（Automatic Differentiation Variational Inference），能擴展到百萬筆資料。代價是：VI 用的「反向 KL」傾向低估後驗變異數、抓不住多峰與重尾，所以它快但通常比 MCMC 不保守。實務上常用 VI 快速探索、再用 MCMC 對最終模型精修。

幾何視角下的 funnel 與重參數化。 階層模型的漏斗幾何之所以毒害 HMC，根因在於後驗的局部曲率隨變異數參數劇烈變化——leapfrog 用的固定步長 $\epsilon$ 在漏斗頸部過大、在開口過小，沒有單一 $\epsilon$ 能同時適配。Riemannian Manifold HMC（黎曼流形 HMC）用 Fisher 資訊矩陣當位置相依的質量矩陣 $M(\theta)$ 來「拉直」幾何，是理論上的解；但計算成本高，所以實務多採非中心化參數化這個便宜的等價變換。

值得追問的方向：當模型本身可能設定錯誤（misspecified）時，標準貝氏更新會「過度自信」。這催生了 generalized / tempered Bayes（在似然上加溫度參數 $\eta$，$p(\theta\mid\mathcal{D})\propto p(\mathcal{D}\mid\theta)^{\eta}p(\theta)$）與 PAC-Bayes 這類對抗模型誤設的理論。把這些放回 Educational Omics 的脈絡——學習資料常是多模態、結構複雜且模型必然簡化的——這正是「校準的不確定性」比「精確的點估計」更值錢的場域。