校准贝叶斯推理

Calibrating Bayesian Inference

校准贝叶斯推理

摘要
尽管贝叶斯统计因其直观的不确定性量化和灵活的决策能力而在心理学研究中广受欢迎，但其在有限样本中的表现可能不可靠。本文揭示了一个关键弱点：当研究者所选的先验分布与真实参数生成过程不匹配时，贝叶斯推断在长期运行中可能产生误导。鉴于真实过程在实践中通常未知，我们提出一种更安全的替代方案：对贝叶斯可信区域进行校准，以实现频率学派的有效性。这一标准更强，能保证无论底层参数生成机制如何，贝叶斯推断均有效。为在实践中解决校准问题，我们提出了一种新颖的随机逼近算法。我们开展并报告了一项蒙特卡洛实验，结果表明，在某些参数生成情形下，未经校准的贝叶斯推断可能过于宽松（liberal），而我们的校准方案始终能保持有效性。

关键词：贝叶斯推断，频率学派推断，统计有效性，可信区域，随机逼近，黎曼优化

引言

近几十年来，使用贝叶斯方法的心理学出版物显著增加（例如，Kruschke, 2021；van de Schoot 等, 2021；van de Schoot, Winter, Ryan, Zondervan-Zwijnenburg & Depaoli, 2017；Volpe 等, 2025）。贝叶斯推断通过后验概率度量提供直观的不确定性量化。借助从后验分布中抽取的（近似）随机样本，可方便地以无需解析推导的方式对模型参数进行推断、预测未来数据并评估模型拟合（Gelman, Carlin, Stern & Rubin, 2013）。

现成的贝叶斯分析软件不仅包括通用的马尔可夫链蒙特卡洛（MCMC）抽样器，如 JAGS（Plummer, 2017）和 Stan（Stan Development Team, 2024），还包括专为特定建模框架设计的程序，如 Mplus（Muthén & Muthén, 1998–2024）和 blavaan（Merkle & Rosseel, 2018）。

贝叶斯推断的哲学与统计基础在于在不确定性下做出一致的决策（DeGroot, 1970；Savage, 1954）：用恰当的概率测度编码先验知识，并在观测新数据后通过贝叶斯公式更新知识。然而，贝叶斯方法在心理学中的应用大多是工具性的，而非哲学性的。将先验明确联系到研究者信念的正当性说明仍然罕见，而依赖默认或共轭先验的做法却十分普遍。通常有三种理由被用来为此辩护。

第一，某些默认先验——尤其是“弱信息”和“客观”先验（例如 Berger, Bernardo & Sun, 2024, 2015；Datta & Mukerjee, 2004；Gelman, Jakulin, Pittau & Su, 2008）——在现有文献中已被证明具有良好的理论性质或强健的实证表现。第二，在适当的正则条件下，随着样本量增大，先验的影响会减弱，所得的后验推断在大样本下通常与频率学派结果相似（例如 Bernstein-von Mises 定理；van der Vaart, 1998，第10章）。第三，通常建议进行敏感性分析，以确保统计结论对不同先验选择的稳健性（例如 Depaoli, 2022；Depaoli, Winter & Visser, 2020；Van Erp, Mulder & Oberski, 2018）。

然而，当应用于现实世界的心理学研究时，上述辩护的合理性常常值得怀疑。首先，“客观性”和“无信息性”在定义默认先验时并未建立在单一、统一的框架之上（例如 Kass & Wasserman, 1996）。此外，哪种默认先验表现最佳往往取决于具体模型（Yang & Berger, 1998）。因此，寻找一个在所有情况下都表现优异的先验很可能是一个难以实现的目标。

第二，由于研究焦点或实际考量，实质性研究者可能不得不处理小样本。例如，由多重社会身份定义的交叉性亚群体通常过于狭窄，难以积累足够数据（例如 Cole, 2009）。基于大样本理论的统计程序在小样本应用中可能数值不稳定，甚至产生误导性推断（例如 van de Schoot & Miočević, 2020）。

第三，先验敏感性分析可能无法得出明确结论。几乎总能找到一种病态的先验分布（例如，其质量集中在远离原始贝叶斯解的区域），从而推翻原有结论。此外，贝叶斯计算可能计算成本过高，难以重复大量次数。因此，先验敏感性分析通常局限于一组有限且任意选择的先验，对先验设定的诊断价值甚微。

为应对现实中“实用主义贝叶斯”（pragmatic Bayes）的普遍性，贝叶斯方法的方法论研究越来越关注在数据和/或参数重复抽样下的表现。事实上，任何在长期运行中系统性地产生错误结论的推断程序都应被摒弃。为此，过去几十年开展了大量蒙特卡洛（MC）实验，其中贝叶斯推断程序（即假设检验和区间估计）在各种数据与参数生成机制及设计因素（如样本量、协变量数量等）下被评估（例如 Finch & French, 2019；McNeish, 2016, 2017a, 2017b；Preacher & MacCallum, 2002；Smid, McNeish, Miočević & van de Schoot, 2020）。然而，MC 研究的一个主要局限在于其结论依赖于特定模型和设计，难以推广到未明确测试的情境之外。因此，目前使用贝叶斯分析的心理学研究所得结论的可信度尚未完全确立。

本文有两个主要目标：一是教学性的，二是方法论的。借鉴统计决策理论（Berger, 1985）和推断模型（Inferential Models, IMs；C. Liu & Martin, 2024；Martin & Liu, 2015）的相关成果，我们强调贝叶斯有效性（Bayesian validity）这一核心概念：即关于参数的不合理陈述不太可能频繁地具有较大的后验概率。这一原则为贝叶斯程序的长期频率表现提供了正当性（例如 Martin, 2022a, 2022b, 2022c；精确的定义见“贝叶斯推断与统计有效性”一节）。在缺乏先验知识、仅将贝叶斯方法作为工具使用的场景下，我们指出应追求更强的频率学派有效性（frequentist validity），因为它必然蕴含对任意先验设定下的贝叶斯有效性。

在方法论层面，我们提出一种计算程序，用于校准基于后验的推断以确保频率学派有效性。该方法结合了无梯度随机逼近（gradient-free stochastic approximation, SA）与流形优化（manifold optimization）（例如 Absil, Mahony & Sepulchre, 2008；Spall, 1992）。我们在一项蒙特卡洛实验中将该方法应用于线性-正态因子分析模型（例如 Jöreskog, 1969），结果表明：若未经适当校准，贝叶斯推断可能是无效的。

本文其余部分结构如下：我们首先回顾统计决策理论、推断模型（IMs）和统计有效性的理论基础。具体而言，我们阐明一个关键事实：频率学派有效性可确保对任意先验选择下的贝叶斯有效性。接着，我们介绍一种实用的计算算法，用于校准贝叶斯推断以实现频率学派有效性。我们开展了一项概念验证型蒙特卡洛（MC）实验，将基于随机逼近（SA）的校准推断与基于MCMC抽样的原始推断进行对比。最后，文章通过讨论本研究的启示、局限性及未来研究方向予以总结。

贝叶斯推理与统计有效性

贝叶斯模型和嵌套置信区域

后验可能性

嵌套置信区域不仅仅是模型参数的一组估计量。它们在贝叶斯推理中的基础作用可以通过可能性理论正式建立（Dubois, 2006; Dubois & Prade, 1988; Zadeh, 1978）。设

后验可能性轮廓可以通过在所有 α 级上拼接嵌套置信区域族来直观地描绘；图1中可以找到图形说明。更深入的可能性理论及其在统计推理中的应用可以在例如Denoeux和Li (2018)，C. Liu和Martin (2024)，以及Martin (2025b)中找到。

后验可能性是贝叶斯框架内推断的基石。可能性等值线函数 ᵧ 为所有嵌套可信区域提供了一个简洁的总结：对于任意给定的 α ∈ [0,1]，可通过取该等值线的上 α-水平集直接获得一个集合估计量。此外，在贝叶斯检验中，零假设的可信度由其后验可能性保守地量化。例如，若简单假设 H = {θ₀} 的可能性 Π̄ᵧ{H} = ᵧ(θ₀) 小于预设的（较小的）α 水平，则该假设被拒绝。

统计有效性

为了使任何统计程序在实践中可靠，最好能确立该程序在各种情境下均能产生可靠的结果。特别是在对模型参数进行推断时，重要的是不应在长期运行中反复将低可能性值赋予频繁发生的事件。越来越多的贝叶斯和频率学派统计学家，尽管对概率代表什么持有不同观点，但仍支持评估推断程序在重复抽样下表现的重要性（例如 Grünwald, 2018；Martin, 2022a, 2022b, 2022c）。接下来，我们将回顾贝叶斯有效性和频率学派有效性的概念。我们强调，基于后验可能性测度的贝叶斯推断，在参数与数据模型均被正确设定的假设下，满足贝叶斯有效性。同时，具有频率学派有效性的程序，在数据模型被正确设定的前提下，对所有先验而言自动满足贝叶斯意义下的有效性。

贝叶斯有效性

令 h : Y × Q → ℝ 为任意 P_{Y,Θ}-可积函数。根据富比尼定理（Fubini’s）

定理（Billingsley, 2012，定理 18.3）指出，我们可以将函数 h 关于数据和参数的联合期望写成如下迭代期望的形式：

校准置信区域

在本节中，我们提出了一种通用的计算策略，通过阈值化观测检验统计量来校准后验可能性轮廓并确保频率论有效性（8）。作为输入，我们定义了一个通过阈值化观测检验统计量来定义的嵌套置信区域族。然后，我们使用无梯度SA算法基于检验统计量的生存函数来校准置信区域。值得注意的是，我们解决的优化程序本质上等同于为IM可能性轮廓（9）找到一个“变分近似”（Cella & Martin, 2024; Martin, 2025a）。我们提案的一个独特贡献是同时扰动SA（Spall, 1992, 2000, 2009）和流形优化的新组合，这增强了校准算法的可扩展性，并促进了其在现实规模模型中的应用。

检验统计量和嵌套置信区域

具体在当前工作中，我们从两种类型的检验统计量构建嵌套置信区域：Wald统计量和后验密度比（PDR）统计量，分别生成椭圆形和HPD置信区域。

Wald 统计量与椭圆区域

与Wald统计量相关的嵌套置信区域族可以表示为：

后验密度比统计量和最高后验密度区域

或者，我们可以根据PDR统计量定义置信区域：

（13）是标准大样本理论中最大似然估计的似然比统计量的推广。还要注意，PDR统计量（13）是对Martin（2022b）在更一般的部分先验背景下定义的“相对合理性排序”的对数变换。类似于（12），可以通过收集具有足够低PDR统计量值的参数值来构建置信区域：

优化问题中，其可行区域形成欧几里得空间的可微分子流形，可以通过黎曼梯度算法有效解决（Absil et al. 2008；另见Y. Liu 2020, 2021，关于黎曼梯度算法及其在心理测量问题中的应用的可访问介绍）。特别地，黎曼梯度上升算法将欧几里得空间中的常规梯度上升算法推广，具有两个显著区别。首先，最陡上升方向只能在当前迭代处局部定义，通过黎曼梯度获得，这是通过将目标函数的周围梯度投影到子流形的切空间来实现的。对于我们的问题（15），

在应用标准黎曼梯度算法求解[15]时，另一个挑战来自于对p值函数的精确梯度∇θπy(θ)的评估。这是因为p值函数是一个积分，其定义域通过检验统计量依赖于θ。尽管有时可以通过广义莱布尼茨法则（也称为雷诺兹输运定理；参见例如Flanders, 1973; Reddiger & Poirier, 2023）获得解析导数，但一般而言，计算这些导数是困难的，因为检验统计量本身的计算可能涉及优化。一种简单的尝试是用有限差分（FD）估计来近似梯度的第r个元素（r = 1, ..., q），记为∇θrπy(θ)：

对于某个小的扰动 c > 0，其中 er 是一个在第 r 个元素处为 1、其余位置为 0 的单位向量。由于方括号内的项在(20)中是其期望值的无偏估计量，我们或许可以期望通过黎曼随机逼近（SA）算法（例如，Bonnabel, 2013; Tripuraneni, Flammarion, Bach, & Jordan, 2018）找到(15)的一个近似解。然而，这种近似存在两个问题：(a) 如果 c 固定，有限差分近似的偏差不会消失；(b) 随着参数数量增加，需要对检验统计量进行更多次评估，从而导致“维度灾难”。

为解决这些问题，我们提出一种SPRSA算法，它结合了针对所有 q 个参数的同时扰动有限差分法（解决了问题(a)），以及沿黎曼梯度迭代 k = 0, 1, ... 进行衰减的有限差分步长序列 {ck}（解决了问题(b)）。设 Y = g(U, θ) 为一个数据生成算法，其中随机分量 U ~ PU，且分布 PU 完全已知。在第 k 次迭代时，∇θπy(θ(k)) 的同时扰动有限差分估计量定义为

使用SPRSA算法实现实际效率，需要仔细地、逐案调整多个方面：学习率序列{ak}、有限差分率序列{ck}以及总迭代次数K。在我们的数值研究中，SPRSA算法的调参细节将在“模拟研究”部分的“抽样与调参细节”中提供。

蒙特卡洛实验

数据生成

协方差结构模型已被广泛用于解释观测响应变量之间的相关模式（Bollen, 1989; Jöreskog, 1970）。令 Z = (Z₁, ..., Zₙ)ᵀ ∈ ℝⁿˣᵐ 为独立同分布（i.i.d.）的样本数据，其中每个 m 维响应向量 Zᵢ ~ N(μ, Σ(θ))，其均值向量 μ ∈ ℝᵐ，协方差矩阵 Σ(θ) ∈ ℝ⁺ᵐˣᵐ 由参数 θ ∈ ℝq 参数化。在我们的数值示例中，我们考虑一个一维公共因子模型（Jöreskog, 1969）：

在(24)中，λ = (λ₁, ..., λₘ)ᵀ ∈ ℝᵐ 是因子载荷向量，ψ ≥ 0 表示公共因子方差，u = (u₁, ..., uₘ)ᵀ ∈ ℝ₊ᵐ 收集了唯一因子方差参数。为了识别该模型，要求公共因子 ηᵢ ∈ ℝ 与唯一因子 εᵢ ∈ ℝᵐ 相互独立，并且第一个因子载荷 λ₁ 固定为1。以下协方差矩阵由公共因子模型推导得出：

为了避免对参数空间施加边界限制，我们对所有方差参数进行了对数变换：

令 ζ = log ψ / 2 表示公共因子方差的对数标准差（SD），并且

ωⱼ = log uⱼ / 2，j = 1, ..., m，表示第 j 个唯一因子方差的对数标准差。我们将这些无界参数收集在参数向量 θ = (ζ, λ₂, ..., λₘ, ω₁, ..., ωₘ)ᵀ 中，其维度为 q = 2m。由于 Zᵢ, i = 1, ..., n 的独立同分布正态性，样本交叉乘积矩阵 Y = Σᵢ₌₁ⁿ (Zᵢ − Z̄)(Zᵢ − Z̄)ᵀ（其中 Z̄ 表示样本均值向量）服从 Wish(Σ(θ), n−1) 分布，即尺度矩阵为 Σ(θ)、自由度为 n−1 的威沙特分布。由于交叉乘积矩阵 Y 是协方差结构 Σ(θ) 的充分统计量，我们在整个模拟研究中将 Y 视为数据。Y 的直接数据生成算法为 Y = g(U, θ) = Σ(θ)¹/² U Σ(θ)¹/²，其中 U ~ Wish(Iₘₓₘ, n−1)。

模拟条件由两个交叉因素决定：响应变量的数量（m = 5 和 15），以及三种参数生成情景（情景1–3）。样本量固定为 n = 100。在情景1中，响应变量的共性在各次重复中随机从 U[.2, .8] 中抽取，涵盖从低到高的共性水平（MacCallum, Widaman, Zhang, & Hong, 1999）。公共因子方差被设定为第一个响应变量的共性值。唯一方差被确定为使得所有响应变量具有单位方差。由此产生的 Θ 的分布作为 P*Θ。在情景2中，所有响应变量具有固定的低共性（.3）。在情景3中，所有响应变量具有固定的高共性（.7）。公共因子和唯一因子方差的确定方式与情景1类似。情景2和3中的参数生成分布 PΘ 是狄拉克测度（即，点质量集中在真实参数上）。我们使用 MATLAB 版本 23.2（MathWorks, 2023）生成数据；每种条件下运行了 512 次重复实验。

抽样与调参细节
我们采用了 Asparouhov 与 Muthén（2010）所推荐的先验设定：对因子载荷使用非正常均匀先验（improper uniform priors），对公共因子方差和唯一因子方差使用逆伽马先验 IG(1, 2)。MCMC 抽样通过 JAGS（Plummer, 2017）完成。由于 JAGS 无法处理非正常先验，我们改用一个弥散的正态先验 N(0, 10¹⁰) 来代替因子载荷的先验，该先验与非正常均匀先验几乎无法区分。在每次重复实验中，我们并行运行 5 条链。每条链的自适应迭代次数、预烧期（burn-in）迭代次数和保留迭代次数分别为 1000、10000 和 10000。利用这 5 条链的保留迭代样本，以 10 的间隔进行抽稀（thinning），共获得 5000 个 θ 的蒙特卡洛样本。在每次重复中，记录 θ 每个分量的潜在尺度缩减因子（PSRF）和有效样本量（ESS）。若所有参数均满足 PSRF ≤ 1.1 且 ESS ≥ 100，则认为该次重复实验已收敛。

为评估后验可能性并进行校准，我们计算了Wald检验统计量[11]和PDR统计量[13]。我们使用MATLAB内置的信赖域算法（通过fminunc）对关于θ的对数后验进行最大化；我们提供了后验的解析梯度和期望Hessian矩阵以加速优化过程。相同的期望Hessian矩阵也用于定义Wald检验统计量，以近似MAP估计量的协方差矩阵。

为调优SPRSA算法，我们进行了初步模拟实验。具体而言，学习率序列由 aₖ = αk⁻ᵝ 确定，其中 α = 0.1，β = 0.65；有限差分率序列设定为 cₖ = γk⁻ᵟ，其中 γ = 0.05，δ = 0.149。可以直观验证，这两个速率序列满足条件[22]。算法的迭代次数设定为 K = 50000。最后一次迭代的有限差分扰动值为 c₅₀₀₀₀ = 0.05 × 50000⁻⁰·¹⁴⁹ ≈ 0.01。校准后的α水平的平均估计值使用公式[23]计算得出。我们在MATLAB中实现了SPRSA算法；若本文被接受发表，我们将在补充材料中提供源代码。

评估标准

我们评估了基于后验推断的贝叶斯有效性，无论是否经过校准。在每次重复实验中，我们使用观测数据 y 和数据生成参数 θ ~ P*Θ 来计算原始后验等高线 ϖy(θ) 和校准后的后验等高线 ̃ϖy(θ)。随后，我们在每种模拟条件下，针对原始后验等高线值和校准后后验等高线值，在所有重复实验中分别构建了经验分布函数（EDFs）。若EDF曲线位于对角线下方，则表明推断是保守的，因而也是有效的；而若曲线位于对角线上方，则表明有效性要求未被满足。为考虑蒙特卡洛误差（MC error），与对角线的比较将参照其95%正态近似蒙特卡洛置信带（即，

所有 θ 的联合推断结果总结于图2中。当响应变量数量较小时（即，m = 5），基于椭圆和HPD可信区域的贝叶斯推断往往无效；一个例外出现在情景3中的椭圆区域。当共性固定且较低时（情景2），椭圆区域表现得最为宽松；而当共性均匀生成时（情景1），HPD区域表现得最为宽松。相比之下，随着响应变量数量增加至 m = 15，有效性通常得以恢复。但在情景3中，当共性固定且较高时，存在一个例外：基于HPD区域的推断仍保持不可接受的宽松性。

同时，使用所提算法进行校准后的后验可能性在所有模拟条件下均实现了有效性。然而，有效性保证会引入保守性，其程度取决于若干因素。特别是，校准后的贝叶斯推断在响应变量较少（即，m = 5）和共性较低（即，情景2）时变得更加保守。校准后，HPD区域通常比椭圆区域更不保守。值得注意的是，当 m = 15 时，基于HPD区域的校准可能性等高线几乎在所有三个参数生成情景下都达到了一致性。

对这一观察的一个可能解释是，与椭圆区域相比，HPD区域能更好地近似 πy 的上层水平集的形状。

示例

为结束本节，我们说明校准结果在实践中如何呈现。当分析单个数据集并选定检验统计量后，我们可以在一组预设的阈值 ξ₁, ..., ξQ 上重复校准程序，从而得到一组校准后的 α 水平，α*(ξ₁), ..., α*(ξQ)。为了突出校准的效果，我们建议选择 {ξq} 作为在一组名义 α 水平下的检验统计量的 (1−α) 后验分位数。然后可将校准后的 α 水平相对于后验分位数的名义 α 水平作图，创建一种类似于百分位-百分位图的可视化图形。

我们以在 m = 5 且参数按情景1生成条件下生成的最后一组数据为例。基于保留的MCMC抽样，Wald统计量和PDR统计量的估计后验密度显示在图3左侧面板中。对于每个统计量，令 Q = 19，ξ₁, ..., ξ₁₉ 为后验分布的 .95, .9, ..., .05 分位数。建议的图形展示呈现在图3右侧面板中。在此数据集中，校准推断对于两种类型的可信区域而言，均比原始后验推断更为保守。更具体地说，对于HPD区域，α 水平调整幅度比椭圆区域更大。

贝叶斯统计在心理学家中广受欢迎，因其能提供直观的不确定性量化、适用于多种建模场景，并且在小样本情况下有时表现优异（例如，Depaoli, 2021；Muthén & Asparouhov, 2012；van de Schoot 等, 2021）。然而，借助推断模型（IM）中关于贝叶斯有效性的关键概念，我们表明：当用于推断的先验与真实参数生成机制（即真实的参数生成先验）不匹配时，贝叶斯方法在重复抽样下可能是不可靠的。由于在实际应用中，数据分析者通常无法获知这一真实生成机制，因此我们提出一种更安全的替代方案：对基于后验的推断进行校准，以实现频率学派的有效性（frequentist validity）——这是一种更强的要求，能够保证在任意参数生成先验下都满足贝叶斯有效性。

为解决校准问题，我们开发了一种SPRSA算法，该算法将流形优化与无梯度随机逼近（SA）相结合。随后，我们报告了一项针对简单单因子模型的蒙特卡洛实验。结果表明，采用一种广泛使用的先验设定的标准贝叶斯推断，其有效性会因可信区域类型、响应变量数量以及真实参数生成机制的不同而失效。相比之下，经过校准的贝叶斯推断在所有模拟条件下均实现了有效性。此外，我们还证明了SPRSA算法可扩展至心理学应用中常见的现实问题规模。文中也提供了校准结果的建议可视化图形展示方式。

本研究存在若干局限，有待未来研究加以解决。首先，我们的模拟仅限于一个简单的单因子模型。鉴于贝叶斯方法的广泛应用，有必要开展更全面的蒙特卡洛实验，以评估常用先验在多大程度上会导致无效推断，并凸显校准的普遍必要性。其次，Wald 统计量和 PDR 统计量的使用要求在 SPRSA 算法的每次迭代中求解两次最大后验（MAP）估计，对于复杂模型而言，这可能带来较大的计算负担。未来的研究可探索替代的有限差分（FD）梯度估计方法或检验统计量，以进一步减轻计算开销。最后，我们的方法假设 p 值函数是可微的，因此无法直接应用于离散数据问题。一种有前景的解决方案是对检验统计量加入随机扰动，从而强制其分布变为连续的。

https://www.researchgate.net/publication/397188720_Calibrating_Bayesian_Inference