附录 K：贝叶斯工作流程¶

【摘要】贝叶斯数据分析方法为使用概率论处理观测、模型参数和模型结构中的所有不确定性提供了一种强大方法。概率编程语言使指定和拟合贝叶斯模型变得更加容易，但这仍然给我们留下了许多关于构建、评估和使用这些模型的选择，以及许多计算方面的挑战。使用贝叶斯推断解决实际问题不仅需要统计技能、学科知识和编程，还需要在数据分析过程中做出决策。所有这些方面都可以理解为应用贝叶斯统计的复杂工作流程的组成部分。除了推断以外，工作流程还包括迭代模型构建、模型检查、计算问题验证及故障排除、模型理解和模型比较等。本文在几个示例背景下，回顾了工作流程的各个方面。需要提醒读者：在实际工作中，应当为任何给定的问题拟合多个模型，即使最终可能只有非常有限的子集能够与结论相关。

【原文】Gelman, Andrew and Vehtari, Aki and Simpson, Daniel and Margossian, Charles C and Carpenter, Bob and Yao, Yuling and Kennedy, Lauren and Gabry, Jonah and Burkner, Paul-Christian and Modrak, Martin; Bayesian workflow; arXiv:2011.01808}; 2020.

1 概述¶

1.1. 从贝叶斯推断到贝叶斯工作流程¶

如果数理统计是『关于如何应用统计学』的理论，那么任何对贝叶斯方法的认真讨论都需要明确其在实践中的使用方式。特别是，我们需要明确地将贝叶斯推断与贝叶斯数据分析的概念分开，特别是将其与完整的贝叶斯工作流程分开。

注解：作者强调贝叶斯推断、贝叶斯数据分析及相应的贝叶斯工作流程的概念需要有非常明确的界定。

贝叶斯推断：只是条件概率（或条件概率密度） \(p(\theta \mid y) \propto p(\theta) p(y \mid \theta)\) 的形式化和计算，可视为一种计算工具。
贝叶斯工作流程：包括模型构建、统计推断和模型检查/改进 三个步骤，以及不同模型之间的比较 （比较不仅是为了模型选择或平均，更是为了深入理解这些模型）。例如：为何某些模型无法预测数据的某些方面？或者一些重要参数的不确定性估计为什么会因模型而不同？即使我们有一个非常喜欢的模型，将其与更简单和更复杂模型推断进行比较，对于理解模型也是很有用的。

图 1 提供了一个贝叶斯工作流程的概要。扩展的贝叶斯工作流程应当还包括（在收集数据和测量数据前的）试验设计以及推断完成后的决策制定，本文重点仅限于对数据进行建模。

图 1：在贝叶斯工作流程中考虑的主要步骤。该图旨在展示某个贝叶斯分析任务可能经历的步骤和路径，当然应该理解任何特定分析都不太可能涉及所有步骤。研究贝叶斯工作流程的目的之一是掌握图中这些步骤和想法是如何组合在一起的，以便更系统地使用它们。

在典型贝叶斯工作流程中，我们会拟合出一组模型，其中某些模型可能比较糟糕（原因可能包括：数据拟合不佳、缺乏与实际理论或实际目标的联系、先验太弱/太强/不合适、编程错误等）；另外一些模型有用但可能存在缺陷（例如：无依据地随意选择混淆因子或变量进行回归、只捕获了部分而非全部功能关系的参数形式等）；另外一些结果则可能有价值，值得采纳。因此，在实践中必须认识到：错误模型、有缺陷模型都会是拟合有用模型无法避免的前驱步骤。理解这一点有助于帮助我们改变设置和使用统计学方法的方式。

1.2 为什么需要贝叶斯工作流程？¶

需要贝叶斯工作流程的原因很多，不仅仅是贝叶斯推断。

贝叶斯方法本身就是不断试错的过程。

计算是贝叶斯方法面临的一项挑战，但是为了获得更为可信的推断，经常不得不完成各种步骤，哪怕需要更多的时间。这其中就包括：拟合更简单的或替代的模型、做不太准确但速度更快的近似计算、摸索拟合过程等。

数据的动态性需要不断调整模型。

在困难的问题中，通常无法提前知道想要拟合什么样的模型。即便在极少数情况下预先选择了可接受的模型，也通常会在增加数据时、或需要对数据进行更详细解释时，调整这些模型。

模型对比和理解的需要。

即便数据是静态的，并且我们知道要拟合的模型，拟合也没有出现问题，我们依然希望能够深入理解模型及其与数据的关系，而这种理解往往必须通过模型比较才能实现。

模型自身不确定性的需要。

有时不同模型会得出不同的结论，但没有一个是明显占优势的。在这种情况下，展示多个模型有助于说明模型选择的不确定性。

1.3 工作流程及其与统计理论和实践的关系¶

工作流程在不同上下文中具有不同的含义。就本文而言，工作流程比『案例分析』更通用一些，但又不如『方法』那样足够精确（见图 2 ）。实际上，我们受到了计算中有关工作流程想法的影响。这些想法包括一些统计上的进展（例如 tidyverse ），有些进展可能不是特别贝叶斯，但具有类似的体验式学习的感觉（ Wickham and Groelmund, 2017 ）。机器学习的许多最新进展都具有类似即插即用的感觉：它们易于使用、易于试验，并且让用户有一种健康的感觉，即拟合模型是一种从数据中学习某些东西的方式，而该方式可能无需表达对某种概率模型或统计假设的承诺。

图 2 显示了我们对统计方法学发展的看法，即『 范例 –> 案例分析 –> 工作流程 –> 方法 –> 理论 』的编码过程。并非所有方法最终都能达到数学抽象级别，但纵观统计学历史，我们已经看到了很多在特定范例、风格化案例、新问题中的工作流程等背景下开发的方法，其中很多在可能的情况下，都做了形式化、编码和理论研究。

图 2：统计方法学的基础工作过程，表示新想法首先出现在范例（ Example ）中，然后被形式化为案例分析（ Case Study ），编制成工作流程（ Workflow ），进一步被赋予一般性的通用实现方法（ Method ），最终形成理论（ Theory ）。

理解图 2 的一种方法是沿着这条路径从左向右地理解统计史上的一些重要思想。有许多想法最初是作为黑客或与统计相邻的工具开始的，最终被形式化为方法并进入统计学的核心。

例如：

分层建模（ Multilevel modeling ）

分层建模扩展了模型，将关于先验的推断纳入到完全的贝叶斯框架中，并被形式化为对先验分布的经验贝叶斯估计。

探索性数据分析（ Exploratory data analysis ）

探索性数据分析可以理解为一种预测性模型检查的形式（ Gelman, 2003 ）。

正则化方法（ Regularization methods ）

诸如 Lasso ( Tibshirani, 1996 ) and Horseshoe ( Piironen et al., 2020 ) 等正则化方法已经取代了回归任务中的临时变量选择工具。

非参数模型（ Nonparametric models ）

高斯过程（ O’Hagan, 1978；Rasumussen and Williams, 2006 ）等非参数模型可以被认为是核平滑等程序的贝叶斯替代品。

在上述每一个案例以及其他许多案例中，这些方法都按照模块化、可用性的路线，被扩充到了统计方法论的框架中。

『工作流程』一词已逐渐在统计和数据科学中使用；参见 Liu et al., 2005；Lins et al., 2008；Long, 2009；Turner and Lambert, 2015。工作流程的相关思想在软件开发和其他信息学领域中广为流传；最近从业者中的讨论包括 Wilson et al., 2014 and 2017。应用统计（不仅仅是贝叶斯统计）变得越来越计算化和算法化，这将工作流程置于了统计实践的中心位置（参见 Grolemund and Wickham, 2017；Bryan, 2017 ； Yu and Kumbier, 2020 ），在应用领域也是如此（例如，Lee et al., 2019 讨论了心理学研究中的建模工作流程）。

贝叶斯工作流程已被 Savage , 2016、Gabry , 2019 和 Betancourt et al., 2020a表达为一个一般性概念。 Gelman , 2011 单独谈到了贝叶斯工作流程的几个组件，但没有以连贯方式进行讨论。此外，也有针对特定问题开发的贝叶斯工作流程，如 Shi and Stevens , 2008 以及 Chiu et al., 2017。

在本文中，我们介绍了贝叶斯工作流程的几个方面，希望这些内容最终能够进入日常实践和自动化软件。我们使用概率编程语言 Stan （ Carpenter et al., 2017 ； Stan Development Team, 2020 ）建立了大部分工作流程，但类似想法也适用于其他计算环境。

1.4 组织贝叶斯工作流程的诸多方面¶

统计工作流程的教科书表达通常是线性的，不同路径对应不同的问题情况。例如：医学临床试验通常从样本量的计算和分析计划开始，然后是数据收集、清理和统计分析，最后是报告 \(p\) 值和置信区间。经济学中观察性研究 可能从探索性数据分析开始，为选择变量的转换形式提供信息，后续是一系列回归分析，然后是一组替代性的模型指定和稳健性研究。

本文中讨论的统计工作流程比教科书和研究文章中介绍的普通数据分析工作流程更加复杂。额外的复杂性主要来自以下几个地方，并且在高级别的工作流程中有许多子工作流程：

（ 1 ）复杂的模型本身就需要更多的计算和评估。

复杂模型的拟合计算本身就很困难，需要进行相当数量的实验来解决从后验分布的计算、逼近或预测模拟问题，同时检查计算算法是否按照预期执行。

（ 2 ）复杂问题意味着更复杂的模型。

对于复杂的问题，我们希望想得到一个更复杂的通用模型（例如：包括相关性、层次结构和随时间变化的参数等特征），因此我们通常从已知的简单模型开始，其在计算上更容易，然后逐渐添加特征，增加模型的复杂性。

（ 3 ）数据的动态性意味着随时要调整模型。

与此相关，我们同时需要考虑数据动态性问题，常见情况有：数据收集正在进行中、或从相关数据集中动态抽取数据，例如民意分析中的新调查或药物试验中其他实验的数据。增加新数据往往需要模型扩展，以允许参数变化或扩展函数形式，例如线性模型一开始可能很适合，但随着新数据的增加而崩溃。

（ 4 ）为了更好地理解模型需要做大量模型比较工作。

除了拟合和模型扩展的挑战之外，通常可以通过与替代模型的推断进行比较来最好地理解模型。因此，我们的工作流程包括用于理解和比较多个能够拟合相同数据的模型的工具。

统计是关于不确定性的。除了数据和模型参数中常见的不确定性之外，我们还常常不确定是否正确地拟合了模型、不确定如何最好地建立和扩展模型、以及对它们理解和解释的不确定。一旦超越了简单的初始设计和分析，工作流程就会变得复杂和混乱，我们工作的重点会在数据探索、实质性理论、计算和结果解释之间来回转移。因此，任何对工作流程步骤的组织尝试都会显得过于简单，而且许多子工作流程也非常复杂，足以作为独立论文或书籍的单独章节。

我们在此讨论了工作流程的许多方面，从实践（尤其是可用时间、计算资源、错误惩罚的严重程度等）角度考虑，从业者可能会选择一些捷径。这些捷径往往会使结果更加难以解释，但必须意识到：近似计算总不会比压根儿不拟合模型更糟糕（近似可以定义为即便没有计算时间限制，也无法给出后验分布的准确估计）。我们在此描述统计工作流程的目的之一，也是希望明确地将各种捷径理解为完整工作流程的近似值，让从业者就如何投资有限时间和精力做出更明智的选择。

1.5. 文章的目的和结构¶

应用统计学中有各种各样的隐性知识，并不总能将其纳入已发表的论文和教科书。本文旨在公开其中一些想法，以改进应用型贝叶斯分析，同时为未来理论、方法和软件研发提出建议。

我们的目标受众是：

(a) 应用贝叶斯统计的从业者，尤其是 Stan 等概率编程语言的用户。

(b) 面向上述用户的方法研发和软件开发人员。

(c) 贝叶斯理论和方法的研究人员，因为我们认为工作流程的许多方面都没有得到充分研究。

在本文其余部分，我们将逐步介绍贝叶斯工作流程的各方面，如图 1 所示，从拟合模型之前要完成的步骤（第 2 节）开始，通过拟合、调试和评估模型（第 3-6 节 )，然后修改模型（第 7 节），并理解和比较一系列模型（第 8 节）。

第 10 节和第 11 节通过两个例子来示范这些步骤，一个例子将特征逐步添加到高尔夫推杆模型中，另一个例子通过一系列调查来解决拟合行星运动的简单模型。前者展示了新数据如何推动模型的优化和改进，同时说明了模型扩展时出现的一些意外挑战。后者展示了计算中存在的挑战给建模带来的困难。两个例子都没有使用贝叶斯工作流程的全部要素，但至少表明，将贝叶斯模型开发的许多方面系统化会有好处。

我们在第 12 节中梳理了一些一般性讨论，并对工作流程可能存在的批评意见做出了回应。

2 在拟合模型之前¶

2.1 选择初始模型¶

几乎所有分析的起点都是从适应先前所做的工作开始的，这些先前的工作包括教科书、案例、相关问题的论文等，有些类似软件工程中设计模式的概念。使用之前的模型并对其进行更新，是有效进行数据分析的捷径。通过对模型模板的结果分析，可以知道模型向哪个方向丰化或简化更有空间。模型的样板可以节省模型构建和计算时间，但同时我们也要承认：那些需要理解结果的人可能存在认知上的负担。快捷方式对人类和计算机都很重要，其能够解释『为什么典型的工作流程是迭代的』（更多信息请参见第 12.2 节）。类似地，如果要尝试对计算机进行编程以自动执行数据分析，必须通过某种算法来构建模型，而这种算法的构建块就表示某种模板。尽管在统计教学领域中，『菜谱式分析』有很强的负面含义（主要讽刺那种只教工具和手段，而不教为什么使用和如何使用的教学方式），我们依然坚持认为模板可以作为更精细分析的起点和比较点。我们应该认识到理论并不是静态的，『科学理论的发展过程』与『统计模型的发展过程』并不相同（ Navarro, 2020 ）。

有时我们的工作流程从一个简单模型开始，其目的就是为了稍后逐步添加特征（以便对不同的参数进行建模，包括测量误差、相关性等）。另外一些时候，我们可能会从一个大的模型开始，并打算在接下来的步骤中逐步剥离一些特征，以努力找到一些简单易懂、但仍能捕获数据关键特征的东西。有时我们甚至会考虑对相同的数据采用『多种完全不同的方法』来建模，因此可能会存在多个可供选择的起点。

注解：这也许解释了大量阅读文献、掌握发展状态的必要性。

2.2 模块化构建¶

一个贝叶斯模型是由若干模块构建而成的，这些模块通常可被视为在必要时能够被替换的占位符。

例如，先用正态分布对数据建模，然后用 长尾分布 或 混合分布 替换它；将回归函数建模为线性函数，然后将其替换为 非线性样条 或 高斯过程 ；可以先将一组观测视为精确的，然后添加一个测量误差的变量构成新模型；可以从弱先验开始，当发现后验推断包含不切实际的参数值时，再调整和增强先验。

将模块视为占位符可以减轻模型构建过程的一些压力，因为你可以随时返回并根据需要概括或添加信息。

模块化构造的想法与统计文献中的长期传统背道而驰。在传统中，整个模型被统一命名，每次对现有模型进行微调时都会给出一个新的模型名称。而给模型的各模块命名，其实更容易地看到不同模型之间的联系，并使其适应给定项目的特定要求。

2.3 对观测数据进行缩放和转换¶

我们希望观测变量都能够被给出实际和道德原因的解释。这导致需要其尽量在自然尺度上描述，并将其建模为相互独立的（如果可能的话），或者需要这些观测变量具有可解释的依赖结构，这有助于使用信息性先验（ Gelman, 2004 ）。

不过剥离出尺度以构造一个与尺度无关的新变量，也是很有帮助的。

例如，在药理学问题中（ Weber et al., 2018 ），预计有一个变量的测量尺度约为 \(50\) ；遵循缩放原则，我们可能会在 \(\log (\theta/50)\) 上建立一个模型，从而在对数尺度上的 \(0\) 对应了一个可解释的值（即原始尺度上的 \(50\) ），而在对数尺度上 \(0.1\) 的差异，则对应于原始尺度上增加或减少了约 \(10 \%\) 。

这种转换不仅为了便于解释；它还为分层建模做好了准备。当我们构建更大的模型（例如：通过合并来自其他患者组的数据或其他药物的数据）时，允许观测变量因组而异非常有意义，并且部分池化在尺度无关的变量上更为有效。

例如，毒理学模型需要研究每个人的肝脏体积。我们没有直接将分层模型拟合到肝脏体积上，而是对肝脏体积做了尺度缩放的变换，将其建模为相对于体重的比例；我们预计这些尺度无关的因子在不同人之间的差异较小，因此与直接对绝对的肝脏体积建模相比，缩放以后拟合的模型可以进行更多的部分池化。缩放变换是一种能够促进有效分层建模的好方法。

在许多情况下，我们可以通过对数变换、 \(\text{\text{logit} }\) 变换、标准化或者归一化来粗略地将参数置于单位尺度上。如果中心和尺度本身是根据数据计算出来的，就像在 rstanarm 中为回归系数的默认先验所设置的那样（ Gabry et al., 2020a ），则我们可以将其视为分层模型的近似，其中中心和尺度是从数据中估计的超参数。

更复杂的转换也可以使参数更易于解释，从而促进先验信息的使用；Riebler et al. (2016) 给出了一类空间相关模型的例子，Simpson et al., (2017) 更普遍地考虑这个想法。

注解：本节建议在做拟合之前，对原始观测数据进行分析，对于能够做尺度变换的变量可以尝试采用对数变化、\(\text{logit}\) 变换、标准化变化或归一化变换等方式，对原始数据进行预处理。在 Martin, 2016 的书中也有标准化和归一化等预处理的类似论述。

2.4 尝试一下先验预测性检查¶

先验预测性检查是在生成模型背景下理解先验分布含义的有用工具（ Box, 1980 ；Gabry et al., 2019；有关如何使用先验分布的详细信息，参见第 7.3 节）。特别是，由于先验预测性检查使用来自模型的模拟数据而不是真实的观测数据，因此能够提供了一种『无需多次使用数据即可优化模型』的方法。

图 3 显示了逻辑回归模型的一个简单的先验预测性检查情况。模拟表明，随着模型中预测变量数量的增加，各个系数的先验（即便是独立的）具有不同的含义。这是回归模型中的普遍现象，随着预测变量数量的增加，如果我们想让模型远离极端预测，就需要模型系数具有更强的先验（或足够多的数据）。

图 3 显示了逻辑回归模型的简单先验预测性检查。模拟数据表明，随着模型中预测变量数量的增加，模型中的各系数表现出不同的含义。这是回归模型中的普遍现象，随着预测变量数量的增加，如果想让模型远离极端预测，我们需要模型系数具有更强的先验（或足够多的数据）。

图 3：演示使用先验预测性检查来了解模型的非显著特征。上图对应于具有 \(100\) 个数据点和 \(2\)、\(4\) 或 \(15\) 个二值预测变量的逻辑回归模型。在每种情况下，所有回归系数都被赋予独立的 \(\text{Normal}(0,1)\) 先验分布。我们对每个模型执行先验预测性检查，首先从先验的系数向量 \(\theta\) 抽取 \(1000\) 次模拟，然后自逻辑回归模型中生成相应的模拟数据集 \(y\)，对模拟数据集求均值得到响应均值统计量 \(\bar y\) 。图中每个柱状条都显示了此统计量的先验预测性分布，即 \(1000\) 次模拟得到的 \(\bar y\) 。当模型中的预测变量数量很少时，这种先验预测性分布会展开，表明该模型与范围广泛的数据体系兼容。但是随着预测变量数量的增加，后验预测性分布变得集中在 \(\bar y = 0\) 或 \(1\) 附近，这表明模型各系数的弱先验意味着这个特定结果变量的强先验。如果我们想要一个 \(\bar y\) 上更温和的先验预测性分布，则系数的先验需要强烈集中在零附近（即具有更强的先验）。

一种有用的作法是先考虑结果变量的先验，然后推导出相应参数的联合先验（ Piironen and Vehtari, 2017 ； Zhang et al., 2020 ）。更一般地说，联合先验允许我们控制较大参数集的整体复杂性，这有助于生成更合理的先验预测，而独立先验很难或不可能实现这样复杂的预测。图 4 显示了对高斯过程模型的三个先验分布选择进行先验预测性检查的示例（ Rasmussen and Willams, 2006 ）。这种模拟数据生成方法和图形比较方法在处理任何模型时都很有用，特别是在模型设置（可理解为上下文）不清楚或模型很复杂时几乎必不可少。

图 4：高斯过程的先验预测性分布。图中分别为平方指数协方差函数具有不同幅度参数 \(τ\) 和长度尺度参数 \(l\) 时的高斯过程先验预测性分布图。（ Gelman et al., 2013 ）。

先验预测性模拟的另一个好处是，有助于从业人员获取或引出关于可测量兴趣量的专家先验知识，这比直接向领域专家征求 “关于不可观测模型参数的意见” 更为容易一些（ O'Hagan et al., 2006 ）。

最后，即使跳过计算先验预测性检查这一步骤，考虑我们选择的先验将如何影响假设的模拟数据集可能也会非常有用。

2.5 生成式和部分生成式模型¶

完全的贝叶斯数据分析需要一个生成模型，即所有可观测数据和模型参数的联合概率分布。这一点与贝叶斯推断有着极其微妙的不同：

贝叶斯推断：贝叶斯推断目的是推断得到后验分布，而根据贝叶斯定理，这只需要指定先验和数据的似然函数即可，不必指定生成模型，不同的生成模型可能具有相同的似然。
贝叶斯数据分析：贝叶斯数据分析要求有生成模型，并且能够进行预测性地模拟和模型检查（见第 2.4、4.1、4.2、6.1 节、6.2 节 ）
贝叶斯工作流程：贝叶斯工作流程的设置会同时考虑一组可能的生成模型，以进行评估、对比和选择等。

举一个简单的例子，假设有数据 \(y \mid \text{Binomial}(n,\theta)\)，其中 \(n\) 和 \(y\) 是可观测变量，我们希望对 \(\theta\) 进行推断。

对于贝叶斯推断而言

如果仅仅做贝叶斯推断，则数据是用通过什么方式生成的无关紧要。相对于估计 \(\theta\) 的目的而言，『固定 \(n\) 次采样（二项式采样）』或者『采样直到出现指定次数成功（负二项式采样）』所对应的似然是等效的，因为两者之间的不同仅由一个 \(y\) 和 \(n\) 决定的乘性因子决定，和 \(\theta\) 无关。

对于贝叶斯数据分析或工作流程而言

如果想模拟来自预测模型的新数据，则上述两种模型是不同的。因为二项式模型产生具有固定值 \(n\) 的重复项，而负二项式模型产生具有固定 \(y\) 值的重复项。这两种不同的生成模型下，先验和后验预测性检查看起来会有所不同（ 第 2.4 和 6.1 节 ）。

上例并不是说贝叶斯数据分析的方法一定更好；对于生成模型的假设有助于提高推断效率，但也很可能出错，而正是这个原因驱动了贝叶斯工作流程的大部分要素。

当然，在贝叶斯分析中使用不完全生成的模型也很常见。

例如，在回归任务中，通常会对给定预测变量 \(x\) 的结果 \(y\) 建模，而没有 \(x\) 的生成模型。另一个例子是带有删失的生存数据，其中通常不对删失过程建模。

在对此类模型执行预测性检查时，要么需要以观测到的预测变量为条件，要么扩展模型以允许对预测变量的新值进行采样。当然，模型的某些部分也可能没有随机生成过程，例如： \(x\) 是通过确定性的实验设计选择而得。

从生成模型角度思考，可以帮助阐明从观测中所学信息的局限性。例如，我们可能想对具有复杂自相关结构的时间过程建模，但如果实际数据在时间上相距很远，可能无法将该模型与具有几乎独立误差的更简单过程区分开来。

此外，使用不正确先验的贝叶斯模型不是完全生成的，因为它们没有数据和参数的联合分布，并且不可能从先验预测分布中采样。当我们确实使用了不正确先验时，可以将其视为 “为了建立一个在参数和数据上具有适当联合分布的完整贝叶斯模型而设置的” 占位符或沿着路径的某个步骤。

在应用工作中，复杂性通常来自于合并了不同数据源。例如，使用州和全国民意调查数据为 \(2020\) 年总统选举拟合贝叶斯模型，部分池化面向基于政治和经济基本面的预测（ Morris、Gelman and Heidemanns, 2020 ）。该模型包括一个针对州和国家舆论中潜在时间趋势的随机过程。使用 Stan 拟合模型能够产生后验模拟，进而用于计算选举结果的概率。基于贝叶斯模型的方法表面上看起来有些类似于 Katz , 2016 描述的投票聚合法，聚合法也通过随机模拟来得到不确定性的度量。而不同之处在于：贝叶斯模型还可以向前运行以生成预测的投票数据；贝叶斯模型不仅是一个数据分析程序，还能够为国家和州级的舆论提供了一个概率模型。

更一般地思考，我们可以考虑一个从最少到最多生成模型的过程。

层次一：完全非生成式方法

其中一个极端是完全的非生成式方法，它们被简单地定义为数据摘要，根本没有数据模型。

层次二：经典频率主义统计模型

其特征在于给定参数 \(\theta\) 时数据 \(y\) 的概率分布 \(p(y; \theta)\) ，但没有 \(\theta\) 的概率分布（即人为 \(\theta\) 是一个确定的固定值）。

层次三：部分生成式的贝叶斯模型

在 \(y\) 和 \(\theta\) 上构建生成模型，但模型中还包括额外的未建模变量 \(x\) 以及样本大小、设计设置和超参数等要素，我们将此类模型写成 \(p(y,\theta \mid x)\) 。

层次四：完全生成式模型

最后一步可以考虑一个完全的生成式模型 \(p(y,\theta,x)\) ，该模型将会对所有数据和模型参数建模，得到一个不会 “遗漏” \(x\) 的生成模型。

在统计工作流程中，我们可以在上述阶梯层次中上下移动，例如从未建模的数据概略算法开始，然后将其形式化为概率模型；或者从概率模型推断开始，将其视为基于数据的估计，并以某种方式对其进行调整以提高性能。在贝叶斯工作流程中，我们可以将数据移入和移出模型，例如采用未建模的预测变量 \(x\) 并允许其具有测量误差，以便模型包含新级别的隐变量（ Clayton, 1992; Richardson and Gilks, 1993 ）。

3 拟合一个模型¶

历史传统中，贝叶斯计算一致局限于使用解析方式和正态近似的组合来执行的。直至 1990 年代，使用 Gibbs and Metropolis 算法对任意模型执行贝叶斯推断成为可能（ Robert and Casella, 2011 ）。

当前用于拟合开放式贝叶斯模型的最先进算法包括变分推断法（ Blei and Kucukelbir, 2017 ）、序贯蒙特卡罗法（ Smith, 2013 ）和 Hamiltonian Monte Carlo ( HMC ) （ Neal, 2011; Betancourt, 2017a ）等。其中变分推断是对期望最大化 (EM) 算法的推广，在贝叶斯语境中，变分推断可被视为提供了对后验分布的一种快速但可能不准确的近似。变分推断（特别是摊销变分推断 ）是用于计算密集型模型（例如深度神经网络 ）的当前标准方法。序贯蒙特卡洛方法是 Metropolis 算法的推广，可以应用于任何贝叶斯计算。 HMC 是 Metropolis 的另外一种推广，由于它可以使用梯度计算在连续概率空间中进行高效移动，因此使用率比较高。

注解：近年出现了 MC Dropout 等直接求解预测结果不确定性的模型，并成为研究热点之一。

在本文中，我们专注于使用 HMC 及其变体来拟合贝叶斯模型（ HMC 在 Stan 等概率编程语言中均有实现）。虽然类似原则也适用于其他软件和算法，但可能存在细节上的差异。

要在贝叶斯工作流程中安全地使用推断算法，算法必须提供足够强大的诊断功能来确定计算的可靠性，或提醒我们何时不可靠？在本文中，我们讨论了一些面向 HMC 方法的一些诊断工具，其中部分也适用于其他推断方法。

3.1 初值、适应和预热¶

除了最简单的情况外，马尔可夫链模拟算法一般会包含多个运行阶段：

阶段一：预热阶段

预热阶段旨在将模拟从可能不具代表性的初始值移动到更接近参数空间的区域。在该目标区域中，对数后验密度接近模型参数的期望值，这与信息论中的典型集有关（ Carpenter, 2017 ）。预热阶段存在初始值选择的问题。初值在渐近极限过程中并不重要，但在实践中很重要，因为错误的初值选择可能会威胁到结果的有效性和收敛的时效性。

阶段二：适应阶段

适应同样发生在预热期间，此时需要有一些程序来设置算法的调整参数；这可以利用从预热运行中收集的信息来完成。

阶段二：采样阶段

第三个阶段是采样，理想情况下一直运行到多个链混合（ Mixed ）为止（ Vehtari et al., 2020 ）。

由此看来，预热主要有三个目的：

目的 1：快速标记有问题的模型。

在拟合正在探索的模型时，可以快速标记出计算上有问题的模型。

目的 2：减少初始值带来的模型偏差。

在拟合正确指定的模型时，运行一段时间以减少初始值导致的偏差；

目的 3：提供参数调整的信息。

为设置调整参数提供有关目标分布的信息。

3.2 运行迭代算法需要多长时间？¶

我们同样希望在更大的工作流程中，考虑运行迭代算法所需的判决依据。

推荐的标准做法是至少运行到所有感兴趣的参数和量的 \(\hat R\) ( 链混合程度的度量）都小于 \(1.01\) （ Vehtari et al., 2020 ），同时可以监测多变量混合的统计量 \(R^*\) （ Lambert and Vehtari, 2020 ）。

有时，在建模早期提前停止是有意义的。例如，运行 MCMC 直到有效样本大小达到数千似乎是一个安全而保守的选择，或者与参数解释所需的精度相比，蒙特卡罗标准差很小，但这会需要很长时间，限制了探索阶段可拟合的模型数量。而且通常情况下，我们的模型会存在一些较大的问题（尤其是程序代码错误），这些问题在仅运行几次迭代后就会变得非常明显，没必要浪费剩余的计算时间。在这方面，为新编写的模型运行多次迭代有些类似于软件工程中的过早优化。对于最终模型，所需迭代次数取决于感兴趣量所需的蒙特卡罗精度。

计算中的另一个选择是利用并行机制，不局限于在多个内核上运行 \(4\) 个或 \(8\) 个独立链的默认设置。除了增加迭代次数外，还可以通过增加并行链的数量来获得有效的方差减少 （参见 Hoffman and Ma, 2020 ）。

3.3 近似算法和近似模型¶

贝叶斯推断通常涉及难以处理的积分，因此需要近似。马尔可夫链模拟是一种近似形式，随着模拟次数增加，理论误差接近于零。如果链是混合的，我们可以很好地估计蒙特卡罗标准差（ Vehtari et al., 2020 ），并且出于实际目的经常将这些计算视为准确的。

不幸的是，随着数据和模型变大，运行 MCMC 至收敛并不总是一个具备可扩展性的解决方案，因此需要更快的近似值。图 5 显示了速度和准确性之间的权衡。此图只是概念性的；在实际问题中，线的位置是未知的。在某些实际问题中，即使在较短的时间尺度上，近似算法的性能也有可能比 MCMC 差。

图 5：贝叶斯计算中近似算法和 MCMC 之间权衡的理想草图。如果任务是获得目标分布的最佳拟合，那么 MCMC 最终应该胜出。但如果考虑的是拟合一组模型，那么使用近似算法以便能够在模型空间中快速切换非常有意义。这些算法中哪一个表现更好取决于用户的时间预算以及两条曲线的交点。

根据在工作流程中所处的位置不同，我们对计算的后验有不同的要求。

在工作流程接近尾声时，我们正在仔细检查最佳尺度和精细特征，因此需要尽量准确的后验分布，这通常需要 MCMC 方法。
在工作流程开始时，我们可以经常根据后验的大规模特征做出建模决策，这些特征可以使用相对简单的方法进行较准确地估计，例如：经验贝叶斯法、线性化法或拉普拉斯近似法、嵌套近似法（如 INLA，Rue et al., 2009 ）等。甚至有时可以采用如期望传播的数据拆分近似方法（ Vehtari、Gelman、Siivola et al., 2020 ）、变分推断的寻找峰值近似方法（ Kucukelbir et al., 2017 ），或受惩罚的最大似然法等。关键是要使用合适的工具来完成这项工作，而不是试图用雕刻家的凿子敲掉挡土墙。

所有上述近似方法背后都有至少十年的实践经验、理论和诊断。没有一刀切的近似推断算法，但是当工作流程包含相对容易理解的组件（例如广义线性模型、分层回归、自回归时间序列模型或高斯过程）时，通常可以构造适当的近似算法。此外，根据所使用的特定近似值，Yao et al.(2018a) and Talts et al. (2020) 描述了通用诊断工具，可用于验证特定近似算法是否重现了特定模型的后验特性。

另一种观点是将近似算法理解为近似模型的精确算法。从这个意义上说，工作流程是抽象计算方案中的一系列步骤，旨在推断出一些最终的、未陈述的模型。更有用的是，我们可以将经验贝叶斯近似之类的东西视为一个依赖于数据点质量先验的先验分布。类似地，拉普拉斯近似可以被视为所需模型的依赖于数据的线性化，而嵌套拉普拉斯近似（ Rue et al., 2009 ; Margossian et al.,2020a ）使用线性化条件后验代替了假定的条件后验。

3.4 快速拟合，快速失败¶

另外一个重要的中间目标是在拟合“坏”模型时能够快速失败。这可以被认为是一种捷径，可以避免浪费大量时间对一个“坏”模型进行接近完美的推断。有大量关于近似算法以快速拟合所需模型的文献，但很少涉及旨在将尽可能少的时间浪费在最终将放弃的模型上的算法。我们认为根据此标准评估方法很重要，特别是“坏”模型通常更难以拟合，从而有可能浪费更多的时间。

对于一个简单的理想化示例，假设你是几个世纪前的天文学家，根据 \(10\) 个带误差的测量数据点来拟合行星轨道椭圆。图 6a 显示了可能出现的数据类型，几乎任何算法都可以很好地拟合。例如，你可以采用不同的 \(5\) 个点集并为每个点拟合精确的椭圆，然后取这些拟合的平均值。或者你可以将一个椭圆拟合到前五个点，稍微扰动它以拟合第六个点，然后稍微扰动它以拟合第七个点，依此类推。或者你可以实现某种最小二乘算法。

图 6：“快速拟合，快速失败”的需求说明：（ a ）理想化代表行星轨道测量值的数据，可以拟合为具有测量误差的椭圆，(b) 被死星扰动的假设轨道的测量值。在第二个示例中，将单个椭圆拟合到数据中将具有挑战性，但在任何情况下，我们都对椭圆拟合没有特别的兴趣。我们希望任何将椭圆拟合到第二个数据集的尝试都会快速失败。

现在假设一些死星出现并改变了轨道，在这种情况下，我们故意选择一个不切实际的例子来创建模型和数据之间的严重差异，因此你的 \(10\) 个数据点看起来像图 6b 。在这种情况下，收敛将更难实现。如果你从拟合前五个点的椭圆开始，则很难采用任何一组小扰动来使曲线拟合系列中后面的点。但更重要的是，即使可以获得最小二乘解，其对应的任何椭圆都与数据拟合的非常糟糕。这说明椭圆不是一个合适的模型，如果将椭圆拟合到这些数据，你应该希望拟合快速失败，以便快速转向更合理的方法。

这个例子具有一种困难统计计算的通用模式，即拟合数据的不同子集会产生非常不同的参数估计。

4 使用人为构造的伪数据发现和理解问题¶

验证计算的第一步是检查模型是否在可接受的时间范围内切实地完成了拟合过程，并且收敛诊断是合理的。在 HMC 上下文中，主要的诊断手段有：不存在发散的转移、\(\hat R\) 诊断指标接近 \(1\) 、聚集趋势的充分有效样本量、尾分位数和能量（ Vehtari et al., 2020 ）。不过，上述诊断方法无法防止能够正确计算的、但代码对应的模型与用户预期不符的概率程序。

我们用于确保统计计算完成得相当好的主要工具是：『将模型拟合到某些实际数据，并检查拟合是否良好』。为此目的，真实数据可能很尴尬，因为建模问题早期可能会与计算问题发生冲突，并且很难判断到底是计算问题还是模型问题。为解决这个问题，可以先将模型拟合到模拟数据来探索模型。

4.1 伪数据模拟¶

在已知真实参数的受控环境中工作，可以帮助我们理解“数据模型和先验”、“可以从实验中学到什么”以及“推断方法的有效性”。基本思想是检查程序在拟合伪数据时，是否恢复了正确的参数值。通常，我们选择在先验上看起来合理的参数值，然后模拟与原始数据具有相同大小、形状和结构的伪数据集。接下来我们将模型拟合到伪数据上来检查几件事：

（ 1 ）用伪数据来检查模型是否能够自证合理性

我们检查的第一件事不是严格地计算，而是在设计方面。对于所有参数，检查观测数据是否能够提供超出先验的额外信息。该过程使用固定的已知参数，从模型中模拟一些伪数据，然后查看我们的方法是否接近重现了已知事实。我们可以通过查看后验的点估计以及后验区间覆盖范围来检查。

如果伪数据检查失败了，即在某种意义上，推断与假设的参数值不接近，或者模型某些组件没有根据数据获得任何信息（ Lindley,1956 ；Goel and DeGroot, 1981 ），则建议分解该模型，让其变得越来越简单，直到模型能够工作为止。然后，在那里尝试识别问题，如第 5.2 节所示。

（ 2 ）用伪数据来检查明显的模型错误或不一致

第二件需要检查的事情是：真实参数是否能够大致恢复到拟合后验分布所隐含的不确定性范围内。原则上，数据不能为参数提供信息的情况不可能发生，如果它在我们的模拟检查中发生了，那一定是哪里出了问题。通常我们不会在运行了单个伪数据模拟并计算了相关后验分布后，就声明一切正常。我们将在下一节看到，此时需要更精细的设置。不过通过单个模拟运行通常会发现非常明显的错误，其中最明显的就是真实参数和推断结果不一致。例如，如果代码中有错误并且拟合了错误的模型，这通常会从参数恢复的灾难性失败中表现出来。

（ 3 ）用伪数据来了解模型的在参数空间中的变化特性

第三件事是使用伪数据模拟来了解模型的行为在参数空间的不同位置会发生什么变化。从这个意义上说，统计模型可以包含许多关于数据如何生成的故事。如第 5.9 节所示，当指数分离得很好时，数据可以提供有关递减指数求和参数的信息，但当两个分量彼此接近时，数据就不那么重要了。这种取决于参数值的不稳定现象在微分方程模型中也比较常见，如第 11 节所示。再举一个例子，分层模型的漏斗型后验分布在瓶颈处与入口处存在较大不同。类似的问题也出现在高斯过程模型中，这取决于高斯过程的长度尺度和数据的分辨率。

所有这些都意味着伪数据模拟在数据的可预测参数空间中可能是特别相关的。这反过来又提出了一个两步的工作程序。在该工作程序中，首先将模型拟合到真实数据中，然后从得到的后验分布中提取参数样本并用于第二步的伪数据检查。这个工作程序的统计特性尚不清楚，但在实践中，我们发现非常有用，既可以揭示计算或模型问题，也可以在基于伪数据的推断确实再现了假定参数值时提供一些保证。

为了进一步实现此想法，可以提供一些能够导致工作程序给出错误答案的伪数据来破坏我们的方法。这种模拟和探索可能是深入理解推断方法的第一步，即使对于仅计划将该方法用于特定应用问题的从业者来说，也很有价值。它也可用于构建一组更复杂的模型以供稍后探索。

伪数据模拟是本文贝叶斯工作流程的关键组成部分，因为它是整个流程中，唯一一处能够直接检查隐变量推断可靠性的地方。根据隐变量的定义，在将模型拟合到真实数据时，不会观测到隐变量。因此，只能评估模型是否拟合了观测数据。如果我们的目标不仅是预测，还需要估计隐变量，那么检查伪数据会有很大帮助。对于过度参数化的模型尤其如此，这些模型中往往不同的参数值能够产生可比较的预测结果（例如第 5.9 节和 Gelman et al., 1996 ）。一般来说，即使模型拟合良好，我们也应该谨慎地对估计的隐变量给出结论。将模型拟合到模拟数据，有助于更好地理解：在已知基本事实的受控环境中，模型究竟可以了解还是不了解潜在过程的哪些因素。一个模型能够对从其生成的伪数据做出很好的推断，并不能保证该模型对真实数据的推断是合理的。但是，如果该模型无法对伪数据进行良好的推断，那么期望该模型能够对真实数据进行合理推断就更加无望了。伪数据模拟提供了了解潜在过程的理想上界。

4.2 基于模拟的较正¶

如图 7 所示，将贝叶斯推断结果（后验分布）与（真实的）点数据进行比较时，存在一个正式的、有时甚至是实际的问题。

图 7：后验分布与假设的真实参数值做比较。在将模型拟合到模拟数据时，我们可以检查后验样本（蓝色直方图）是否与真实参数值（红线）相符。在场景 1 （左图）中，后验以真值为中心，这表明拟合是合理的。在场景 2 （中图）中，真实参数出现在后验分布的尾部。目前尚不清楚这是否表明我们的推断存在错误。在场景 3（右图）中，后验存在多峰，很明显，将后验与单个点进行比较无法验证推断算法。更全面的方法（比如基于模拟的校准）可以教会我们很

即使计算算法工作正常，使用单个伪数据模拟来测试模型也不一定就有效。这不仅是因为在一次模拟中任何事情都可能发生（随机抽取有 \(5\%\) 的机会位于 \(95\%\) 不确定性区间之外），而且还因为贝叶斯推断通常只会在对先验求平均时才会进行校准，单一值无法进行校准。此外，参数恢复失败可能不是因为算法失败，而是因为观测数据无法提供足以更新某些参数先验的不确定性信息。在先验和后验近似单峰时，如果选择的参数值来自先验的中心附近，则可以预期基于伪数据的后验区间会出现过度覆盖。

比我们在第 4.1 节中介绍的方法更全面的方法是基于模拟的校准（ SBC )； Cook et al., 2006 ; Talts et al., 2020 ）。在该方案中，模型参数取自先验；然后根据这些参数值模拟数据；然后模型拟合数据；最后将获得的后验值与用于生成数据的模拟参数值进行比较。通过多次重复此过程，可以检查推断算法的一致性。其想法是：通过在一系列在从先验中抽取的参数模拟的数据集上进行贝叶斯推断，我们应当能够恢复先验。基于模拟的校准对于评估近似算法与理论后验的匹配程度非常有用，即使在后验不可处理的情况下。

虽然在许多方面优于针对真实点数据进行的基准测试，但基于模拟的校准需要多次拟合结果，这会产生大量的计算成本，特别是在不使用并行化的情况下。在我们看来，基于模拟的校准和基于真实值的基准测试是光谱的两端。粗略地说，一个单一的真实值基准测试可能会标记出严重的问题，但并不为任何事情做保证。随着做更多实验，才有可能越来越精细地发现计算中存在的问题。如何通过更少的抽样来理解 SBC 是一个开放的研究命题，我们预期：放弃从更有设计的先验中随机抽样，将使该方法更为有效，尤其是在参数数量相对较少的模型中。

SBC 的一个严重问题是，它与大多数建模者习惯于将先验指定的比所需宽度更宽相冲突。弱信息先验的保守性质可能导致 SBC 期间模拟的数据集变得极端。加布里 et al., 2019 模拟了一个假空气污染的数据集，其中污染甚至比黑洞还密集。这些极端数据集可能导致在真实数据上运行良好的算法出现严重失败。但这并不是计算问题，而是先验问题。解决此问题的一种可能方法是确保先验收紧。然而，一个更务实的想法可能是保持先验并使用真实数据计算合理的参数值，这可以通过粗略估计或计算真实后验来实现。我们建议稍微扩大一下上述合理参数值的估计范围后，再将其用做 SBC 的先验。这将确保所有模拟数据在模型允许的情况下尽可能真实。

4.3 使用构造的数据进行实验¶

理解模型的一个好方法是将其拟合到从不同场景模拟的数据中。

一个简单的例子，假设我们对拟合自其他类型分布数据的线性回归模型的统计特性感兴趣。首先，我们可以指定数据 \(x_i, i = 1,...,n\) ，从先验分布中抽取系数 \(a\) 、 \(b\) 以及残差标准差 \(\sigma\) ，并模拟来自 \(y_i \sim \text{Normal}(a + bx_i ,\sigma)\) 的数据，然后将模型拟合到上述模拟数据。重复这个过程 \(1000\) 次后，我们可以检查区间估计的覆盖范围（可以看出，这是一个基于模拟的校准版本）。

我们可以使用不同的假设模拟数据来拟合相同的模型，例如从 \(t_4\) 分布而不是正态分布中抽取独立的 \(y_i\) 数据点。不过这将使基于模拟的校准出现失败，因为正在拟合一个错误的模型，但有趣的问题是：这些推断会有多糟糕？例如，可以使用 SBC 模拟来检查模型参数后验分布的 \(50\%\) 可信区间和 \(95\%\) 可信区间覆盖范围。

对于更缜密的示例，我们会执行一系列模拟，以便更好地了解观察性研究中的假设和偏差校正。例如：我们可以从一个 \(500\) 个学生参加的期中和期末考试场景开始。首先从 \(\text{Normal}(50,20)\) 分布中抽取学生的真实能力 \(η_i\) 来模拟数据，然后将两场考试的分数 \(x_i,y_i\) 绘制为独立的 \(\text{Normal}( η_i,10)\) 分布。这导致两个分数的相关性为 \(\frac{20^2}{20^2+5^2} = 0.94\) ；我们设计如此高相关值的模拟，主要是为了让图形中的模式表现地更明显。图 8a 显示了结果以及基础回归线 \(E(y \mid x)\) 。

然后，我们构建了一个假设的随机实验，在期中考试之后的教学改进会导致任意学生的期末考试成绩增加 \(10\) 分。我们为每个学生提供接受改进的平等机会。图 8b 显示了模拟数据和基础回归线。在这个例子中，我们可以通过简单地取两场考试之间的分数差异来估计改进效果。对于这组模拟数据，差异的估计值为 \(10.7\)，标准差为 \(1.8\) 。或者可以对期中考试的分数进行回归调整，得出 \(9.7 ±0.6\) 的估计值。

图 8：关于 500 名学生的教育干预效果模拟数据。

图 9：替代的模拟，其中教学改进工作分配不平衡，表现较差的学生更有可能接受改进。

接下来我们考虑一种不平衡分配机制，其中接受教学改进的概率取决于期中分数：\(\text{Pr}(z = 1) = \text{\text{logit} }^{−1}((x −50)/10)\)。图 9a 显示了 \(200\) 个学生的模拟治疗分配，图 9b 显示了据此模型模拟的考试分数。隐含的回归线与之前相同，因为模拟仅更改了 \(z\) 的分布，并没有更改模型 \(y \mid x,z\) 。在新设计下，优先对待表现较差的学生，因此按照之前的模型，对期末考试成绩做简单的比较，将会给出一个较差的估计。对于这个特定的模拟数据，之前的模型会给出两场考试的平均成绩差异为 \(−13.8 ±1.5\)，这是一个可怕的推论，因为我们构造时的真正效果是 \(+10\)。好在对 \(x\) 做调整后的线性回归恢复了教育改进工作的真实效果，产生了 \(9.7 ±0.8\) 的估计值。

这个新的估计值对 \(x\) 的调整函数形式很敏感。我们可以通过模拟来自改进模型的数据来观察这一点，其中真实的改进效果为 \(10\)，但函数 \(E(y \mid x,z)\) 不再是线性的。在此情况下，通过像以前一样从 \(\text{Normal}(η_i,10)\) 中抽取出给定真实能力的期中考试分数，并转换期末考试的分数，如图 10 所示。我们再次考虑两个假设实验：图 10a 显示了完全随机分配的数据，图 10b 显示了使用图 9a 中的不平衡教育分配规则的结果。两个图均显示了隐含的回归曲线。

当我们现在拟合线性回归来估计教育效果时会发生什么？图 10a 中设计的估计是合理的：即使线性模型是错误的，并且产生的估计在统计上也并不完全有效，但设计中的教育平衡概率确保了平均而言，特定误差将会被抵消，估计值为 \(10.5±0.8\)。但教育不平衡的设计存在问题：即使在线性回归中调整 \(x\) 后，估计值仍为 \(7.3 ±0.9\)。

在本文的上下文中，此示例的重点是展示在不同条件下模拟统计系统如何让我们深入了解计算问题，而且更广泛地了解数据和推断。在这个特定例子中，可以通过考虑不同的教育效应、对不可观测项的选择、以及其他选项来进一步研究，并且通常情况下，可以无限期地考虑此类理论探索以解决可能出现的任何问题。

图 10：再次比较两种不同的教育分配方案，这次是在期末和期中考试分数之间的关系为非线性。在本文的上下文中，此示例的重点是展示在不同条件下模拟统计系统如何让我们深入了解计算问题，而且更广泛地了解数据和推断。在这个特定的例子中，人们可以通过考虑不同的教育方案效果、对不可观测项的选择和其他选项来更进一步，而且通常情况下，可以无限地考虑此类理论探索以解决可能出现的任何问题。

5 解决统计计算中存在的问题¶

5.1 统计计算的通俗理论¶

当你遇到计算问题时，通常是模型出现了问题（ Yao、Vehtari and Gelman, 2020 ）。许多收敛不佳的情况对应于参数空间中没有实质性意义的区域，甚至对应于一个毫无意义的模型。图 6 给出了参数空间不相关区域中的病态示例。另外一个基础的问题模型案例是在代码中存在错误、或者是对高斯或逻辑回归中的每个观测量分别使用不同的高斯分布截距，在这些情况下，无法从数据中获取任何有用信息。面对问题模型时，第一直觉不应该是在模型上投入更多的计算资源（例如，通过运行采样器进行更多迭代或减少 HMC 算法的步长），而是检查模型是否包含本质上的病态。

5.2 从简单和复杂的模型开始并在中间相遇¶

图 11 说明了一种常用的调试方法。起点时模型表现不佳，可能无法收敛、或无法在伪数据模拟中重现真实参数值、或无法很好地拟合数据、或产生不合理的推断结果。诊断问题的途径大致是向两个方向移动：

一是逐渐简化性能极差的模型，剥离不必要的特征直到其能够工作；

二是从一个简单易懂的模型开始，逐渐添加特征直到问题出现。

类似地，如果模型有多个组件（例如，一个微分方程以及关于其参数的线性预测器），先使用模拟数据确保每个组件可以单独拟合并完成单元测试，通常是明智的选择。

我们可能永远无法成功拟合最初打算拟合的复杂模型，要么是因为当前可用计算算法拟合太困难，要么是因为现有数据和先验信息不足，无法从模型中进行有用的推断，或者仅仅是因为模型探索的过程将我们引向了与原计划不同的方向。

图 11：建议的调试图。右下方的星号代表在尝试拟合所需的复杂模型时出现问题的场景。左上角的点代表在拟合各种简单版本时成功，右下角的点代表在拟合完整模型的各种简化时失败。虚线表示可以在适合的简单模型和不适合的复杂模型之间的某处识别问题的想法。来自 Gelman and Hill, 2007 。

5.3 应对需要很长时间才能拟合的模型¶

我们通常使用 HMC 来拟合模型，但它可能由于各种原因运行缓慢，其中包括：昂贵的梯度估计、高维参数空间、在后验中部分区域有效的 HMC 在其他区域却不合适等。计算缓慢通常是其他问题的先兆，它不仅表明 HMC 的性能不佳，也同时意味着模型更难调试。

例如，最近在 Stan 用户组中收到一个关于分层逻辑回归的问询，其中包含 \(35,000\) 个数据点、 \(14\) 个预测变量和 \(9\) 个不同截距的批量数据，使用 rstanarm 的默认设置在 Stan 中运行几个小时后仍然未能完成运行。

我们给出了以下建议，没有特定的顺序：

（ 1 ）利用伪数据尽早发现问题。

从模型中生成模拟伪数据，并尝试将模型拟合到伪数据（第 4.1 节），错误指定的模型此时通常会很慢，而模拟的伪数据让我们不用考虑问题是否出在拟合不足，而是直接定性为模型问题。

（ 2 ）从小模型开始性价比更高一些。

大模型计算会很慢，所以应该从较小的模型开始入手构建。例如：首先拟合没有变化截距的模型，然后添加一批不同的截距，然后是下一批，依此类推。

（ 3 ）早期少一些迭代次数。

运行 \(200\) 次迭代，而不是使用 \(2000\) 次的默认值。通常问题比较明显时， \(200\) 次迭代足以提醒你发现问题。你当然也可以运行 \(2000\) 次迭代，但是当你仍在试图弄清楚发生了什么问题时，迭代次数太多没有意义，不如将其用在测试更多模型上。当 \(200\) 次迭代后 \(\hat R\) 仍然很大时，可以尝试运行更长的时间，但完全没有必要一开始就使用 \(2000\) 次迭代。

（ 4 ）尝试设置信息性先验。

在回归系数和分组级别的方差参数（第 7.3 节）上放置适度的信息性先验。

（ 5 ）考虑预测变量之间可能的交互作用。

有 \(14\) 个预测变量但其间毫无交互作用的加法模型让人感觉非常奇怪。此处的这个建议，与用户关心的性能问题似乎无关（事实上，添加交互作用只会增加计算时间），但应该清楚，我们的最终目标是做出预测或了解底层过程，而不仅是获得一组随意选择的预测变量，我们应当选择能够收敛的模型。

（ 6 ）先在在数据子集上拟合模型。

在将模型投入到所有数据之前，了解拟合过程并使其正常工作的一般性建议是现在数据子集上做拟合。由于数据量减小，计算时效性会得到提升，当在数据子集上顺利拟合后，在尝试在全数据上做进一步的工作，这往往会起到事半功倍的效果。

所有上述技巧的共同主题是将任何特定的模型选择视为一种临时考虑，并且一定要认识到：数据分析需要拟合许多模型而不是一个，以便对特定应用问题的计算和推断过程进行控制。

5.4 监控中间量¶

诊断模型问题的另一种有用方法，是在计算过程中保存中间量，并将该中间量与其他 MCMC的输出一起绘制可视化图件，例如使用 bayesplot （一个 R 包，Gabry et al., 2020a ）或 ArviZ （一个 Python 包， Kumar et al., 2019 ）。这些可视化方法可用于取代在代码中插入打印语句的陈旧方法，而且通常我们能够从可视化图表中学到更多的内容。

有时会出现的一个问题是马尔科夫链会卡在参数空间中的某个后验密度比较低的偏僻区域、或 MCMC 算法中无法回到对数后验密度接近期望的区域、或无法接近大部分后验质量所在的区域。正如我们将在第 11 节中说明的那样，在给定这些参数值的情况下，查看模型的预测以了解出了什么问题会很有帮助。但最直接的方法还是绘制出 “基于这些链中的参数值” 的期望数据，进而将参数的梯度转换为期望数据的梯度。这样做应该能够让我们进一步理解：参数的变化是如何被映射到后验相关区域所对应的期望数据上的。

5.5 堆叠以重新加权混合不良的链¶

在实践中，通常 MCMC 算法混合得很好。其他时候，模拟会快速移动到参数空间的不合理区域，表明可能存在错误的模型指定问题、或者观测数据缺乏信息或仅包含弱信息、或只是仅仅因为后验几何形态过于复杂。

但是中间情况也很常见，例如多个链混合起来很慢，但它们都在一般合理的范围内。此时可以考虑使用堆叠来组合模拟，并使用交叉验证为不同的链分配权重（ Yao、Vehtari and Gelman, 2020 ）。这会达到丢弃卡住的链的近似效果。堆叠的结果不一定等价（甚至是渐近等价）于完全贝叶斯推断，但它服务于许多相同的目标，特别适合在模型探索阶段，让我们向前迈进，花更多时间和精力在其他贝叶斯工作流程步骤上，而不会被拟合一个特定模型所困扰。

此外，当与轨迹图和其他诊断工具一起使用时，非均匀堆叠权重有助于以迭代方式了解应当将精力聚焦在什么地方。

5.6 如何应对具有复杂几何形态的后验分布？¶

在后验分布的复杂几何形态中，最常见的情况是存在多个峰。通常可以粗略地将 MCMC 中的多峰后验几何形态分为四类情况：

（ 1 ）类型一：可以显著区分的多峰后验形态

特点：多个峰中只有一个峰的质量比较大，其他峰的质量都接近于零。第 11 节 中有一个相应的例子。

处理方法：通过选择初值、添加信息性先验、对参数设置硬约束等来规避次要峰，或者通过近似估计每个峰的质量来修剪它们。

（ 2 ）类型二：可以显著区分的微对称多峰后验形态

特点：多个峰之间分割地比较清楚，各峰的概率质量相似，例如混合模型中的标签。

处理方法：标准做法是以某种方式限制模型以识别感兴趣的峰；参见 Bafumi et al., 2005 和 Betancourt, 2017b 。

（ 3 ）类型三：可以显著区分但质量不同的多峰后验形态。

特点：多个峰之间隔离地比较清晰，但各峰的概率质量大小不同。例如，在基因调控模型 ( Modrák, 2018) 中，一些数据承认两种不同的调控机制具有相反的效应符号，而接近于零的效应具更低的后验密度。这个问题相较前面两者更具挑战性。

处理方法：在某些情况下，可以使用堆叠（预测模型平均 ）作为近似解决方案，但这并不完全通用，因为它需要定义一个预测性的兴趣量；一种更完整的贝叶斯替代方法是引入具有混合组份的强信息性先验将模型分成若干部分，然后在给定先验的每一个组件的情况下分别拟合模型；其他时候，可以使用具有排除某些可能峰效果的强先验来解决问题。

（ 4 ）类型四：具有算术不稳定尾部的高概率质量单一后验形态。

特点：概率质量比较集中的单一峰，但在尾部存在算数不稳定的低矮峰。

处理方法：如果分布的质量附近初始化，则大多数推断应该不会出现问题。如果对极其罕见的事件特别感兴趣，那么应该对问题进行重新参数化，因为从通常默认的几百到几千的有效样本量中，能够学到的东西非常有限。

5.7 通过重新参数化来适应复杂的后验分布¶

通常，基于 HMC 的采样器在如下情况下效果是最好的：质量矩阵得到适当调整，并且联合后验分布的几何形态比较平滑，不存在犄角、尖头等不规则之处。这对于许多经典模型来说很容易满足，例如 Bernstein-von Mises 定理等表明：当有足够数据时，后验会变得相当简单。

不幸的是，一旦模型变得稍微复杂一点，我们就不能再保证有足够数据来实现这个乌托邦了。对于这些模型，通过选择使后验几何更简单的参数化方案，可以极大地改善 HMC 的行为。

例如，当分组级别的方差参数接近零时，分层模型可能会在极限情况下存在困难的漏斗病例（ Neal, 2011 ），但如果遵循 Meng and van Dyk, 2001 讨论的原则，许多此类情况的计算困难问题都可以通过重新参数化来解决，；另见 Betancourt and Girolami, 2015。

5.8 通过对联合概率分布的结构化和边缘化降低计算难度¶

后验分布中具有挑战性的几何形态通常是来自于参数之间的相互作用。一个例子是上面提到的漏斗形状，我们可以在绘制分组级别的尺度参数 \(\phi\) 和个体级别的均值参数 \(\theta\) 的联合概率密度时清楚地观察到。相比之下， \(\phi\) 的边缘概率密度表现会好很多，可以很高效地从边缘后验中抽取 MCMC 样本，

\[ p(\phi \mid y) = \int_\Theta p(\phi,\theta \mid y) \ \ \rm{d}\theta \]

为了利用 MCMC 得到后验图，贝叶斯规则告诉我们只需要边缘似然 \(p(y \mid \phi)\) 。然后，通过从条件分布 \(p(\theta \mid \phi,y)\) 中精确采样来恢复 \(\theta\) 的后验。这种区分参数的边缘化策略特别适用于具有正态似然的高斯过程（例如 Rasmussen and Williams, 2006 ；Betancourt, 2020b ）。

一般而言，我们并不知道 \(p(y \mid \phi)\) 和 \(p(\theta \mid \phi,y)\) 的概率密度，但可以使用拉普拉斯近似等方法来近似这些分布，特别是对于隐高斯模型 （例如，Tierney and Kadane, 1986 ；Rasmussen and Williams, 2006 ；Rue et al., 2009 ； Margossian 等，2020b ）。当与 HMC 结合使用时，如果允许以偏差为代价，这种边缘化方案可以得到比重新参数化更有效的结果；有关该主题的讨论，请参阅 Margossian et al., 2020a。

5.9 增加先验信息¶

通常，计算中的问题可以通过包含已经可用但尚未包含在模型中的先验信息来解决，例如，因为领域专家的先验获取成本很高（ O'Hagan et al., 2006 ； Sarma and Kay, 2020 ），我们通常会从一些模型模板和先验模型开始（第 2.1 节）。在某些情况下，运行更多迭代也有帮助。但是当添加合理的先验信息时，许多拟合问题会消失，但这并不是说先验的主要用途是解决拟合问题。

我们可能已经假设（或希望）数据对模型的所有部分都足够有用，但仔细检查或作为计算诊断的副产品，我们会发现情况并非如此。在经典统计中，模型有时被归类为可识别的或不可识别的，但这会产生误导（即便在添加了部分或弱可识别性等中间类别之后），因为可以从观测中学到的信息量同时还取决于实际获得的数据。此外，“识别”在统计学中被正式定义为渐近属性，但在贝叶斯推断中，我们关心对有限数据的推断，特别是考虑到模型通常会随着数据增加而增大和复杂。渐近结果可以提供对有限样本性能的一些见解，但我们通常更愿意考虑后验分布。 Lindley , 1956 以及 Goel and DeGroot , 1981 讨论了如何测量实验提供的信息，即量化后验与先验之间的不同。如果数据没有提供模型某些方面的信息，那么我们可以通过先验来提供更多信息以改善情况。

此外，我们通常更喜欢一个参数能够被数据更新的模型，而不是一个更接近事实但数据无法提供足够信息的模型。在第 6.2-6.3 节中，我们将讨论用于评估数据点或超参数信息量的工具。

我们用估计一个 “衰减率未知的下降指数求和” 任务来做说明。这项任务是数值分析中众所周知的病态问题之一，也出现在药理学等应用中 ( Jacquez, 1972 )。

假设数据

\[ y_i = (a_1e^{−b_1x_i} + a_2e^{−b_2x_i} ) \times e^{\epsilon_i} ，\mathrm{for} \ i = 1,...,n, \]

其中，独立误差 \(\epsilon_i \sim \text{Normal}(0,\sigma)\) ；系数 \(a_1\)、\(a_2\) 和残差标准差 \(\sigma\) 被限制为正数；参数 \(b_1\) 和 \(b_2\) 也是正的（假设呈指数下降而不是指数增长），并且被限制为有序，即 \(b_1 < b_2\) ，以便唯一地定义两个模型组件（注解：这里可能存在可识别性问题，因此做了强制的排序）。

图 12：来自模型 \(y = (a_1e^{−b_1x} + a_2e^{−b_2x}) ∗\mathrm{error}\) 的模拟数据：(a) 设置 \((b_1,b_2) = (0.1,2.0)\) ，和 (b) 设置 \((b_1,b_2) = (0.1,0.2)\) 。右图中绘制的虚线代表曲线 \(y = 1.8e^{-0.135x}\)，它几乎与这些数据范围内的真实曲线完全一致。很容易将模型拟合到左图中的数据并恢复真实的参数值。然而，对于右图中的数据，观测并没有提供充分减少不确定性的信息，模型只能使用提供了必要缺失信息的先验分布来稳定地拟合。

我们从能够清楚区两条曲线的模型设置开始模拟伪数据。

第一条曲线采用设置 \(b_1 = 0.1\) 和 \(b_2 = 2.0\)，两个下降指数的尺度通过比例因子 \(20\) 被区分开。我们模拟了 \(1000\) 个数据点，其中预测变量 \(x\) 在 \(0\) 到 \(10\) 之间均匀分布，并且随意地设置 \(a_1 = 1.0, a_2 = 0.8, \sigma = 0.2\) 。图 12a 显示了真实曲线和相应的模拟数据。然后我们使用 Stan 从数据中拟合模型。模拟运行平稳，后验推断恢复了模型的五个参数，鉴于数据是从模型中模拟出来的，这并不奇怪。

第二条曲线使问题变得稍微困难了一些。该模型在构造上仍然是正确的，但是我们没有将 \((b_1,b_2)\) 设置为 \((0.1,2.0)\)，而是将其设为 \((0.1,0.2)\)，现在只有 \(2\) 的比例因子将两个下降指数的尺度分开。模拟数据如图 12b 所示。现在尝试在 Stan 中拟合模型时，会得到可怕的收敛。两个下降指数变得非常难以区分，图中还包含了一条曲线为 \(y = 1.8e^{−0.135x}\) ，你会发现该曲线基本上不可能与真实模型区分开来。

遇到这种情况时，我们可以使通过添加先验信息使计算更加稳定。例如，将所有参数的先验都默认为 \(\text{Normal}(0,1)\) ，就可以做充分正则化。但如果模型已经建立，且参数大致在单位尺度上，那么先验信息仍然很弱，如第 2.3 节所述。此时还可以通过与 \(\text{Normal}(0,0.5)\) 和 \(\text{Normal}(0,2)\) 等替代方案进行比较，来检查后验推断对先验选择的敏感性；或者通过相对于先验尺度做推断区分，如第 6.3 节中的讨论。

使用信息性正态分布先验会增加后验密度的对数凹尾（ tail-log-concavity ），进而有利于加快 MCMC 的混合时间。但是信息性先验并不代表这用模型准确度来换取计算效率；相反，在这种情况下，随着计算负担减轻，模型拟合的效果也会得到改进，这是第 5.1 节中讨论的通俗定理的一个实例。更一般地，也可以有厚尾的强先验，但此时后验密度无法保证对数凹尾。另一方面，先验在似然的上下文中也可能很弱，但同时仍然能保证对数凹尾。这与 “模型的信息量取决于提出的问题” 有关。

在抽象层次上，可以阶梯性地考虑以下四个步骤：

（ 1 ）MCMC 混合不良；
（ 2 ）困难的几何形态作为对上述的数学解释；
（ 3 ）模型某些部分的弱信息数据作为对上述的统计解释；
（ 4 ）实质性先验信息作为上述问题的解决方案。

从该阶梯的起点开始，我们有计算的故障排除；从该阶梯的终点开始，我们有计算的工作流程。

作为另一个例子，当分组级别的方差参数接近零，而我们试图在极限范围内避免分层模型的漏斗病态时，可以使用避零先验（例如，对数正态分布或逆伽马分布 ）以避开似然的高曲率区域； Chung et al., 2013, 2014 讨论了惩罚边缘似然估计的相关想法。当此类先验信息可用时，避零先验是有意义的，例如对于高斯过程的长度尺度参数（ Fuglstad et al., 2019 ）。不过我们在此使用此类约束，只是为了使算法运行速度更快。最终要理解，如果使用约束性先验来加速计算，其实质是在向模型添加信息。

更一般地说，我们发现，统计拟合算法的不良混合通常可以通过更强的正则化来解决。然而，这并不是免费的：为了在不模糊待估计模型各要素的情况下实施有效的正则化，我们需要一些主题相关的知识（即真实的先验信息、专家只是、经验信息等）。随意调整模型直到计算问题消失非常危险，可能会威胁到推断的有效性，不过现在还没有一种好方法来诊断这个问题。

5.10 增加数据¶

与添加先验信息类似，可以通过添加在模型内处理的新数据源来约束模型。例如，一次校准实验可以得到结果变量的标准差信息。

在其他情况下，对于较大数据集表现良好的模型在小数据范围内可能存在计算问题；图 13 展示了一个例子。虽然在这种情况下，后验的漏斗形态看起来与分层模型中的漏斗相似，但这种病态更难避免，通常只能认为：数据未能提供支撑完整模型的信息，需要使用更简单的模型。 Betancourt , 2018 进一步讨论了这个问题。

图 13：通过添加数据改变后验分布几何形态的示例。(A) 正态似然作为平均值 \(\mu\) 和标准差 \(\sigma\) 的函数，随着从标准正态分布中抽取的可用观测值数量 \(N\) 增加而不断增加。对于 \(N = 1\)，似然无界的增长，对于 \(N = 2\)，几何形态仍然是漏斗状，会导致计算问题。对于 \(N = 8\) ，漏斗状形态被有效抑制。 (B) 使用 Lotka-Volterra 人口动态模型（总共 \(8\) 个参数）的实际例子，如 Carpenter , 2018 所述，图中显示了在 Stan 中拟合模型得出的两个参数的后验样本散点图。具有 \(6\) 个数据点的拟合显示出漏斗形状。红点表示发散的转移，表明采样器遇到了困难的几何形态，其结果不可信。当使用 \(9\) 个数据点时，Stan 能够拟合模型，尽管几何形态略有不均匀。当拟合到 \(21\) 个数据点时，该模型表现良好。

6 评估和使用拟合后的模型¶

模型拟合好后，评估拟合效果的工作流程会更加复杂，因为可以检查许多不同的东西，而这些检查中的每一个都可以指向许多方向。统计模型可以针对多个目标进行拟合，并针对不同用户群体开发相应的统计方法。模型需要被检查的要素主要取决于应用。

6.1 后验预测性检查¶

后验预测性检查类似于先验预测性检查（第 2.4 节），只不过使用的参数实例抽取自后验分布而不是先验分布。先验预测性检查是一种理解模型和指定先验含义的方法，而后验预测性检查允许人们检查模型对真实数据的拟合情况 （ Box, 1980 ， Rubin, 1984 ， Gelman、Meng and Stern, 1996 ）。

将来自后验预测性分布的模拟数据集与实际数据集进行比较时，如果实际数据集不能代表后验预测性分布，则表明模型无法描述数据的某些方面。最直接的检查就是“对比来自预测性分布的模拟数据与来自完整分布的真实数据”，或 “比较两者之间在数据集（或数据子集）上计算的汇总统计量，特别是对于模型中未包含的分组” （参见图 14 ）。

图 14：利用 bayesplot R 包实现的一个通过后验预测性检查诊断出的模型未拟合示例。在所有图中，\(y\) 是基于输入数据绘制的，而 \(y_{rep}\) 是基于后验预测性分布绘制的。 (A) 概率密度检查的例子。一个正态分布的后验预测性检查，用于判断模型是否拟合了从对数正态分布模拟生成的数据：\(y_{rep}\) 的尾部表现与 \(y\) 非常不同。 (B) 统计量检查的例子。一个面向二项分布的后验预测性检查，用于判断模型是否拟合了从 \(\beta-\) 二项分布模拟生成的数据。将 \(y_{rep}\) 数据集的标准差 (sd) 直方图与 \(y\) 的标准差进行比较。检查显示数据具有比模型更大的标准差。 (C) 概率质量检查的例子。离散数据的 “直方图” 检查（注意 \(y\) 和 \(y_{rep}\) 之间的颜色切换）。这个检查看起来不错：\(y\) 中个体计数的频率分布完全落在 \(y_{rep}\) 中。 (D) 相同的模型和数据，但检查按预测变量进行了分组，该预测变量未包含在模型中但实际上影响了结果变量。其中分组 High 比较系统地偏离了 \(y_{rep}\) 的范围，表明存在模型尚未能捕获的额外变异源。

没有通用方法来选择应该对模型执行哪些检查，但运行上述的直接检查可以防止比较明显的错误指定。同时，也没有通用方法来判定这种检查何时会失败，需要调整模型。根据分析目标以及特定情况的成本和效益，我们可能会容忍模型未能捕获数据的某些方面，或者必须着力于改进模型。总的来说，我们试图找到 “严格的测试” （ Mayo, 2018 ）：如果模型对我们最关心的问题给出误导性的答案，那么这些检查可能会失败。

图 15 显示了来自应用的更复杂的后验预测性检查项目。这个例子展示了预测性模拟与图形显示相结合的方式，它也说明了预测性检查的实际挑战，因为通常需要针对特定问题提出独特的可视化。

图 15：后验预测性检查示例。左列显示从心理学实验中得到的观测数据，显示为 \(15×23\) 的二值响应数组，来自 \(6\) 个参与者中的每一个，按行、列和人排序。右边的列显示了来自拟合模型的后验预测性分布的 \(7\) 个复制数据集。每个复制的数据集都已完成排序，因为这是显示的一部分。检查揭示了观测数据中有而复制数据中没有的一些模式，表明模型与数据之间在某个方面未能拟合。引自 Gelman et al., 2013。

6.2 单数据点和数据子集的交叉验证和影响¶

后验预测性检查通常足以揭示模型的错误拟合，但由于它既使用数据做模型拟合，又将其用于评估错误拟合，因此结果可能过于乐观。

在交叉验证中，部分数据点会被保留，模型仅拟合剩余的数据点，并在保留的数据上检查预测性能。这改进了预测性检查诊断，对于灵活模型尤其如此（例如，参数多于观测值的超参数化模型）。我们发现对进一步评估模型有用的三种交叉验证诊断方法是：

（ 1 ）使用交叉验证预测新分布做校准检查；
（ 2 ）识别哪些观测或观测分组更难于预测；
（ 3 ）识别特定观测的影响程度，即在其他观测之外，它提供了多少信息。

在所有三种情况下，使用重要性采样对留一法交叉验证做高效估计，通过在每个数据点被保留时消除再拟合模型的需要来促进实际使用（ Vehtari et al., 2017 ； Paananen et al., 2020 ）。

尽管预测性分布的完美校准不是贝叶斯推断的最终目标，但查看留一法交叉验证 (LOO-CV) 预测性分布的校准程度，可以发现改进模型的一些机会。虽然后验预测性检查一般都是将预测的边缘分布与数据分布进行比较，但留一法交叉验证预测性检查着眼于条件预测性分布的校准。在良好校准下，给定保留观测值情况下，预测性分布的条件累积概率分布（也称为概率积分变换，PIT ）应当是均匀的。对均匀性的偏离可以显示出预测性分布的不足或过度分散。图 16 显示了 Gabry et al., 2019 的一个例子。其中留一法交叉验证的概率积分变换值 (LOO-PIT) 在中间附近过于集中，表明与实际条件观测值相比，预测性分布过于分散。

图 16：留一法交叉验证概率积分变换 (LOO-PIT) 图评估拟合模型的预测校准。在完美校准下，LOO-PIT 值 将是均匀的。当前情况下的值集中在中间附近，表明预测性分布过于分散 ( Gabry et al., 2019 )。

除了查看条件预测性分布的校准之外，还可以查看哪些观测值更加难以预测，并查看是否存在模式或解释说明为什么有些观测值比其他观测值更难预测。这种方法可以揭示数据或数据处理中的潜在问题，或者指出模型改进的方向（ Vehtari et al., 2017 ； Gabry et al., 2019 ）。我们通过对孟加拉国一个地区居民的调查进行分析来说明，该地区受到饮用水中砷的影响。井中砷含量升高的受访者被问及他们是否有兴趣从邻居的井中换水，并且一系列模型适用于在给定家庭信息的情况下预测这种二元响应（ Vehtari et al., 2017 ）。

图 17 比较两个模型的逐点对数评分。图 17a 左侧和图 17b 右下方的分散蓝点对应于其中一个模型拟合特别差的数据点，即对预期对数预测密度的负贡献较大。我们可以对所有逐点差异求和，以产生预期对数预测密度 \(\rm elpd_{loo}\) 的估计差异 \(16.4\) ，标准差仅为 \(4.4\) ，但除此之外，我们可以使用此图来找出哪些数据点给模型带来了问题，在此案例中，是 \(10-15\) 个现有砷水平非常高的非交换者。

在这个特定的例子中，我们没有跟进建模问题，因为即使是更适合数据的精细模型也不会改变结论，因此不会改变孟加拉国的任何建议行动。Gabry et al., 2019 提供了一个例子，其中 LOO-CV 表明问题促使人们努力改进统计模型。

图 17：逻辑回归示例，依据逐点对留一法交叉验证误差的贡献来比较两个模型：（ a ）直接比较 LOO 的贡献； (b) 将 LOO 差异作为关键预测因子（现有砷含量）的函数。为了帮助理解，我们根据二值输出对数据进行了着色，红色对应于 \(y = 1\) ，蓝色代表 \(y = 0\) 。对于任何给定数据点，一个模型将比另一个拟合的更好，但对于本例，图表显示，两个模型之间的 LOO 差异是由于线性模型对 \(10-15\) 个特定数据点的预测不佳造成的。引自 Vehtari et al., 2017 。

上述两种方法侧重于预测，但我们也可以看看当每个数据点被遗漏时参数推断会发生什么变化，这提供了对每次观测所带来影响的感觉。与更普遍的交叉验证一样，如果数据被聚类或以其他方式结构化，因此需要删除多个点才能产生效果，这种方法有局限性，但仍然可以作为贝叶斯工作流程的一部分，因为它的计算成本很低并且在许多应用环境中都很有价值。遵循这种交叉验证的想法，可以根据逼近 LOO-CV 时计算的重要性权重分布的性质来总结单个数据点 \(y_i\) 的影响（有关逼近和相应实现的详细信息，请参见 Vehtari et al., 2017 ，以及 loo R 包 (Vehtari et al., 2020)。

重要性加权的另一种方法是将数据点的删除作为更大模型空间中的梯度。假设我们有一个简单的独立似然，3\(\prod^n_{i=1} p(y_i \mid \theta)\) ，我们使用更一般的形式，\(\prod^n_{i=1} p(y_i \mid \theta)^{\alpha_i}\) ，当对于所有 \(i\) ，\(\alpha_i = 1\) 。留一法交叉验证对应于一次为一个观测设置 \(\alpha_i = 0\) 。但是另一种选择，由 Giordano et al. , 2018讨论并由 Giordano , 2018 实施，是计算作为 \(\alpha\) 函数的增广对数似然的梯度：这可以解释为一种差分交叉验证或影响函数。

多级（分层）模型的交叉验证需要更多的思考。留一法仍然是可能的，但它并不总是符合我们的推断目标。例如，在执行多级回归以调整政治调查时，我们通常对估计州级的意见感兴趣。模型可以在州级显示真正的改进，而这在个人观测的交叉验证级是无法检测到的（ Wang and Gelman, 2016 ）。 Millar , 2018、Merkle、Furr and Rabe-Hesketh, 2019 and Vehtari , 2019 在分层模型中展示了不同的交叉验证变体及其近似值，包括留一单元输出和留一组输出。在应用问题中，我们进行了混合，保留一些个体观测和一些群体，然后评估两个级别的预测（ Price et al., 1996 ）。

不幸的是，使用重要性抽样来近似这种交叉验证程序往往比在留一法。这是因为一次忽略了更多的观测，这意味着从完整模型到子集模型的后验分布发生了更强的变化。因此，我们可能不得不依赖更昂贵的模型改装来获得留一单元和留一组出交叉验证结果。

6.3 先验信息的影响¶

复杂模型可能难以理解，因此需要探索性模型分析（ Unwin et al., 2003 ； Wickham, 2006 ）和可解释的人工智能（ Chen et al., 2018 ； Gunning, 2017 ； Rudin, 2018 ），它们对方法进行了补充用于使用一系列相互关联的方法评估、比较和平均模型，包括交叉验证、堆叠、提升和贝叶斯评估。

在本节中，我们将讨论了解拟合模型下的后验推断如何依赖于数据和先验的方法。

统计模型可以通过两种方式来理解：生成式和推断式。从生成的角度来看，我们想了解参数如何映射到数据。我们可以执行先验预测性模拟来可视化模型中可能的数据（如图 4 所示）。从推断的角度来看，我们想了解从输入（数据和先验分布）到输出（估计和不确定性）的路径。

了解先验影响的最直接方法是通过用多个先验重新拟合模型来运行敏感性分析，然而，如果模型需要很长时间才能拟合，这可能会非常昂贵。但是，有一些捷径。

图 18：针对毒理学问题的统计模型的静态敏感性分析示例。每个图都显示了在两种条件下代谢百分比的后验模拟图（因此每个图顶部和底部的点簇），绘制与模型中的两个参数的关系。这些图显示，在一种条件下对代谢百分比的推断对任一参数几乎没有敏感性，但对另一种情况具有很强的敏感性。这种图表可用于评估对替代先验分布的敏感性，而无需重新拟合模型。来自 Gelman、 Bois and Jiang, 1996。

一种方法是计算先验和后验之间的收缩，例如，通过比较每个参数的先验标准差或通过比较先验和后验分位数。如果相对于先验的数据对特定参数提供信息，则该参数的收缩应该更强。这种类型的检查在文献中得到了广泛的发展。例如，参见 Nott et al., 2020 。

另一种方法是使用重要性采样来使用旧模型的后验来近似新模型的后验，前提是两个后验足够相似以便重要性采样桥接（ Vehtari et al., 2019 ； Paananen et al.。 , 2020)。如果不是，这本身也是有价值的信息（参见第 6.2 节）。

另一种方法是执行静态敏感性分析，这是一种研究后验推断对先验扰动的敏感性的方法，而不需要模型使用替代先验分布重新拟合（ Gelman、Bois and Jiang, 1996 ；示例见图 18 ）。图 18 中的每个图都显示了后验模拟，揭示了作为模型中参数的函数的感兴趣量（在本例中，是身体代谢的毒素的百分比）的后验依赖性。

将图 18 视为四个散点图，因为这两个图表中的每一个实际上都是叠加的两个图，一个用于低暴露于毒素的条件，另一个用于高暴露。这四个图的每一个都可以用两种方式解释。首先，直接解释显示了感兴趣的预测量（体内代谢的百分比）与特定参数（例如，毒素代谢的比例系数）之间的后验相关性。其次，可以间接读取每个散点图，以揭示绘制在 \(y\) 轴上的数量对绘制在 \(x\) 轴上的参数的先验的敏感性。解释如下：绘制在 \(x\) 轴上的参数的 35 先验分布的变化可以通过根据新的先验与分析中使用的比率重新加权图表上的点来执行。使用这些图，重要性权重可以隐式地可视化：可以根据散点图中的依赖性看到先验分布变化的影响。

也可以更正式地研究先验和数据到后验的映射，如第 1 节所述 6.2.

6.4 对推断结果进行总结¶

贝叶斯推断非常适用于具有隐变量和其他不确定性无法解决的设置的问题。此外，我们经常使用分层模型，其中包括表示变化的成批参数。例如，在报告我们的选举预测模型的结果时，我们对预测投票的不确定性以及各州之间的变化感兴趣。

不幸的是，显示贝叶斯推断的常用方法并没有完全捕捉到我们的多层次变化和不确定性。推论。参数估计、不确定性和标准差的表格甚至图表仅显示一维边缘，而边缘后验分布图对于具有许多参数的模型来说是笨拙的，并且也无法捕捉层次结构中不确定性和变化之间的相互作用模型。

首先，我们应该遵循良好的统计实践和图形数据和拟合模型的一般原则，既为了“探索性数据分析”的目的，即发现数据中的意外模式，也更直接地了解模型与模型的关系。用于拟合它的数据。

我们从 Ghitza and Gelman , 2020 对 2016 年美国总统大选投票偏好的分析中说明了图形模型探索和总结的一些用途。图 19 显示了支持两位候选人的估计性别差距，首先显示为地图，然后显示为散点图。该地图是可视化这些估计的自然第一步，它揭示了一些有趣的地理模式，然后在散点图中以更集中的方式进行探索。图 20 通过与更简单的县级估计进行比较来评估模型。此示例演示了探索性图形中的一般工作流程，其中推断总结的结果激发了未来的探索。

Gabry et al., 2019 提出了我们关于贝叶斯工作流程图形的一些想法，其中一些已经在 R 包 bayesplot 中实现（ Gabry et al.,2020b ，另见 Kay,2020ab ； Kumar,2019 ）。概率编程最终有可能允许像任何其他数据对象一样操纵随机变量，所有绘图和计算中都隐含着不确定性（ Kerman and Gelman, 2004, 2007 ），但需要做更多的工作来将这种似然转化为现实，超越点估计和区间估计，以便我们可以充分利用我们拟合的模型。

7 修改模型¶

模型构建是一项类似于语言的任务，其中建模者将现有组件（线性、逻辑和指数函数；加法和乘法模型；二项式、泊松和正态分布；变系数等）组合在一起，以包含新的数据和现有数据的附加特征，以及与感兴趣的底层过程的链接。

如第 2.2 节所述，模型的大部分部分可以被视为允许替换或扩展的占位符。或者，我们可能会发现需要使用近似模型或近似算法，如第 3.3 节所述。

图 19：从模型拟合到 \(2016\) 年美国总统竞选期间的调查数据：(a) 估计之间的性别差距白人选民，(b) 绘制的估计性别差距与奥巴马 \(2012\) 年白人选票中估计的县级选票份额。每个圆圈的面积与该县的选民人数成正比，地图上的颜色代表从深紫色（比较白人男性和白人女性的投票模式没有差异）到浅灰色（支持克林顿的白人女性）的范围比白人男性多 \(7.5\) 个百分点）到深绿色（白人性别差距为 \(15\) 个百分点）。显著的地理模式——南部大部分地区的白人性别差距较低，西部以及东北部和中西部的大部分地区较高——激发了散点图，这表明白人性别差距往往在白人人口较多的县中最高。投票接近平分。来自 Ghitza and Gelman , 2020 。

图 20：对图 19 所示模型的一个方面的评估，比较了根据两个不同模型估计的县级对克林顿的支持。我们将其包含在这里是为了说明可以使用图形方法来总结、理解、评估和比较拟合模型的方式。来自 Ghitza and Gelman , 2020 。

模型扩展能够响应新数据、模型对已有数据拟合失败、现有拟合过程计算困难等一系列问题。对于 Gelman、Hullman et al., 2020 中描述的选举预测。我们从 Linzer, 2013 的民意调查聚合模型开始，我们在 \(2016\) 年对其进行了调整，但在某些摇摆州的预测明显失败，我们将其归因于州之间投票摇摆的相关性建模不佳，以及民意调查中存在非抽样错误（ Gelman and Azari, 2017 ）。在我们的修订中，扩展了模型以包含这两个功能。第 10 节和第 11 节给出了迭代模型构建和评估的扩展示例。

7.1 为数据构建模型¶

在统计学的教科书中，给定参数的数据分布通常只是给定的。但在应用中，我们希望基于对数据拟合情况和领域专业知识的组合，来建立数据模型。如果模型是从一个有限菜单中做选择，那我们至少希望它是开放的。通常，数据模型最重要的方面不是其概率分布形式，而是数据如何与感兴趣的基础参数相关联。例如，在选举预测中，民意调查模型包括个人民意调查和民意调查平均值的非抽样误差项（ Shirani-Mehr et al., 2018 ）。

此时一个相关挑战出现了，因为数据在到达我们之前通常会经过预处理，因此任何生成式模型都必然是一个近似值。这可能出现在元分析中、或者使用机器学习算法和降维技术将大量预测变量组合成一个或两个汇总变量的场景中。与往常一样，我们需要关注数据质量，数据模型最重要的方面可能是偏差，而不是传统上被认为是测量误差的那个部分。理解这一点会影响贝叶斯工作流程，因为扩展模型以允许系统误差非常有价值，我们将在 10.5 节给出一个例子。

7.2 合并额外的数据¶

有时统计方法最重要的方面不是它对数据做了什么，而是使用了什么数据。贝叶斯工作流程的一个关键部分是扩展模型以利用更多数据。这可以像添加回归预测变量一样简单，但当添加更多参数时，可能有必要假设并非所有参数都能同时对模型产生重大影响。看到这一点的一种方法是将参数的添加视为先前集中在零的先验分布的松弛。例如，通过在回归中添加交互项来扩展选举模型以解释政治两极分化，从而允许近年来国民经济预测因子的系数较低。有时我们有两种类似数据的测量形式，因此需要为两个数据源创建一个生成模型。

有时这会带来技术挑战，例如当将样本的直接测量与总体汇总统计相结合，或整合不同质量的测量时（例如， Lin et al., 1999 ），或当信息可用于表格的部分边缘时（ Deming and Stephan, 1940 ）。在 Weber et al., 2018 中，我们拟合了一个药理学模型，该模型具有一组服用某药物的患者的直接数据，但只有一组接受过竞争对手产品的患者的平均数据。为了避免将所有平均患者的结果建模为潜在数据的计算成本，我们设计了一种分析方法来近似相关积分，以便可以将平均数据包含在似然函数中。

7.3 使用先验分布¶

在贝叶斯统计中，我们经常谈论非信息性或完全信息性先验，但这两种先验通常都不存在：均匀先验包含了一些信息，因为它取决于模型的参数化方案；参考先验取决于用于收集新数据的假设渐近机制（ Berger et al., 2009 ）；一个信息丰富的先验也很少能够包含所有的可用知识。

相反，我们更愿意考虑似然的阶梯：

-> （不适当的）平坦先验；

—-> 模糊但适当的先验；

——-> 信息非常弱的先验；

———-> 通用的弱信息先验；

————-> 特定的信息性先验。

例如，我们的选举模型包括随机游走条款，以允许州和国家层面的态度变化。这些随机游走中的每一个都有标准差参数，对应于一天内无法解释的变化（在 \(\text{logit}\) 尺度上）。

这个新创的标准差可以被赋予：

一个均匀先验
一个模糊但适当的先验，例如，\(\text{Normal}^+(0 ,1000)\)，其中使用符号 \(\text{Normal}^+\) 来表示被截断为正的正态分布
一个弱先验，例如，百分比尺度上的 \(\text{Normal}^+(0,1)\) ，这将允许不切实际的大日常支持候选人的概率会发生 \(0.25\) 个百分点的变化，但仍会使分布远离极端参数值
或者更多信息的先验，例如 \(\text{Normal}^+(0,0.1)\) ，它不包含我们所有的先验知识，但确实将此参数软限制在合理范围内。

我们的观点是：选择先验的同时，也在选择要包含在分析任务中的主题相关信息的量。

理解先验分布的另一种视角是将其作为约束。例如，如果我们拟合线性模型加上样条或高斯过程，\(y = b_0 + b_1x+ g(x) + \text{error}\) ，其中非线性函数 \(g\) 有界，然后在 \(g\) 上具有足够强的先验，我们正在拟合一条接近线性的曲线。在此示例中，先验分布可以表示先验信息，也可以被其视为近似线性曲线模型的声明的一部分。Simpson et al., 2017 进一步讨论了使用先验分布向更简单的模型收缩。这也导致了更普遍的观点，即先验就像统计模型的任何其他部分一样，可以用于不同的目的。模型和先验之间的任何明确区别在很大程度上是任意的，并且通常主要取决于创建者的概念背景区别。

获得合理推断所需的先验信息量很大程度上取决于参数在模型中的作用以及参数在层次结构中的深度（ Goel and DeGroot, 1981 ）。例如，主要控制中心量（例如均值或中位数）的参数通常比尺度参数更能容忍模糊先验，这比控制尾部量的参数更能容忍模糊先验，例如广义极值分布的形状参数。当模型具有层次结构时，与数据更接近的参数通常比层次结构更下方的参数更愿意容忍模糊的先验。

在贝叶斯工作流程中，一系列模型需要先验。通常，随着模型变得更加复杂，所有先验都需要变得更紧密。以下多级数据的简单示例（例如，参见 Raudenbush and Bryk, 2002 ）说明了为什么这可能是必要的。

考虑数据为 \(y_{ij}, i = 1,...,n_j, j = 1, ..,J\) 的工作流程。这里 \(i\) 索引了观测，\(j\) 索引了组。想象一下，对于工作流程中的第一个模型，我们选择忽略组结构并使用简单的正态分布来表示与均值的偏差。在这种情况下，一些信息预算将用于估计总体均值，而其中一些是花在观测标准差上。如果我们有适量的数据，无论先验有多弱，均值将近似为 \(\bar y= \sum^n_{i=1} y_i/n\) 。此外，新观测的预测性分布将近似为正态 \((\bar y,s)\)，其中 \(s\) 是样本标准差。同样，这对于观测标准差的大多数合理先验都是正确的，无论它们有多模糊。

如果工作流程的下一步是使用多级模型以允许平均值随组变化，那么信息预算仍然需要在标准差和均值之间进行划分。但是，该模型现在有 \(J + 1\) 个额外参数（每个组一个，组间标准差一个），因此需要进一步划分均值的预算以对组均值建模，而标准差的预算需要在组内偏差和组间偏差之间分配。但是我们仍然拥有相同数量的数据，因此需要小心确保这种模型扩展不会破坏估计。这意味着除了在新参数上放置适当的先验之外，可能还需要收紧总体均值和观测标准差的先验，以免信息缺乏导致无意义的估计。

一个相关的问题是在高维空间中测量的汇集。例如，在预测变量的数量增加的回归中，如果希望将大部分先验质量保持在峰附近（因为远离峰的体积增加得更快），则系数向量的先验需要在峰附近具有更高的密度，参见 Piironen and Vehtari, 2017 。考虑线性或逻辑回归，以及如果权重的边缘先验是固定的， \(R^2\) 上的先验会发生什么。如果我们希望 \(R^2\) 上的先验保持不变，则权重的先验需要变得更紧。图 21 使用先验预测性检查（参见第 2.4 节）来显示线性回归中 \(26\) 个权重的通常弱先验如何对应于强烈支持更高 \(R^2\) 值的先验，也影响后验。从这个角度来看，如果它们大量出现，弱信息但独立的先验可能共同具有强信息。

必须为工作流程中的每个模型指定先验。扩展模型可能需要额外考虑参数化。例如，当从 \(\text{Normal}( \mu ,\sigma)\) 到具有 \(ν\) 自由度的 \(t\) 分布 \(t_ν( \mu ,\sigma)\) 时，我们必须注意 \(\sigma\) 上的先验。尺度参数 \(\sigma\) 对于两者看起来相同模型，但 \(t\) 分布的方差实际上是 \(\frac{ν}{ν−2} \sigma^2\) 而不是 \(\sigma^2\) 。因此，如果 \(ν\) 很小，则 \(\sigma\) 不再接近残差标准差。

一般来说，我们需要考虑模型中所有参数的联合先验，这个可以在生成式模型的上下文中评估，以免不幸的取消或共振导致比建模者实际想要的更不稳定或信息量更多的先验（ Gelman et al., 2017 ； Kennedy et al., 2019 ）。正如第 2.4 节所讨论的，先验预测性检查是在特定数据模型的上下文中探索和理解先验分布的一种很好的通用方法。

上述示例包含关于先验的特定信息，但也包含关于我们如何看待工作流程的元信息在构建统计模型时。诸如“信息预算仍然需要划分”之类的短语代表了我们做出的关于我们为包含先验信息付出了多少努力的重要但有些非正式的决定。在不承认已经做出的权衡和选择的情况下，在按原样呈现最终模型的文章或教科书中，这种担忧并不总是很清楚。

图 21：贝叶斯 \(R^2\) 的先验和后验分布用于回归预测学生成绩的 \(26\) 个预测变量，使用三个不同的先验系数：(a) 默认弱先验，(b) 正常先验与预测变量的数量成比例，以及（ c) 规范化的马蹄形先验。来自 Gelman、Hill and Vehtari , 2020 的第 12.7 节。

7.4 模型拓扑¶

考虑一类模型，为了简单起见在某些特定的受限域中，例如自回归移动平均 (ARMA) 模型、二元分类树或具有一些固定输入变量集的线性回归。任何这些框架中的模型都可以构造为偏序：例如，\(\text{AR(1)}\) 比 \(\text{AR(2)}\) 简单，\(\text{AR(2)}\) 比 \(\text{ARMA(2,1)}\) 简单， \(\text{MA(1)}\) 也比 \(\text{ARMA(2,1)}\) 简单，但 \(\text{AR(1)}\) 和 \(\text{MA(1)}\) 本身不是有序的。类似地，树分裂形成偏序，线性回归中包含或排除的 \(2^k\) 种似然可以构造为立方体的角。正如这些例子所示，这些模型框架中的每一个都有自己的拓扑或网络结构，由类中的模型及其偏序决定。

我们称之为模型的拓扑而不是概率空间，因为我们不一定感兴趣为各个模型分配概率。我们在这里的兴趣不是对模型求平均，而是在模型之间导航，拓扑是指模型之间以及网络中相邻模型中参数之间的连接。

这个想法的一个示例实现是 Automatic Statistician ( Hwang et al., 2016 ； Gharamani et al., 2019 )，它在特定但开放的类别（例如，时间序列模型和线性回归模型）中搜索模型，使用推断和模型批评来探索模型和数据空间。我们相信，通过对统计建模语言诱导的模型拓扑结构的更正式理解，可以更好地理解并最终改进此类算法。另一个方向是基于菜单的软件包，例如 Prophet （ Taylor and Lethem, 2018 ），它允许用户从一组构建块中组合模型（在这种情况下，用于时间序列预测）。在此类包中，重要的是不仅要能够构建和拟合这些模型，而且要了解每个模型与适合相同数据的更简单或更复杂的变体相比。

自动加法模型就足够了，这里每个模型本身都是一个高维对象。不同模型的输出，作为概率随机变量，可以相加、相乘、线性混合、对数线性混合、逐点混合等，这在我们需要指定的模型拓扑选择范围内。

此外，每个模型在一个框架内有它自己的内部结构，涉及可以从数据中估计的参数。而且，重要的是，网络中不同模型中的参数可以“相互交谈”，因为在模型本身的范围之外具有共享的、可观测的意义。在机器学习和应用统计学中，预测和因果推断是两个熟悉的模型共享推断量示例。在预测中，一组越来越复杂的程序可用于特定的预测目标。在因果推断中，可以使用一系列回归来估计治疗效果，从简单的差异开始，然后转向复杂的交互模型，调整观测到的治疗组和对照组之间的差异。回想一下，因果推断是涉及反事实的预测的一种特殊情况；例如，参见 Morgan and Winship , 2014 。

因此，统计或机器学习模型的拓扑结构包括模型的偏序，以及更大框架内跨不同模型的参数或参数函数之间的连接。另一个扭曲是先验分布为结构添加了一个连续的维度，在模型之间架起了桥梁。

8. 理解和比较多个模型¶

8.1 可视化彼此相关的模型¶

贝叶斯工作流程的关键方面超越了贝叶斯数据分析，是我们在处理单个问题的同时拟合多个模型。我们在这里谈论的不是模型选择或模型平均，而是使用一系列拟合模型来更好地理解每个模型。用 Wickham、Cook and Hofmann , 2015 的话来说，我们寻求“探索模型拟合的过程，而不仅仅是最终结果”。我们拟合多个模型有几个原因，包括：

拟合和理解一个大模型可能很困难，所以我们将从简单的模型构建它。
在构建我们的模型时，我们犯了很多错误：拼写错误、编码错误、概念错误（例如，没有意识到观测值不包含模型某些部分的有用信息）等。
随着我们获得更多数据，我们通常会相应地扩展我们的模型。例如，如果我们在做药理学研究并获得一组新患者的数据，我们可能会让某些参数因组而异。
通常我们拟合一个在数学上明确指定的模型，但是一旦我们将其拟合到数据中，我们就会意识到我们可以做更多的事情，所以我们扩展它。
与此相关的是，当我们第一次拟合模型时，我们经常将它与各种占位符放在一起。我们经常从弱先验开始并增强它们，或者从强先验开始并放松它们。
我们将检查模型，发现问题，然后扩展或替换它。这是“贝叶斯数据分析”的一部分；额外的“工作流程”部分是我们仍然保留旧模型，不是为了求平均值，而是为了理解我们在做什么。
有时我们拟合简单模型作为比较。例如，如果你要进行因果推断的大回归，你还需要进行简单的未调整比较，然后查看调整后的结果。
上述想法被列为出于统计考虑，但有时由于计算问题，我们被迫采取行动。

图 22：比较多个模型中感兴趣的数量的推断的假设图。这里的目标不是执行模型选择或模型平均，而是了解当我们从简单的比较（在图的左侧）到最终模型（在图的右侧）时，对感兴趣量的推断如何变化图），经过各种中间选择。

鉴于我们正在拟合多个模型，我们还必须关注研究人员的自由度（ Simmons et al., 2011 ），如果选择了单个最佳模型，则最直接来自过度拟合，或者更巧妙的是，如果我们不小心，我们可以考虑从一组拟合模型中得出的推论来包含一些总的不确定性，而无需认识到还有其他模型我们可以拟合。这种担忧出现在我们的选举预测模型中，最终我们只有少数过去的总统选举用于校准我们的预测。

图 22 说明了基本思想：该图可以表示，例如，用一系列越来越复杂的模型估计的因果效应，从简单的治疗控制比较开始，然后进行一系列调整。即使最终的兴趣只在最终模型中，了解随着调整的增加推断如何变化也很有用。

遵循建议的工作流程并探索模型的拓扑结构通常可以引导我们找到通过所有检查的多个模型。我们可以执行多元宇宙分析并使用所有模型，而不是仅选择一个模型，并查看模型之间的结论如何变化（ Steegen et al., 2016 ； Dragicevic et al., 2019 ； Kale、Kay and Hullman, 2019 ）。使用多元宇宙分析还可以减轻验证模型和先验的工作量：如果结论没有改变，那么决定哪个模型“最好”就不太重要了。图 23 显示了一个可能的输出示例。其他分析选择（数据预处理、响应分布、评估指标等）也可以进行多元宇宙分析。

图 23：来自 Niederlová et al. 的补充材料的多元宇宙分析结果。 (2019)。热图显示了关于表型（ PD、CI、REN、HEART、LIV ）与 BBSome 各种基因突变（ BBS01 - BBS8、cLOF - 功能完全丧失）之间关系的选定结论的统计评估，使用一组贝叶斯分层逻辑回归模型和成对频率论测试。基于后验预测性检查，贝叶斯模型被分类为“主要” （通过所有检查）、“次要” （某些检查中的小问题）和“有问题的拟合”，尽管我们看到大多数结论适用于所有可能的模型。模型在包含的预测变量和先验变量（默认/宽/窄 /非常窄）方面有所不同。

8.2 交叉验证和模型平均¶

我们可以使用交叉验证来比较适合相同数据的模型（ Vehtari et al., 2017 ）。在进行模型比较时，如果比较中存在不可忽略的不确定性（ Sivula et al., 2020 ），我们不应简单地选择具有最佳交叉验证结果的单一模型，因为这会丢弃交叉验证中的所有不确定性过程。相反，我们可以维护此信息并使用堆叠来组合推断，使用权重设置以最小化交叉验证错误（ Yao et al., 2018b ）。我们认为堆叠比传统的贝叶斯模型平均更有意义，因为后者可以强烈依赖于对预测影响最小的模型方面。例如，对于数据充分了解并且其参数在单位尺度上的模型，将参数先验从 \(\text{Normal}(0,10)\) 更改为 \(\text{Normal}(0,100)\) 会将边缘似然除以大约 \(10^k\) （具有 \(k\) 个参数的模型），同时保持所有预测基本相同。此外，堆叠考虑了联合预测，当候选模型列表中有大量相似但弱的模型时效果很好。

在概念上，堆叠有时可以被视为逐点的模型选择。当有两个模型并且第一个模型在 \(20%\) 的时间内优于第二个模型时，堆叠权重将接近 \((0.2,0.8)\) 。有鉴于此，堆叠填补了面向独立错误的机器学习验证与现代大数据的分组结构之间的空白。因此，模型堆叠也是模型拟合异质性的一个指标，这表明我们可以通过分层模型进一步改进聚合模型，因此堆叠是模型改进的一个步骤，而不是其本身。在极端情况下，还可以进行模型平均，以便不同的模型可以应用于不同的数据点（ Kamary et al., 2019 ； Pirš and Štrumbelj, 2019 ）。

在贝叶斯工作流程中，我们将拟合许多不感兴趣的模型包括在任何平均值中；这种 “脚手架 ” 包括那些故意过于简单的模型（包括只是为了与感兴趣的模型进行比较）和纯粹出于实验目的构建的模型，以及存在重大缺陷甚至编码错误的模型。但即使在消除了这些错误或故意的过度简化之后，在进行预测时也可能有几个模型可以对其进行平均。在我们自己的应用工作中，通常没有很多机会执行这种模型平均，因为我们更喜欢持续地模型扩展，但是在某些场景中，用户会合理地希望对竞争的贝叶斯模型进行平均预测，例如 Montgomery and Nyhan, 2010。

8.3 比较大量模型¶

有很多问题，例如在具有多个潜在相关预测变量的线性回归中，有许多候选模型可用，所有这些都可以描述为单个扩展模型的特殊情况。如果候选模型的数量很大，我们通常有兴趣找到一个相对较小的模型，该模型与我们的扩展模型具有相同的预测性能。

这导致了预测器（变量）选择的问题。如果我们有许多模型做出类似的预测，那么根据最小化交叉验证误差来选择其中一个模型会导致过度拟合和次优模型选择（ Piironen and Vehtari, 2017 ）。相比之下，投影预测变量选择在寻找具有良好预测性能的较小模型方面已被证明是稳定可靠的（ Piironen and Vehtari, 2017 ； Piironen et al., 2020 ； Pavone et al., 2020 ）。虽然搜索大模型空间通常与过度拟合的危险相关，但投影预测方法通过仅检查基于扩展模型的预测的投影子模型而不是将每个模型独立地拟合到数据来避免这种情况。

除了变量选择，投影预测模型选择可用于广义加性多级模型中的结构选择（ Catalina et al., 2020 ）以及为复杂的非参数模型创建更简单的解释（ Afrabandpey et al., 2020 ）。

9 把建模当做软件开发¶

用概率编程语言开发统计模型意味着编写代码，因此是一种软件开发形式，有几个阶段：编写和调试模型本身；将数据转化为合适的建模形式所必需的预处理；以及随后的理解、交流和使用由此产生的推论的步骤。开发软件很难。很多事情都可能出错，因为有很多活动部件需要仔细同步。

软件开发实践旨在减轻由编写计算机程序的固有复杂性引起的问题。不幸的是，许多方法转向教条、豆数或两者兼而有之。我们可以推荐的一些参考文献是 Hunt and Thomas, 1999 以及 McConnell , 2004 ，它们为开发人员提供了可靠、实用的建议。

9.1 版本控制使与他人和你过去自己的协作顺利进行¶

版本控制软件，例如 Git，应该是项目的第一个基础设施。学习版本控制似乎是一项巨大的投资，但能够键入单个命令以恢复到以前的工作版本或获取当前版本和旧版本之间的差异是非常值得的。当你需要与他人共享工作时甚至是在纸上共享工作时效果更好 - 工作可以独立完成然后自动合并。虽然版本控制会跟踪一个模型中较小的更改，但将明显不同的模型保存在不同的文件中以方便比较模型是很有用的。版本控制还有助于记录迭代模型构建中的发现和决策，从而提高过程的透明度。版本控制不仅仅针对代码。它也适用于报告、图表和数据。

版本控制不仅仅针对代码。它也适用于报告、图表和数据。版本控制是确保所有这些组件同步的关键部分，更重要的是，可以将项目回滚到以前的状态。版本控制特别有用，因为它能够打包和标记与里程碑报告和出版物相对应的模型和数据的“候选发布”版本，并将它们存储在同一目录中，而无需求助于可怕的 _final_final_161020.pdf 样式命名约定。

在处理用于指导政策决策制定的模型时，公共版本控制存储库提高了特定报告使用了哪些模型、数据、推断参数和脚本的透明度。一个很好的例子是帝国理工学院的模型和脚本存储库，用于估计 COVID-19 的死亡和病例数（ Flaxman et al., 2020 ）。

9.2 边运行边测试¶

理想情况下，软件设计从最终用户的目标到实现它所需的技术机制自上而下地进行。对于贝叶斯统计模型，自上而下的设计至少涉及数据输入格式、数据的概率模型和先验，但也可能涉及模拟器和模型测试，如基于模拟的校准或后验预测性检查。理想情况下，软件开发从经过良好测试的基础功能到更大的功能模块自下而上地工作。这样，开发将通过一系列经过充分测试的步骤进行，在每个阶段都只在经过测试的部分上进行构建。与构建大型程序然后调试它相比，以这种方式工作的优势与增量模型开发相同——更容易跟踪开发出错的地方，并且你在使用经过良好测试的基础工作的每一步都更有信心。

无论是初始开发还是修改代码，计算开发的关键是模块化。大杂乱的功能难以记录、难以阅读、异常难以调试，并且几乎不可能维护或修改。模块化意味着从较小的可信部分构建更大的部分，例如低级功能。每当代码片段重复时，都应该将它们封装为函数。这导致代码更易于阅读和维护。

作为低级函数的示例，可以通过实施标准化 \(z(v) = (v −\text{mean}(v))/\text{sd}(v)\) 为广义线性模型重新调整预测变量。虽然这个函数看起来很简单，但从 \(\text{sd}\) 函数开始，它有时被定义为 \(\text{sd}(v) =\sqrt{\sum^n_{i=1}(v_i −\text{mean}(v))^2/n}\) ，有时被定义为 \(\text{sd}(v) = \sqrt{\sum^n_{i=1}(v_i - \text{mean}(v))^2/(n -1)}\) 。如果在标准化函数层面不解决这个问题，在推断过程中会产生神秘的偏差。不依赖于 \(\text{sd}()\) 函数的简单测试将在函数开发过程中解决这个问题。如果选择是除以 \(n -1\) 的估计值，则需要决定当 \(v\) 是长度为 \(1\) 的向量时要做什么。在存在非法输入的情况下，在输入-输出例程中进行检查会有所帮助让用户知道什么时候出现问题，而不是让错误渗透到以后神秘的除以零错误。

微分方程的三次样条或欧拉求解器的实现是高级函数的一个示例，应在使用前对其进行测试。随着函数变得越来越复杂，由于边界条件组合学的问题、更一般的输入（例如要积分的函数）、数值不稳定性或不精确性，其域的区域可能会或可能不会接受，因此它们变得更难测试，这取决于应用，需要稳定的衍生品等。

9.3 使其基本上可重现¶

任何项目的一个崇高目标是使其完全可重现，因为另一台机器上的另一个人可以重新创建最终报告。这不是科学领域考虑的可重复性类型，科学领域希望确保影响得到未来新数据的证实（现在通常称为“可重复性”，以便更好地区分不同概念）。相反，这是确保始终如一地完成一项特定分析的更有限（但仍然至关重要）的目标。特别是，我们希望能够生成与原始文档基本相同的分析和图表。位级可再现性可能是不可能的，但我们仍然会在实际水平上将等效性进行比较。如果这种类型的复制改变了论文的结果，我们会争辩说原始结果不是特别可靠。

与其在运行模型时在命令行上输入命令（或直接将命令输入到 R 或 Python 等交互式编程语言中），不如编写脚本来通过模型运行数据并生成你需要的任何后验分析。可以为 shell、R、Python 或任何其他编程语言编写脚本。

脚本应该是独立的，因为它应该在一个完全干净的环境中运行，或者最好在不同的计算机上运行。这意味着脚本不得依赖于已设置的全局变量、正在读入的其他数据或脚本中没有的任何其他内容。脚本是很好的文档。如果运行项目只是一行代码，这似乎有点过分，但脚本不仅提供了一种运行代码的方法，而且还提供了一种关于正在运行的内容的具体文档的形式。对于复杂的项目，我们经常发现构建良好的一系列脚本比一个大的 R Markdown 文档或 Jupyter notebook 更实用。

根据长期的再现性需求，为手头的工作选择合适的工具很重要。为了保证位级的可重复性，有时甚至只是为了让程序运行，从硬件到操作系统，再到每一个软件和设置，都必须用它们的版本号来指定。随着脚本的初始编写和复制尝试之间的时间流逝，即使环境随脚本一起提供，也几乎不可能实现位级可复制性，就像在 Docker 容器中一样。

9.4 使其具有可读性和可维护性¶

像对待其他形式的写作一样对待程序和脚本为观众提供了如何使用代码的重要视角。不仅其他人可能想要阅读和理解一个程序或模型，开发人员自己以后也会想要阅读和理解它。 Stan 设计的动机之一是让模型在变量使用（例如，数据或参数）、类型（例如，协方差矩阵或无约束矩阵）和大小方面进行自我记录。这使我们能够理解 Stan 代码（或其他静态类型概率编程语言的代码），以便在没有应用它的数据上下文的情况下也能被理解。

可读性的很大一部分是一致性，特别是在命名和布局方面，不仅是程序本身，还有存储它们的目录和文件。编码的另一个关键原则是避免重复，而是将共享代码提取到可以重用的函数中。

代码的可读性不仅仅与注释有关——它还与可读性的命名和组织有关。事实上，注释会使代码的可读性降低。最好的方法是编写可读的代码，而不是带有注释的不透明代码。例如，我们不想这样写：

real x17; // 氧气水平，应该是正的

当我们可以这样写时：

real<lower = 0>oxygen_level;

同样，我们也不想这样做：

target += -0.5 * (y - mu)^2 / sigma^2; // y 分布 normal(mu, sigma),

当我们可以写：

target += normal_lpdf(y  \mid  mu, sigma);

好的做法是尽量减少内联代码注释，而是编写可读代码。如上述示例所示，编程语言为用户提供了他们需要使用的工具，从而促进了干净的代码。

面向用户的函数应该在函数级别记录其参数类型、返回类型、错误条件和行为——这是用户看到的应用程序编程接口 (API)，而不是代码内部。针对开发人员的内联代码注释的问题在于，它们在开发过程中很快就会变得陈旧，最终弊大于利。相反，与其在内联记录实际代码，不如将函数减少到可管理的大小，并应选择名称以便代码可读。较长的变量名并不总是更好，因为它们会使代码结构更难以扫描。编写代码文档时应该假设读者很好地理解编程语言；因此，只有在代码偏离了语言的惯用用法或涉及复杂算法时才需要文档。当试图对长表达式或代码块进行注释时，请考虑将其替换为一个命名良好的函数。

Python贝叶斯分析(中文)

附录 K：贝叶斯工作流程

Contents

附录 K：贝叶斯工作流程¶

1 概述¶

1.1. 从贝叶斯推断到贝叶斯工作流程¶

1.2 为什么需要贝叶斯工作流程？¶

1.3 工作流程及其与统计理论和实践的关系¶

1.4 组织贝叶斯工作流程的诸多方面¶

1.5. 文章的目的和结构¶

2 在拟合模型之前¶

2.1 选择初始模型¶

2.2 模块化构建¶

2.3 对观测数据进行缩放和转换¶

2.4 尝试一下先验预测性检查¶

2.5 生成式和部分生成式模型¶

3 拟合一个模型¶

3.1 初值、适应和预热¶

3.2 运行迭代算法需要多长时间？¶

3.3 近似算法和近似模型¶

3.4 快速拟合，快速失败¶

4 使用人为构造的伪数据发现和理解问题¶

4.1 伪数据模拟¶

4.2 基于模拟的较正¶

4.3 使用构造的数据进行实验¶

5 解决统计计算中存在的问题¶

5.1 统计计算的通俗理论¶

5.2 从简单和复杂的模型开始并在中间相遇¶

5.3 应对需要很长时间才能拟合的模型¶

5.4 监控中间量¶

5.5 堆叠以重新加权混合不良的链¶

5.6 如何应对具有复杂几何形态的后验分布 ？¶

5.7 通过重新参数化来适应复杂的后验分布¶

5.8 通过对联合概率分布的结构化和边缘化降低计算难度¶

5.9 增加先验信息¶

5.10 增加数据¶

6 评估和使用拟合后的模型¶

6.1 后验预测性检查¶

6.2 单数据点和数据子集的交叉验证和影响¶

6.3 先验信息的影响¶

6.4 对推断结果进行总结¶

7 修改模型¶

7.1 为数据构建模型¶

7.2 合并额外的数据¶

7.3 使用先验分布¶

7.4 模型拓扑¶

8. 理解和比较多个模型¶

8.1 可视化彼此相关的模型¶

8.2 交叉验证和模型平均¶

8.3 比较大量模型¶

9 把建模当做软件开发¶

9.1 版本控制使与他人和你过去自己的协作顺利进行¶

9.2 边运行边测试¶

9.3 使其基本上可重现¶

9.4 使其具有可读性和可维护性¶

10. 包含模型构建和模型扩展的工作流程示例：高尔夫推杆¶

10.1 第一个模型：逻辑回归¶

10.2 从第一原则建模¶

10.3 在新数据上测试拟合后的模型¶

10.3 在新数据上测试拟合后的模型¶

10.4 一种解释击球力度的新模型¶

也就是说，除了角度满足要求之外，必须足够用力才能到达球洞，但又不能走得太远。条件（ a ）我们之前考虑到了；现在需要添加条件（ b ）。¶

10.5 通过加入一个模糊因子来扩展模型¶

10.6 高尔夫示例的一般教训¶

11 具有不可预期的多峰后验的工作流程示例：行星运动¶

11.1 运动的机械模型¶

11.2 拟合一个简化模型¶

11.3 坏马尔可夫链，慢马尔可夫链？¶

11.4 建立模型¶

11.5 行星运动示例的一般经验教训¶

12 讨论¶

12.1 统计建模和预测的不同视角¶

12.2 迭代式模型构建的必要性¶

12.3 模型选择和过拟合¶

12.4 更大的数据集需要更大的模型¶

12.5 预测、泛化和后分层化¶

12.6 继续往前走¶

参考文献¶

引文信息¶

5.6 如何应对具有复杂几何形态的后验分布？¶