21.5 后验计算、近似与诊断 · 概率论基础教程

后验分布定义了推断目标，数值算法负责计算这个目标。本节从后验积分出发，推导层次正态模型的 Gibbs 条件分布，说明重要性权重的可积条件，再讨论 Laplace 近似与 Bernstein–von Mises 定理。通过多峰分布和小样本伯努利模型，区分后验不确定性、积分误差与近似误差，并说明有效样本量和收敛诊断的适用范围。

先修与阅读位置

主要先修：14.3 蒙特卡洛误差与方差缩减、14.4 马尔可夫链蒙特卡洛、21.4 层次模型与经验贝叶斯。

技术部分另用：19.3 最大似然的渐近理论，适用段落见下述说明。

计算段接21.4，并按14.2–14.4补生成、重要性抽样与MCMC；后验渐近正态段另需19.3。模拟次数与样本量对应不同极限，不能把两段当作同一个收敛问题。查看全书依赖图与分支路线。

21.5.1 后验积分与三种误差

给定观测\(y\)，令\(\widetilde\pi(\theta)=L(\theta;y)\pi_0(\theta)\)为未归一化后验密度，相对于固定的支配测度积分。首先必须确认\(0<Z=\int\widetilde\pi(\theta)\,d\theta<\infty\)。后验均值、事件概率和预测概率都可以写成

\[\mu_g=E_\pi[g(\theta)]=\frac{\int g(\theta)\widetilde\pi(\theta)\,d\theta}{Z}.\]

(21.5-1)

定义21.5.1 计算目标与误差

后验不确定性指给定模型、先验和数据后，参数或预测量在后验下的分散程度；数值积分误差来自网格、截断或有限次模拟；近似误差来自用另一分布代替原后验。模型设定是否合适还需另外检查，不能由积分精度判断。

增加模拟次数可降低满足相应条件的蒙特卡洛误差，不会自动缩小后验方差，也不会修复错误的似然。低维积分可用数值求积，但须控制尾部截断并检验网格细化的影响；多维时可用14.2节的随机变量生成、14.3节的方差缩减与14.4节的 MCMC。以下不再重推 Metropolis–Hastings 接受概率。

21.5.2 层次正态模型的 Gibbs 更新

考虑固定的\(J\ge2\)组，每组\(n_j\ge2\)个观测，令\(N=\sum_jn_j\)。在21.4节的层次模型中，以组均值\(\theta_j\)为参数，使用精度\(\lambda=1/\sigma^2,\omega=1/\tau^2\)：

\[y_{jr}\mid\theta_j,\lambda\sim N(\theta_j,\lambda^{-1}),\qquad \theta_j\mid\mu,\omega\sim N(\mu,\omega^{-1}).\]

(21.5-2)

第一层给定各组均值与精度后独立，第二层给定\(\mu,\omega\)后组间独立。取相互独立的适定先验\(\mu\sim N(m_0,v_0)\)、\(\lambda\sim\operatorname{Gamma}(a_\sigma,b_\sigma)\)、\(\omega\sim\operatorname{Gamma}(a_\tau,b_\tau)\)，其中\(v_0>0\)，四个 Gamma 超参数均正；Gamma 使用形状／率，密度正比于\(x^{a-1}e^{-bx}\)。

组内不同观测的协方差是\(\tau^2\)，单个观测方差是\(\tau^2+\sigma^2\)，这里把\(\mu\)作为固定的边际模型参数，故重复观测可区分两层方差。若改写\(\theta_j=\mu+\alpha_j\)，本模型用\(E(\alpha_j\mid\mu,\omega)=0\)确定总体位置。不能再把所有\(\alpha_j\)当作不受约束的固定效应，否则平移\(\mu\)并反向平移它们保持似然不变。固定效应的和为零约束是另一种建模方案，不能悄悄加到本模型的独立正态随机效应上。

命题21.5.1 四类完整条件分布

以下“其余量”包括数据及全部未更新的未知量，\(\bar y_j\)为组内均值：

\[\theta_j\mid\text{其余量}\sim N\!\left(\frac{\lambda n_j\bar y_j+\omega\mu}{\lambda n_j+\omega},\frac1{\lambda n_j+\omega}\right).\]

(21.5-3)

\[\mu\mid\text{其余量}\sim N\!\left(\frac{m_0/v_0+\omega\sum_j\theta_j}{1/v_0+J\omega},\frac1{1/v_0+J\omega}\right).\]

(21.5-4)

\[\lambda\mid\text{其余量}\sim\operatorname{Gamma}\!\left(a_\sigma+\frac N2,b_\sigma+\frac12\sum_{j,r}(y_{jr}-\theta_j)^2\right).\]

(21.5-5)

\[\omega\mid\text{其余量}\sim\operatorname{Gamma}\!\left(a_\tau+\frac J2,b_\tau+\frac12\sum_j(\theta_j-\mu)^2\right).\]

(21.5-6)

完整推导与后验适定性证明。关于\(\theta_j\)的对数密度核为\(-[\lambda\sum_r(y_{jr}-\theta_j)^2+\omega(\theta_j-\mu)^2]/2\)。展开后平方项系数是\(\lambda n_j+\omega\)，一次项系数是\(\lambda n_j\bar y_j+\omega\mu\)，配方即得式(21.5-3)。对\(\mu\)合并先验平方项与全部组均值平方项，系数分别是\(1/v_0+J\omega\)、\(m_0/v_0+\omega\sum_j\theta_j\)，得到式(21.5-4)。

关于\(\lambda\)的核为\(\lambda^{a_\sigma-1+N/2}\exp\{-\lambda[b_\sigma+\sum_{j,r}(y_{jr}-\theta_j)^2/2]\}\)，关于\(\omega\)同理，得到两项 Gamma 分布。所有完整条件分布适定尚不能单独证明联合后验适定；这里另用似然上界\(L\le(\lambda/2\pi)^{N/2}\)。在归一化联合先验下积分，得\(Z\le(2\pi)^{-N/2}E_0\lambda^{N/2}<\infty\)，因为 Gamma 的所有正阶矩有限；似然在先验支持上正，故\(Z>0\)。证毕。

Gibbs 每步从完整条件分布随机抽取，后续更新使用最新状态。只代入条件均值是确定迭代，不是后验抽样。上述推导确定了目标与更新核，遍历性和误差极限定理仍需另行讨论。

例21.5.1 两组观测的条件分布核算

两组数据为\((0,2)\)与\((3,5)\)，取\(m_0=0,v_0=4\)，两项精度先验均为\(\operatorname{Gamma}(2,2)\)。在固定状态\((\theta_1,\theta_2,\mu,\lambda,\omega)=(1,3,2,1,1)\)核查各条件分布：

\[\theta_1\mid\cdots\sim N(4/3,1/3),\quad \theta_2\mid\cdots\sim N(10/3,1/3),\quad \mu\mid\cdots\sim N(16/9,4/9).\]

(21.5-7)

观测残差平方和是\(1+1+0+4=6\)，组均值相对总体位置的平方和是\(1+1=2\)，故精度条件分布为\(\operatorname{Gamma}(4,5)\)和\(\operatorname{Gamma}(3,3)\)。这些数字是同一状态下的条件分布检查，不是一轮顺序更新的抽样结果，更不是联合后验均值。

图21.5-1　式(21.5-3)的精确条件均值，两组条件方差均为\(1/3\)。图表示固定精度与总体位置时的收缩，不表示链的轨迹。

21.5.3 重要性权重与自归一化积分

从可计算密度\(q\)独立抽取\(\theta^{(1)},\ldots,\theta^{(M)}\)。须有\(q>0\)覆盖后验为正的区域，均按同一支配测度理解；否则漏掉的区域无法靠权重补回。定义

\[w_s=\frac{\widetilde\pi(\theta^{(s)})}{q(\theta^{(s)})},\qquad \widehat\mu_g=\frac{\sum_{s=1}^Mw_sg(\theta^{(s)})}{\sum_{s=1}^Mw_s}.\]

(21.5-8)

该比值仅在分母正时定义。若\(E_\pi|g|<\infty\)，大数定律给分母平均趋于\(Z>0\)、分子平均趋于\(Z\mu_g\)，因而估计一致；有限次模拟一般有偏。若进一步有\(E_q[w^2(1+g^2)]<\infty\)，独立样本 CLT 与比值展开给

\[\sqrt M(\widehat\mu_g-\mu_g)\Rightarrow N\!\left(0,\frac{E_q[w^2(g-\mu_g)^2]}{Z^2}\right).\]

(21.5-9)

推理是将差写成\([M^{-1}\sum_sw_s(g_s-\mu_g)]/[M^{-1}\sum_sw_s]\)，分子均值为零，分母收敛到\(Z\)。细节见14.3节及文献[1]。常用权重有效样本量\((\sum_sw_s)^2/\sum_sw_s^2\)衡量权重集中程度，但不保证支持覆盖，也不是每个\(g\)的误差上界。提议分布尾部过轻时，即使支持完全覆盖，二阶矩仍可能无限。

21.5.4 Laplace 近似：曲率与全局质量

令\(\ell(\theta)=\log L(\theta;y)+\log\pi_0(\theta)\)。若内部后验模态\(\widetilde\theta\)附近可作二阶展开，且负 Hessian 矩阵\(J=-\nabla^2\ell(\widetilde\theta)\)正定，则局部近似为

\[\ell(\widetilde\theta+u)=\ell(\widetilde\theta)-\tfrac12u^{\mathsf T}Ju+r(u),\qquad \pi(\theta\mid y)\approx N_d(\widetilde\theta,J^{-1}).\]

(21.5-10)

近似整个后验还要求：该模态邻域承载主要质量，在后验集中尺度上余项一致小，邻域外积分可忽略。仅验证 Hessian 正定不够。边界、重尾或多个重要模态可破坏这些条件；变换到无约束参数时还须在后验密度中计入 Jacobian。

例如精确密度\(\pi(\theta)=\{\phi(\theta-3)+\phi(\theta+3)\}/2\)，其中\(\phi\)为标准正态密度，均值为0、方差为10。右模态满足\(m=3\tanh(3m)\)，数值为\(m\approx2.99999991\)，负对数曲率\(J=1-9/\cosh^2(3m)\approx0.99999945\)。只保留右峰得到近似\(N(3,1)\)，其正半轴概率约0.99865，原分布的概率却恰为0.5。

图21.5-2　蓝线由精确混合公式计算；虚线用右模态及其数值曲率计算，点间连线。两条密度各自积分为1，图不是模拟产生的密度估计。

21.5.5 后验渐近正态的技术范围

Bernstein–von Mises 定理讨论数据量趋于无穷时的后验形状，与模拟次数趋于无穷不同。取正确设定、可识别的独立同分布模型，固定维数\(d\)，真实参数\(\theta_0\)为内部点，单个观测 Fisher 信息\(I_0\)正定。先验固定且适定，在真值附近有连续且正的密度。还须控制全局后验质量，不能仅由最大似然渐近正态推出结论。

命题21.5.2 在局部展开与集中前提下的后验正态极限

设最大似然估计\(\widehat\theta_n\)满足19.3节的正则渐近理论。令\(k_n=\log L_n\)、\(h=\sqrt n(\theta-\widehat\theta_n)\)，假设对每个固定\(K\)，

\[\sup_{\|h\|\le K}\left|k_n(\widehat\theta_n+h/\sqrt n)-k_n(\widehat\theta_n)+\tfrac12h^{\mathsf T}I_0h\right|\xrightarrow{P_{\theta_0}}0.\]

(21.5-11)

再假设缩放后验在该尺度上集中：对每个\(\varepsilon>0\)，\(\lim_{K\to\infty}\limsup_{n\to\infty}P_{\theta_0}\{\Pi_n(\|h\|>K\mid Y)>\varepsilon\}=0\)。那么

\[\left\|\Pi_n(h\in\cdot\mid Y)-N_d(0,I_0^{-1})\right\|_{\mathrm{TV}}\xrightarrow{P_{\theta_0}}0.\]

(21.5-12)

上述前提下的完整证明。总变差距离定义为对可测事件的概率差绝对值取上确界。固定球\(B_K\)，先验连续正、估计一致及式(21.5-11)，使缩放后的未归一化密度在球上一致趋于\(\pi_0(\theta_0)e^{-h^{\mathsf T}I_0h/2}\)，收敛按数据概率理解。球的体积有限且极限积分正，归一化后，后验条件于球内的密度与正态条件密度之差的积分趋于零。

将两项分布各自分成球内和球外，总变差距离不超过这项条件距离加两项球外概率。先选大\(K\)使正态尾部小；后验集中前提使后验尾部也以高数据概率小；再令\(n\)增大消去球内条件距离。对任意误差容限重复此步骤，得到式(21.5-12)。证毕。这里证明的是所列技术前提的蕴含；在一般模型中验证后验集中还涉及检验或尾部控制，不能用一次 Taylor 展开代替。文献[2]给出后验正态理论与具体模型。

对于正则标量参数，或真值处梯度非零的光滑标量函数，结合估计量的非退化频率渐近正态，通常的等尾可信区间可有渐近频率覆盖。该结论不是有限样本精确保证，也不覆盖任意形状的可信集。边界、先验不支持真值、模型错设或维数随样本增长，均须另外分析。

例21.5.2 伯努利后验：精确积分与小样本近似

独立伯努利样本8次成功6次，先验\(p\sim\operatorname{Beta}(2,2)\)，于是精确后验是\(\operatorname{Beta}(8,4)\)。其均值、方差和内部模态分别为\(2/3\)、\(2/117\approx0.017094\)、\(7/10\)。对数密度核\(7\log p+3\log(1-p)\)的负二阶导数在模态处为

\[J=\frac7{0.7^2}+\frac3{0.3^2}=\frac{1000}{21},\qquad \pi(p\mid Y)\approx N(0.7,0.021).\]

(21.5-13)

正态近似给区间\([0,1]\)外约0.01922的概率，精确后验在那里概率为零。精确等尾95%可信区间由 Beta 分布函数反演得\([0.390257,0.890737]\)，端点是数值求积与求根结果。若固定先验并令\(n\to\infty\)、\(p_0\in(0,1)\)，后验的一阶正态近似才以样本比例为中心、方差约\(\widehat p(1-\widehat p)/n\)；不能把这项极限当作8次试验的精度证明。

21.5.6 MCMC 误差与诊断的边界

14.4节说明不变分布与遍历性。后验适定和更新保持不变性并不单独保证链快速混合。遍历大数定律需合适的不可约性、Harris 遍历性等条件以及\(g\)可积；常规误差条还需要更强的 Markov 链 CLT。设链已平稳，\(g\)平方可积，自协方差绝对可和，且该 CLT 成立，定义

\[\gamma_k=\operatorname{Cov}_\pi(g(\theta^{(0)}),g(\theta^{(k)})),\quad \eta_g^2=\gamma_0+2\sum_{k=1}^{\infty}\gamma_k,\quad \operatorname{MCSE}\approx\sqrt{\eta_g^2/M}.\]

(21.5-14)

若\(\eta_g^2>0\)，相应有效样本量为\(\operatorname{ESS}_g=M\gamma_0/\eta_g^2\)，随所计算的函数改变。后验标准差是\(\sqrt{\gamma_0}\)，MCSE 则衡量后验期望估计的数值波动。自相关估计和长程方差估计本身也有误差。

应比较分散初值的多条链，检查轨迹、秩归一化的分裂\(\widehat R\)、主体与尾部 ESS，并针对所报告均值或分位数检查精度[3]。这些诊断发现部分问题，没有证明后验已被完整探索：图21.5-2中，所有链都停在右峰，也可能得到接近1的\(\widehat R\)。预热长度不能普遍保证去除初值影响；若预热时调节更新核，可在正式采样时固定该核，继续自适应则须满足相应自适应遍历条件。

21.5.7 习题与解答

习题21.5-1　更换总体均值的先验方差

保持例21.5.1的状态和全部数据，仅把\(v_0\)改为1。求\(\mu\)的完整条件分布，并解释此时是否需要对独立组随机效应施加和为零约束。

先验精度从\(1/4\)变成1，两组条件精度仍合计2，故总精度为3。一次项为\(0+1+3=4\)，配方得到\(\mu\mid\cdots\sim N(4/3,1/3)\)，条件均值更靠近先验均值0。独立正态随机效应已经通过条件分布的零均值确定总体位置；实现的组偏差之和不必等于零。再施加精确和为零会改变联合分布，必须重推条件分布，不能保留原更新公式。

习题21.5-2　支持覆盖却没有有限方差

后验为\(N(0,1)\)，重要性提议为\(N(0,1/4)\)，目标\(g(\theta)=\theta\)。求权重，判断自归一化均值是否一致，以及式(21.5-9)是否适用。

两项密度之比为\(w(\theta)=\tfrac12e^{3\theta^2/2}\)，支持均是实数轴。\(E_qw=1\)且\(E_q|w\theta|=E_\pi|\theta|<\infty\)，故分子均值趋于0、分母均值趋于1，比值一致趋于0。但\(w^2q\)正比于\(e^{\theta^2}\)，因此\(E_qw^2\)和\(E_q(w^2\theta^2)\)均无限。有限方差 CLT 的条件失败，不能使用该渐近方差或通常的正态误差条；有限批次样本恰巧集中在原点附近也不能证明尾部误差小。

习题21.5-3　相关抽样的 MCSE

平稳链满足\(\theta_{s+1}=0.9\theta_s+\sqrt{1-0.9^2}\,\varepsilon_{s+1}\)，创新独立标准正态且独立于过去。计算\(g(\theta)=\theta\)、\(M=1900\)时的长程方差、ESS 与渐近 MCSE。

平稳分布为\(N(0,1)\)。迭代递推式，创新与\(\theta_0\)独立，得\(\gamma_k=0.9^k\)，几何级数给\(\eta_g^2=1+2(0.9/0.1)=19\)。该高斯 AR(1)链满足 CLT，故\(\operatorname{ESS}_g=1900/19=100\)、MCSE 近似为\(\sqrt{19/1900}=0.1\)；后验标准差仍是1。有限次均值的精确方差为\([M+2\sum_{k=1}^{M-1}(M-k)0.9^k]/M^2\)，所以0.1是长程渐近误差尺度，不是有限次精确标准差。

习题21.5-4　局部诊断不能验证先验支持

独立伯努利数据的真值为\(p_0=0.8\)，先验只支持\((0,0.5)\)。若多条后验链的\(\widehat R\)接近1且 ESS 很大，能否援引式(21.5-12)，或断言95%可信区间覆盖真值？

对任意数据量，后验在\(p>0.5\)的概率仍为零，因为先验密度在那里为零。后验不可能在真值附近集中，先验在真值附近正的前提已经失败，局部后验集中前提也失败。取通常位于后验支持内的等尾可信区间，其端点都不超过0.5，因而频率覆盖\(p_0=0.8\)的概率为零。诊断良好最多支持链已有效探索这个指定后验，不能检验模型支持是否包括真值。

关键术语

后验积分：在归一化后验分布下计算函数的期望。
完整条件分布：给定全部其余未知量及数据时，一个未知量的条件分布。
Gibbs 抽样：轮流从完整条件分布抽取状态的 MCMC 方法。
重要性权重：目标未归一化密度与提议密度之比。
自归一化：以样本权重总和消去未知归一化常数。
Laplace 近似：利用模态附近对数密度的二次展开构造局部正态近似。
总变差距离：两个概率分布对同一可测事件的概率差绝对值的上确界。
有效样本量：按特定方差或权重标准折算的独立抽样数量，须注明定义。
MCSE：蒙特卡洛标准误，描述数值估计量因有限模拟产生的波动。
分裂 Rhat：将各链分段后比较链间与链内变异的诊断量，秩归一化版本也检查尺度等差异。

参考文献

Owen, A. B. Monte Carlo Theory, Methods and Examples，第9章，尤其第9.2、9.3、9.7节。作者官方版本：自归一化重要性抽样、二阶矩条件、权重诊断与模态曲率。本节的支持与正态提议反例独立计算。
Kempthorne, P. (2016). MIT 18.655, Lecture 18: Asymptotics III: Bayes Inference and Large-Sample Tests，后验渐近正态及伯努利例。大学官方讲义。命题21.5.2明确把一致局部展开和后验集中作为前提，证明它们如何推出总变差正态极限，不证明任意模型自动满足这些前提。
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. and Bürkner, P.-C. (2021). Rank-Normalization, Folding, and Localization: An Improved Rhat for Assessing Convergence of MCMC, Bayesian Analysis 16(2), 667–718。作者所在大学版本：秩归一化分裂 Rhat、主体及尾部 ESS，并讨论诊断不完备的限制。