第六章 近似概率推理 《贝叶斯认知模型  逆向工程思维》

第六章：近似概率推理 《贝叶斯认知模型 逆向工程思维》

https://dokumen.pub/bayesian-models-of-cognition-reverse-engineering-the-mind-1nbsped-9780262049412-9780262381048-9780262381055.html

总结

在许多方面，近似推断算法是贝叶斯建模的隐藏秘密，因为在几乎所有有趣的概率模型中，精确推断都是不可能的。理解不同的推断算法及其适用情境是研究这些模型的研究人员工具包的重要组成部分。通常，推断的便捷性是驱动模型选择的一个因素：写出一个极其复杂的生成模型很容易，但无法保证能够应用贝叶斯规则来计算该模型中的后验分布。最近的创新，例如高效概率编程语言的发展（参见第18章），可以使处理表达力强的概率模型变得更加容易。使用这些语言的好处之一是，用于指定概率模型的代码可以自动应用于该模型的不同推断算法。然而，基本了解这些算法是什么以及它们如何工作，对于理解它们何时会失败至关重要。

近似推断算法还为人类心智和大脑如何可能执行贝叶斯推断这一问题提供了答案。正如研究人员在贝叶斯模型中有多种算法可选以近似不同假设的后验概率一样，我们的心智和大脑也可能有许多可能的策略可用于近似贝叶斯推断。

我们将在本书后面重新探讨这一主题。第11章考虑蒙特卡洛方法如何解释人类概率推理的某些方面，第13章通过考虑心智和大脑如何最有效地利用此类计算资源来扩展这一想法，而第12章则将本章介绍的一些算法与经典和近期关于人工神经网络的工作联系起来。

当我们应用贝叶斯推理时，需要评估的概率分布很快就会变得非常复杂。图形模型提供了一些加速概率推理的工具，但这些工具在大多数变量仅直接依赖于相对较少的其他变量时效果最佳。然而，在许多情况下，计算后验分布或条件概率仍然是一项计算上的挑战。例如，当我们无法计算贝叶斯定理分母所需的积分时，如何对连续变量进行贝叶斯推理？或者，如何在庞大的离散假设空间（例如，定义在一组变量上的所有有向图构成的空间）上计算后验分布？

概率推理的计算复杂性带来了两类问题。第一类是工程问题：作为建模者，我们如何计算我们所构建的概率模型的预测结果？第二类是逆向工程问题：如果概率推理如此难以计算，人类行为又如何可能与这些模型的预测结果保持一致？在本章中，我们将介绍一些近似概率计算的方法，这些方法对上述两个问题都具有相关性，但我们将重点关注工程问题。我们将在第11章详细探讨逆向工程问题。

近似概率推理的方法多种多样，包括通用的数值积分或优化方法，以及专门为概率推理开发的方法，例如环状信念传播（loopy belief propagation；Pearl, 1988）和变分推理（variational inference；Jordan, Ghahramani, Jaakkola & Saul, 1999）。所有这些方法对于解决使用概率模型时的工程问题都具有重要意义，其中许多方法也可能有助于阐明人类心智和大脑如何执行概率推理。本章的重点将放在蒙特卡洛方法（Monte Carlo methods）上——这类方法通过用来自概率分布的样本来替代该分布本身，从而近似概率计算。蒙特卡洛方法是近似推理的基本工具之一，可适用于各种各样的概率模型。它们还能为人类心智和大脑如何应对概率推理中严峻的计算挑战提供清晰的假设来源。在本章末尾，我们将简要讨论变分推理，该方法近年来日益趋向于优化。关于蒙特卡洛方法的更多细节，我们推荐Neal（1993）、MacKay（1998）以及Robert和Casella（1999）；关于变分推理的入门介绍，可参见Blei、Kucukelbir和McAuliffe（2017）。

6.1 简单蒙特卡洛

蒙特卡洛方法的基本思想是用一组来自某概率分布的样本来表示该分布。这些样本可以指示哪些取值具有较高的概率（因为高概率的取值更有可能被抽样得到），并且在进行各种计算时，可以用这些样本代替原始分布本身。在处理认知的概率模型时，我们通常感兴趣的是理解关于参数化模型的后验分布——例如，一个带有因果强度参数的因果网络——或者关于一类模型的后验分布——例如，定义在一组变量上的所有因果网络结构空间，或定义在一组对象上的所有分类树结构空间。来自后验分布的样本有助于发现模型的最佳参数值或模型类中的最佳模型，并可用于估计后验分布在这些最佳假设上的集中程度（即学习者应对这些假设抱有多大信心）。

抽样也可用于近似后验分布上的平均值。例如，在计算给定数据下某个参数化模型的后验概率时，需要计算该模型的边缘似然（marginal likelihood），即数据在模型所有参数设置下的平均概率。在进行理想的贝叶斯预测（例如计算加权硬币的后验预测分布）时，同样需要对所有参数设置进行平均。最后，我们可能还希望对模型结构空间进行平均，从而对那些无论哪种结构正确都可能成立的模型特征做出预测。例如，我们可以估计在一个结构未知的复杂因果网络中，变量A导致变量B的可能性有多大——方法是计算在从网络结构后验分布中抽取的高概率样本中，存在A → B连接的概率（例如，Friedman & Koller, 2000）。

6.1.1 一个例子：计算因果支持度

为了说明蒙特卡洛方法如何用于近似数值积分并有助于评估贝叶斯模型，我们回顾第4章中介绍的因果结构学习模型——因果支持度（Griffiths & Tenenbaum, 2005）。为了计算一组共现数据 d 对某一因果关系的支持证据，我们需要计算如下积分：

6.2 简单蒙特卡洛何时会失效？

简单蒙特卡洛方法的一个局限在于，它并不容易自动地从大多数概率分布中生成样本。使用蒙特卡洛方法要求我们能够从分布 p(x) 中生成样本。然而，情况往往并非如此。对于均匀分布、高斯分布、指数分布、泊松分布、伽马分布和贝塔分布等分布，存在直接的算法用于生成随机样本，但这些分布属于行为特别良好的分布（事实上，它们都属于指数族分布，如第3章所述）。其他一些分布可能需要巧妙的方法才能从中采样，开发针对此类分布的抽样方法本身就可以成为一个研究课题。

当我们想用蒙特卡洛方法从通过应用贝叶斯法则计算出的后验分布中采样时，往往会遇到一类非常具体的问题，使得简单蒙特卡洛方法难以适用。

在处理后验分布时，我们可能会陷入这样一种境地：我们实际上并不知道某个特定假设的概率——我们只知道该概率的值（最多相差一个常数）。如果我们有一个先验分布 P(h)、一个似然函数 p(d|h)，以及一个观测数据点 d，则后验分布由下式给出：

其中分母中的求和遍历所有假设。由于分母的作用是确保 P(h|d) 是“归一化”的——即它是一个定义良好的概率分布，所有假设的概率之和为1——这个求和通常被称为归一化常数。当假设空间很大且离散时——例如，n 个对象的所有划分，或 n 个变量上的所有有向图——计算归一化常数会迅速变得不可行。如果分子容易计算，我们会陷入这样一种情形：我们能够计算与 P(h|d) 成比例的某个量，并可比较不同假设下这些概率的相对大小，但我们只知道它们的实际概率（最多相差一个常数）。

当我们对连续变量应用贝叶斯法则时，同样的问题也会出现。在这种情况下，我们有：

其中分母中的积分可能没有解析解。在这种情况下，我们可以尝试对该积分进行数值近似——但随着 θ 的维度增加，这会变得越来越困难——或者我们可以依赖那些只需要知道 p(θ|d) 最多相差一个常数的方法。

在这些情况下，简单蒙特卡洛不是一个可行的解决方案——我们甚至不知道想要从中抽样的概率分布是什么，更不用说如何从中抽样了。然而，确实存在一些更复杂的蒙特卡洛方法，即使我们只知道概率分布最多相差一个常数，也可以应用。

6.3 拒绝采样

假设我们有一个难以从中采样的分布 p(x)——可能是因为我们只知道它的值（最多相差一个常数）。如果另一个分布与 p(x) 接近但易于采样，那么通常可以使用拒绝采样（rejection sampling）来从 p(x) 中生成样本。

改变采样问题 假设存在一个分布 q(x)，我们可以很容易地从中采样，我们将称其为“提议分布”（proposal distribution）。此外，假设我们能够计算概率 p(x)（最多相差一个常数），即我们可以计算 p*(x) ∝ p(x)。再进一步假设，我们知道一个常数 c，使得对所有 x 都满足 c q(x) ≥ p*(x)。那么，我们可以通过执行以下步骤来生成一个样本：

1. 从 q(x) 中抽取样本 x；
2. 从区间 [0, c q(x)] 中均匀抽取一个随机实数 r；
3. 如果 r < p*(x)，则保留 x；否则，拒绝该样本并回到第1步。

最终得到的结果是从 p(x) 中抽取的一个样本。

为什么这种方法有效？当 p(x) 是一维（1D）分布时，这一点可以直观理解（参见图6.1，其中 q(x) 被取为均匀分布）。两个随机值 x 和 r 构成了二维（2D）平面上某点的坐标。该点是从包含 p*(x) 曲线的一个区域中均匀随机抽取的（实际上，该区域位于曲线 c q(x) 之下）。如果我们只保留位于 p*(x) 曲线下方的点，那么我们抽取每个值 x 的概率就与 p*(x) 成比例。由于 p*(x) ∝ p(x)，这等价于从分布 p(x) 中进行采样。

通过拒绝采样进行贝叶斯推理 拒绝采样可用于从后验分布中采样，从而近似贝叶斯推理。我们将从数据和假设均为离散的情形开始，因此我们希望从 P(h|d) 中采样。

第一步是构造一个用于采样的分布。我们可以将后验写成关于假设 h 和数据 d′ 的联合分布：P*(h, d′) = P(h) P(d′|h) δ(d′, d)，其中 δ(·, ·) 在其参数匹配时为1，否则为0。这看起来有点奇怪，但实际上它相当于这样的假设：我们首先从先验中采样假设 h，然后从似然函数 P(d′|h) 中采样数据 d′，但只保留那些与观测数据 d 匹配的数据。这定义了一个在观测数据上集中全部质量的分布，虽然技术上它是关于 h 和 d′ 的联合分布。

如果我们使用提议分布 Q(h, d′) = P(h) P(d′|h)，那么显然有 P*(h, d′) ≤ Q(h, d′)，因为 δ(·, ·) 只能取值 0 或 1。因此，c = 1。此外，由于当 d′ = d 时 P*(h, d′) 等于 Q*(h, d′)，而在其他情况下为零，因此当且仅当 d′ ≠ d 时我们会拒绝一个样本。

这提供了一种通过拒绝采样从后验分布中抽样的简单方法：从先验中抽取假设 h，然后从与该假设相关的似然函数 P(d′|h) 中抽取数据 d′，如果 d′ 与观测数据 d 不匹配，则拒绝该假设。拒绝某个假设 h 的概率即为 P(d|h)，而我们保留下来的假设将是来自后验分布的样本。

扩展与局限性 该方法自然可推广到连续假设的情形：我们从先验分布 p(θ) 中采样，并以概率 P(d|θ) 保留这些样本。然而，当数据是连续的时候，该方法不可用，因为这些数据的概率会变得无穷小。在这种情况下，研究人员使用了一种相关的方法，称为近似贝叶斯计算（ABC；Beaumont, Zhang, & Balding, 2002），其中引入第二个函数来定义数据集之间的接近程度度量。然后从先验中抽取样本，用以模拟数据，只有当模拟数据与观测数据足够接近时才保留样本。

尽管拒绝采样很简单，但它通常效率不高。因为我们只以概率 P(d|h) 保留每个样本，所以保留样本的比例由 P(d|h) 关于 P(h) 的期望给出，即 Σₕ P(d|h) P(h) = P(d)。因此，该方法仅在观测数据的概率较高时才高效——通常意味着观测数据量很小。然而，拒绝采样是近似贝叶斯推理的良好直觉来源，并可作为定义在复杂模型中执行贝叶斯推理策略的起点。

6.4 重要性采样

拒绝采样能够精确地从目标分布 p(x) 中生成样本，但它可能非常低效。另一种从相同思想出发的方法——即先从错误的分布中生成样本，再对其进行校正——被称为重要性采样（importance sampling）。

改变采样问题 假设我们希望计算某个函数 f(x) 关于概率分布 p(x) 的期望值，但我们有另一个更容易从中采样的分布 q(x)。此外，假设当 p(x) > 0 时，总有 q(x) > 0。

我们可以通过乘以并除以同一项将 q(x) 引入我们的期望表达式中，从而得到：

尽管存在偏差，但在简单蒙特卡洛方法可以应用的情况下，使用重要性采样器仍有若干合理的原因。首先，它允许使用同一组样本，通过为每个分布赋予不同的权重，来评估关于一系列分布的期望值。其次，重要性采样器所产生的估计值可能比简单蒙特卡洛方法产生的估计值具有更低的方差。如果函数 f(x) 在 p(x) 较小的区域取到其最极端的值，则简单蒙特卡洛估计的方差可能很大。若选择 q(x) 与 f(x) 互补，重要性采样器的方差可低于简单蒙特卡洛方法。特别地，当指定 q(x) 为：

时，采样器的渐近方差达到最小。

这不是一个实用的程序，因为找到该分布需要计算 Eₚ₍ₓ₎[f(x)]；但这一事实表明，最小方差采样器不必是 p(x)，从而说明重要性采样有时能够提供比简单蒙特卡洛更低方差的期望估计值。

通过重要性采样进行贝叶斯推理 我们可以使用重要性采样来进行贝叶斯推理。我们希望近似的分布是后验分布 p(h|d)。我们可以通过选择另一个在假设空间上的分布 q(h) 来生成样本，并为所得样本分配与 p(h|d)/q(h) 成比例的权重来实现这一点。这种方法的一个吸引人的特性是，由于权重仅与该比率成比例，我们也可以等价地分配与 p(d|h)p(h)/q(h) 成比例的权重。这为我们提供了一种无需计算归一化常数即可从后验分布生成样本的方法——显著降低了计算成本。

此处对重要性采样的偏差和方差的分析清楚地表明，当 q(x) 和 p(x) 足够相似、使得权重不会变化太大时，该方法最为有效。因此，使用重要性采样进行贝叶斯推理在以下情况下最有效：如果我们能找到一个接近后验分布 p(h|d) 的提议分布 q(h)。在某些情形下——当观测数据对代理信念的影响较小时，对应于似然函数 p(d|h) 不取极端值——我们可以直接使用先验分布 p(h) 作为我们的提议分布。在这种情况下，分配给每个样本的权重就是 p(d|h)p(h)/p(h) = p(d|h)，即似然值。由此得到的用于近似贝叶斯推理的算法极为简单：从先验分布 p(h) 中生成样本，并为这些样本分配与似然 p(d|h) 成比例的权重。出于显而易见的原因，该算法被称为“似然加权”（likelihood weighting）。

6.4.1 一个例子：计算空间巧合

Griffiths 和 Tenenbaum（2007a）提出的关于空间巧合的贝叶斯模型（在第3章中已讨论）要求计算一个复杂的积分。一般来说，应用贝叶斯奥卡姆剃刀的挑战在于：要计算某个特定模型下所有观测数据的边缘概率，需要对该模型的所有参数进行积分。在空间巧合的情形下，这意味着需要对混合模型的参数进行积分。

在空间巧合场景中所假设的模型相对简单：要么所有炸弹随机坠落（h₀），要么其中一些炸弹围绕某个目标坠落（h₁）。如果 nB 颗炸弹在某个区域 ℛ 内随机坠落，则任何一组炸弹位置 x 的概率仅为 p(x|h₀) = 1/|ℛ|ⁿᴮ。如果所有炸弹都围绕目标坠落，那么评估起来也相对容易——我们希望计算在高斯分布下 x 的边缘概率，即对高斯分布的均值 μ 和协方差矩阵 Σ 进行积分：p(x|h₁) = ∫∫ p(x|μ, Σ) p(μ) p(Σ) dμ dΣ。尽管这听起来令人望而生畏，但如果对 μ 和 Σ 使用共轭先验（在此情况下，μ 可以是正态分布或均匀分布，Σ 则使用逆威沙特先验），实际上操作起来是直接的（详见 Minka, 2001）。

棘手的部分出现在某些炸弹被瞄准、而其他炸弹随机坠落的情况下。此时我们面对的是一个混合模型，其中一个成分在 ℛ 上服从均匀分布，另一个成分是具有未知均值和方差的高斯分布。如果我们知道哪些炸弹是被瞄准的，问题仍然简单——我们可以对被瞄准炸弹的概率在 μ 和 Σ 上积分，而其余炸弹则服从 ℛ 上的均匀分布。令 z 表示将炸弹分配到混合成分的向量，我们便可以轻松计算 p(x|z)。但问题是，我们并不知道 z，因此我们需要对 z 所有可能取值求和。

写出这一表达式，我们有

其中 P(z) 是关于分配的先验分布。问题是，如果有 nB 颗炸弹，则 z 有 2ⁿᴮ 种可能取值——这个数目会迅速变得非常大。使用简单蒙特卡洛方法是一个诱人的解决方案——我们可以直接从 P(z) 中抽取 z 的样本，然后对得到的 p(x|z) 值求平均。不幸的是，在某些炸弹确实被瞄准的情形下，这种方法效果不佳，因为大多数 z 的选择不会挑出被瞄准的炸弹，从而导致 p(x|z) 的值非常小；而少数选择会接近真实情况，从而使 p(x|z) 的值非常大。简单蒙特卡洛方法得出的结果将更多地取决于抽样的运气，而非实际答案。

对于混合模型等复杂模型计算边缘概率，是重要性采样极为有效的场景之一。诀窍在于选择一个提议分布 Q(z)，使其在 p(x|z) 较大的那些 z 取值上赋予比 P(z) 更高的概率。一种已提出的实现方法是：使用期望最大化（EM）算法等推理算法来估计混合模型的参数，然后利用这些参数定义 Q(z)。在此情形下，我们可以通过尝试找到每个轰炸场景中似乎对应目标位置的值，来估计 μ 和 Σ，然后计算 P(z|x, μ, Σ) 并将其作为 Q(z)。

这一方法存在一些风险——它会选出一个可能的目标，但可能存在其他同样合理的可能性——因此标准做法是将 Q(z) 设为该分布与先验 P(z) 的混合。这就是用于计算与人类行为进行比较的概率所采用的方法（想出如何实现并运行该算法所需的时间甚至比执行实验和分析数据还要长！）。

6.5 序贯蒙特卡洛

在认知的概率模型中，一个基本问题是如何解释学习者如何应对概率推理所带来的严峻计算挑战。我们已经讨论了在静态情境下——即所有数据都可用、学习者只需评估一组假设——如何用蒙特卡洛方法解决其中部分问题。然而，当需要动态概率推理时（例如由于模型内在动力学特性，或因更多数据逐渐可用），挑战似乎变得更加严峻。这些问题可以使用序贯蒙特卡洛方法加以解决。

6.5.1 应对动态世界

序贯蒙特卡洛方法是从一系列概率分布中生成样本的方案，通常利用连续分布之间的关系，使用前一个分布的样本来生成下一个分布的样本。此类方法的标准例子是粒子滤波器（particle filter），它是重要性采样的序贯版本。粒子滤波器最初是为了在动态环境中对变量进行推理而开发的——“滤波”问题是指根据一系列观测值推理当前世界状态的问题。

假设我们有一个结构类似于隐马尔可夫模型（HMM）的概率模型，其中一串潜在变量 z₁, …, zₙ 负责生成一串观测变量 x₁, …, xₙ。在最简单的情况下，每个 zᵢ 从仅依赖于 zᵢ₋₁ 的分布中生成，而每个 xᵢ 从仅依赖于 zᵢ 的分布中生成。假设我们希望计算 p(zₙ | x₁, …, xₙ)。虽然当 zᵢ 离散时（如 HMM 的前向算法），或当 zᵢ 连续且 p(zᵢ | zᵢ₋₁) 是均值为 zᵢ₋₁ 线性函数的高斯分布时（第5章讨论的卡尔曼滤波器），已有高效的算法可用于计算该概率，但粒子滤波器提供了一种能够推广到这些情形之外的推理方法。

粒子滤波器背后的基本思想是，我们将使用重要性采样来近似 p(zₙ | x₁, …, xₙ)。具体而言，我们可以进行重要性采样，其中我们以“先验”p(zₙ | x₁, …, xₙ₋₁) 作为提议分布。如果我们能从 p(zₙ | x₁, …, xₙ₋₁) 生成样本，那么我们可以通过赋予每个样本与 p(xₙ | zₙ) 成比例的权重来构建一个重要性采样器。我们如何从 p(zₙ | x₁, …, xₙ₋₁) 生成样本呢？嗯，我们知道：

因此，如果我们能从分布 p(zₙ₋₁ | x₁, …, xₙ₋₁) 中生成带权重 wⱼ 的样本 zₙ₋₁⁽ʲ⁾，我们就可以直接从以下分布中采样：

这是一个简单的混合分布。或者，我们也可以对采样进行分层处理，确保样本具有更高的多样性，即对每个 zₙ₋₁⁽ʲ⁾，我们从 p(zₙ | zₙ₋₁⁽ʲ⁾) 中生成一个样本。

前一段所述的程序揭示了一个有趣的递归关系——只要我们能够从 p(zₙ₋₁ | x₁, …, xₙ₋₁) 中采样，它就为我们提供了一种从 p(zₙ | x₁, …, xₙ) 中采样的方法。这使我们能够引入支撑粒子滤波器的序贯蒙特卡洛方案（参见图6.3）。首先，我们从 p(z₁ | x₁) 中生成样本。然后，对于每个样本 z₁⁽ʲ⁾，我们从 p(z₂ | z₁⁽ʲ⁾) 中生成 z₂⁽ʲ⁾，并为每个所得样本分配一个与 p(x₂ | z₂⁽ʲ⁾) 成比例的权重 wⱼ。此步骤等价于前面介绍的似然加权过程。我们对所有 i = 3, …, n 重复这一过程，将每个样本 zᵢ⁽ʲ⁾（现称为“粒子”）的旧权重乘以 p(xᵢ | zᵢ⁽ʲ⁾)。在每一步都不需要对粒子权重进行归一化——这可以在整个过程结束时再进行。

这种简单的递归方案被称为序贯重要性采样。然而，它存在一个重大问题：随着时间推移，由于某些粒子最终可能获得一系列极不可能的 zᵢ 值，其权重可能会显著发散。在某种程度上，这是一种计算资源的浪费，因为那些权重极小的粒子对 zₙ 分布的最终表示贡献甚微。

为解决这一问题，我们可以采用一种称为“序贯重要性重采样”（sequential importance resampling）的替代方法。在此方法中，我们监控归一化粒子权重的方差。当该方差超过预设阈值时，我们从与归一化权重相对应的概率分布中抽取一组新的粒子。这增加了处于空间良好区域内的粒子数量。随后，权重被重置为在所有粒子上均匀分布，过程继续进行。重采样会降低粒子集合中的多样性，因此如果由 p(zᵢ|zᵢ₋₁) 决定的内在动力学并未产生大量变异性，则可能需要采用其他引入多样性的方案（例如扰动粒子的取值）。

尽管标准粒子滤波算法遵循这一方案，但这仅触及了可能的序贯蒙特卡洛技术的表面（有关全面综述，参见 Doucet, de Freitas, & Gordon, 2001；Chopin & Papaspiliopoulos, 2020）。即使在粒子滤波框架内，也有许多选项可用于提升特定问题上的性能。例如，在某些情况下，可以直接计算 p(zᵢ|x₁, …, xᵢ₋₁)，此时从该分布采样效果更好。在其他情况下，我们可以在提出 zᵢ 的取值时考虑观测值 xᵢ 的价值，从而提高提议的准确性并降低权重的方差。与其他蒙特卡洛方法一样，针对特定概率模型定制算法是一门艺术。

6.5.2 累积证据

虽然处理动态模型很重要，但一个更为基础的情境是：当我们希望根据新证据更新关于假设的后验分布时，必须执行具有挑战性的概率计算。这要求重新计算每个单独假设的概率，当假设数量庞大时，计算量可能极其巨大。然而，序贯蒙特卡洛在此情境下同样有用。

想象我们有一组假设 h ∈ ℋ 和一系列观测数据 d₁, …, dₙ。我们希望设计一个高效的方案来更新我们的信念，使我们能够在每次获得新观测 dᵢ 后逐步计算 p(h|d₁, …, dₙ)（即得到每次观测后的后验分布），而无需每次都在全部 |ℋ| 个假设上求和。再次想象执行重要性采样，其中我们以 p(h|d₁, …, dₙ₋₁) 作为提议分布。如果 dᵢ 在给定 h 的条件下相互独立，那么我们可以直接赋予所得 h 样本以与 p(dₙ|h) 成比例的权重。因此，与之前类似，只要我们能从 p(h|d₁, …, dₙ₋₁) 中获取样本，我们就有了从 p(h|d₁, …, dₙ) 中获取样本的方法。这使我们能够建立类似于粒子滤波器中所使用的那种递归结构。

最简单的序贯蒙特卡洛方案如下：首先，从先验分布 p(h) 中生成样本。然后，每当新的观测 dᵢ 到来时，按 p(dᵢ|h) 重新加权这些样本。当我们需要计算后验分布时，可以对权重进行归一化。这潜在地为我们提供了一组带权重的假设，作为每次 dᵢ 之后后验分布的近似。

尽管该方案因其简洁性而具有吸引力，但它存在一个基本问题：我们在近似中使用的假设将始终相同，无论我们观察到的数据如何——对应于我们从先验中抽样的那组假设。如果经过 d₁, …, dₙ 之后的后验分布与先验分布差异很大，这将导致近似效果很差。因此，我们需要引入诸如扰动粒子之类的技术，以便在每次迭代中生成不同的假设。

在扰动粒子时，我们的目标是不改变粒子所来自的分布——如果它们在扰动前是从 p(h|d₁, …, dᵢ) 中抽取的样本，那么扰动后它们仍应是从该分布中抽取的样本。实现这一点的一种方式是对分布 p(h|d₁, …, dᵢ) 应用一次或多次马尔可夫链蒙特卡洛（MCMC）算法的迭代（详见 Fearnhead, 2002）。我们将在考虑一个序贯蒙特卡洛在实际应用中的例子之后，更详细地讨论这些算法。

6.5.3 一个例子：多目标跟踪推理

在 Vul 等人（2008）提出的多目标跟踪模型（第5章已总结）中，多个目标的存在带来一个问题：尽管每个目标都遵循带有高斯噪声的线性动力学，因而可以用粒子滤波器追踪，但当目标占据相似位置时，它们可能会彼此混淆。Vul 等人（2008）通过引入一个额外变量 γₜ 对此进行了建模，γₜ 是时间 t 上各目标索引的一个排列。例如，若有八个目标，初始分配 γ₁ 将是向量 (1, 2, 3, 4, 5, 6, 7, 8)。然而，若目标1和2在下一步交叉路径并彼此混淆，γ₂ 将变为 (2, 1, 3, 4, 5, 6, 7, 8)。因此，多目标跟踪不仅涉及计算目标位置的后验分布，还涉及计算 γₜ 的后验分布。

由于该模型的定义方式，如果 γₜ 已知，则目标将按照各自的线性-高斯动力学演化。这意味着推理问题简化为跟踪 γₜ 的问题。对于每个 γₜ₊₁ 的取值，基于 γₜ 及其先前时刻的位置，计算目标在 t+1 时刻位置的概率是直接的。这意味着我们为粒子滤波器设置了一个完美的框架。

Vul 等人（2008）假设每个 γₜ 相互独立——即目标索引（以及由此产生的不同目标之间的潜在混淆）在每个时间点 t 都从相同的分布 P(γₜ) 中抽取，并且与前一时刻 t−1 的情况无关。数据包含一组对每个目标 mₜ 的“运动”观测，其中包含每个目标在时间 t 的位置和速度。粒子滤波器采用一组 γ₁, …, γₜ 的样本，捕捉基于 m₁, …, mₜ 的分配后验分布，并在观测到 mₜ₊₁ 后更新这些样本，以反映 γₜ₊₁ 的后验分布。这可以通过为每个样本增补一个从先验 P(γₜ₊₁) 中抽取的 γₜ₊₁ 值，并为每个结果样本分配一个与 p(mₜ₊₁|m₁, …, mₜ, γ₁, …, γₜ₊₁) 成比例的权重来实现——该权重可通过将模型所假设的线性动力学和高斯噪声应用于由 γₜ₊₁ 指示的目标来计算。在给定 γₜ₊₁ 的条件下，每个目标的位置演化方式如同前一章讨论的卡尔曼滤波器中那样。

6.6 马尔可夫链蒙特卡洛（Markov Chain Monte Carlo）

从概率分布中生成样本的最灵活方法之一是马尔可夫链蒙特卡洛（MCMC）。该方法可用于构建任意概率分布的采样器，即使这些分布的归一化常数未知亦可使用。MCMC算法最初是为了求解统计物理学中的问题而发展起来的（Metropolis、Rosenbluth、Rosenbluth、Teller 和 Teller，1953年），如今已广泛应用于物理学、统计学、机器学习及相关领域（更多细节参见 Newman & Barkema，1999；Gilks、Richardson & Spiegelhalter，1996；Mackay，2003；Neal，1993）。与其他蒙特卡洛方法一样，MCMC算法为处理概率分布以及关于人类如何近似实现贝叶斯推理的过程假设提供了工具。在第10章中，我们还将探讨这样一种观点：MCMC算法可用于设计行为实验，使我们能够从人们的主观概率分布中进行采样。

6.6.1 马尔可夫链

顾名思义，马尔可夫链蒙特卡洛基于马尔可夫链理论——即随机变量序列，其中每个变量在给定其直接前驱变量的条件下，与所有先前变量条件独立。在马尔可夫链中，一个变量取某一特定值的概率，取决于前一变量的值，该概率由该马尔可夫链的转移核（transition kernel）决定。转移核 K(x⁽ⁱ⁺¹⁾ | x⁽ⁱ⁾) 表示马尔可夫链在第 i 步处于状态 x⁽ⁱ⁾ 时，在第 i+1 步转移到状态 x⁽ⁱ⁺¹⁾ 的概率。当状态是离散的时，这些信息可以编码在一个矩阵中，其中行对应下一次迭代的状态值，列对应由核指定的分布。这被称为转移矩阵（transition matrix）。

马尔可夫链的一个众所周知的性质是它们倾向于收敛到一个平稳分布（stationary distribution）：随着马尔可夫链长度的增加，链中某个变量取某一特定值的概率会收敛到一个固定量，该量由转移核的选择所决定。平稳分布被定义为满足以下等式的分布 π(x)：

这意味着，如果我们从 π(x′) 中采样一个状态 x′，然后根据分布 K(x|x′) 选择下一个状态 x，则状态上的概率分布仍保持为 π(x)。这就是该分布被称为“平稳”的含义：一旦状态上的概率分布达到平稳分布，它就不再改变。对平稳分布的收敛性可以通过正式定义 p(x⁽ᵐ⁾ | x⁽¹⁾) 来描述，该式表示马尔可夫链从初始状态 x⁽¹⁾ 开始，经过 m 次迭代后处于状态 x⁽ᵐ⁾ 的概率。只要马尔可夫链满足遍历性（ergodicity）条件（详见 Norris, 1997），当 m 趋于无穷大时，该概率将收敛到平稳分布。也就是说，

对于所有 x⁽¹⁾ 和 x⁽ᵐ⁾ 成立。因此，如果我们从马尔可夫链中通过选取某个初始值，然后反复根据转移核所指定的分布进行采样，最终我们将生成来自平稳分布的样本。出于同样的原因，当 m 足够大时，对函数 f(x) 在 x⁽ᵐ⁾ 取值上的平均值将逼近该函数在概率分布 π(x) 上的平均值。¹

6.6.2 Metropolis-Hastings 算法

在 MCMC 中，构造一个马尔可夫链，使其平稳分布为我们希望从中生成样本的目标分布。如果目标分布是 p(x)，则转移核 K(x⁽ⁱ⁺¹⁾ | x⁽ⁱ⁾) 将被选择为使得 π(x) = p(x) 成为平稳分布。幸运的是，存在一种简单的程序，可用于为任意给定的 p(x) 构造具有该性质的转移核。这一程序被称为 Metropolis-Hastings 算法（Hastings, 1970；Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller, 1953）。

Metropolis-Hastings 算法将核 K(x⁽ⁱ⁺¹⁾ | x⁽ⁱ⁾) 定义为两个概率步骤的结果。第一步使用一个任意的建议分布 q(x* | x⁽ⁱ⁾)，生成一个针对 x⁽ⁱ⁺¹⁾ 的候选值 x*。第二步是决定是否接受该候选值，依据是该候选值是如何生成的，以及它在目标分布下的概率。确保最终的马尔可夫链收敛到目标分布的关键在于接受函数（acceptance function）的定义方式。

Metropolis-Hastings 算法中所使用的接受函数是

其中 p(x) 是我们的目标分布。如果从 [0, 1] 上的均匀分布生成的一个随机数小于 A(x*, x⁽ⁱ⁾)，则提议值 x* 被接受为 x⁽ⁱ⁺¹⁾ 的取值。否则，马尔可夫链保持在前一个值上，即 x⁽ⁱ⁺¹⁾ = x⁽ⁱ⁾。直观地讲，公式 (6.21) 表明：若候选值 x* 在目标分布下具有较高概率，则我们更有可能接受该候选值；而若提议分布强烈偏好候选值 x* 而非当前值 x⁽ⁱ⁾，则我们接受它的可能性较低。这两个因素意味着我们的马尔可夫链会在 p(x) 值较高的区域花费更多时间，且不会过度受提议分布波动的影响。如果你好奇公式 (6.21) 的来源，我们将在下一节中详细解释。

综上所述，算法流程如下：从初始值 x⁽¹⁾ 开始。定义一个提议分布 q(x* | x⁽ⁱ⁾)，并根据 x⁽ⁱ⁾ 采样得到 x*。通过应用公式 (6.21) 并代入目标分布 p(x)，决定是否接受 x*。若接受，则 x⁽²⁾ = x*；否则，x⁽²⁾ = x⁽¹⁾。重复此过程多次，直到马尔可夫链明显收敛到其平稳分布，此时对所得 x⁽ᵐ⁾ 取平均即可近似该量在 p(x) 上的期望值。

为了说明 Metropolis-Hastings 算法的应用，我们可以用它从高斯分布中生成样本。为简化起见，我们假设高斯分布均值为零、方差为 1，但同样的程序也可用于更有趣的情形，例如从高斯均值的后验分布中采样。一般而言，从高斯分布中采样很容易，因此我们并不需要使用 MCMC，但它提供了一种方式来演示所涉及的步骤。我们马尔可夫链的状态空间是一维的，范围从 -∞ 到 ∞。我们需要选择一个提议分布，以便从一个状态转移到另一个状态。通常，最好选择一个既不太大（意味着绝大多数提议都会被拒绝）也不太小（意味着提议会被接受，但在状态空间中几乎没有任何实际移动）的提议分布。我们将使用一个高斯分布作为提议分布，其均值对应于前一个值，

且方差为目标分布方差的 0.2 倍。重复“生成提议并决定是否接受”的过程 500 步，并采用三个不同的起始点，所得结果如图 6.4 所示。

Metropolis-Hastings 算法的一个优点是，它只需要对概率分布 p(x) 有有限程度的了解。观察公式 (6.21) 可知，实际上即使我们只知道 p(x) 的值仅至一个常数倍，Metropolis-Hastings 算法仍然可以应用，因为只有这些量的比值会影响算法。这在进行贝叶斯推理时极为有用——正如我们之前所展示的，应用贝叶斯规则最困难的部分通常是计算归一化常数。使用 Metropolis-Hastings 算法，我们无需计算这个归一化常数即可从后验分布中生成样本——我们只需知道似然函数和先验分布。这使得处理概率模型变得容易得多，因为定义似然函数和先验分布只需考虑相关的生成过程和归纳偏置，而无需考虑解析可处理性的问题。因此，我们可以超越共轭先验的受限世界，开始探索更丰富、更具表达力的模型。

6.6.3 分析 Metropolis-Hastings 算法

Metropolis-Hastings 算法看起来有点像魔法——它允许我们定义一个马尔可夫链，使其收敛到任意我们想要的平稳分布。这怎么可能？关键在于由公式 (6.21) 给出的接受函数。我们可以利用马尔可夫链的另一性质推导出这一接受函数，并解释它为何有效。

假设我们有一个转移核为 K(x⁽ⁱ⁺¹⁾ | x⁽ⁱ⁾) 的马尔可夫链，它满足以下方程：

对所有 x 和 x′ 成立，其中 π(x) 是该马尔可夫链的平稳分布。这一性质被称为细致平衡（detailed balance）——直观地讲，这意味着一旦马尔可夫链达到其平稳分布，从 x 到 x′ 或从 x′ 到 x 的转移概率是相同的（具有此性质的马尔可夫链被称为可逆的）。细致平衡是一个比存在平稳分布更强的条件——你可以验证，只要对公式 (6.22) 两边同时关于 x⁽ⁱ⁾ 积分，即可得到公式 (6.19)。

Metropolis-Hastings 算法中使用的接受函数是通过细致平衡推导出来的。与之前一样，令 A(x*, x⁽ⁱ⁾) 表示在状态 x⁽ⁱ⁾ 下接受候选值 x* 的概率。每次提议被拒绝时，都会向 x⁽ⁱ⁾ 添加一小部分概率质量，因此“未接受任何提议”的概率（记作 A = 0）是

如果我们能定义一个满足该关系的接受函数，那么我们就定义了一个以 π(x) 为平稳分布的马尔可夫链。任何与 q(x⁽ⁱ⁾, x*)p(x) 成比例的接受函数 A(x*, x⁽ⁱ⁾) 都将导致一个以 p(x) 为平稳分布的马尔可夫链。可以很容易验证，公式 (6.21) 中给出的接受函数（即在 Metropolis-Hastings 算法中所使用的函数）满足这一条件，因为它最大化了提议被接受的概率（从而最小化了 P(A=0 | x⁽ⁱ⁾)）。然而，我们将在第 10 章有机会探索其他一些有效的接受函数。

6.6.4 Gibbs 采样

如果我们能够从与 p(x) 相关的分布中采样，我们就可以使用其他 MCMC 方法。特别是，如果我们能够从给定其余变量条件下每个变量的条件概率分布 p(xⱼ | x₁, …, xⱼ₋₁, xⱼ₊₁, …, x_d) 中采样，我们就可以使用另一种流行的算法——Gibbs 采样（Geman & Geman, 1984；Gilks 等, 1996），它在统计物理学中被称为热浴算法（heatbath algorithm）（Newman & Barkema, 1999）。针对目标分布 p(x) 的 Gibbs 采样器是由如下定义的马尔可夫链：对每个 xⱼ，从其条件分布 p(xⱼ | x₁, …, xⱼ₋₁, xⱼ₊₁, …, x_d) 中抽取样本。你实际上可以从 Metropolis-Hastings 算法推导出 Gibbs 采样器：如果你提议仅从其余元素取值给定的条件下对 x 的一个元素进行采样，则该提议总是会被接受。这正是 Gibbs 采样所做的。你可以依次循环遍历所有 xⱼ，或在每次迭代中随机选择一个进行采样。

Gibbs 采样与 EM 算法有很多共同之处，但它通常用于从后验分布中生成样本，而 EM 算法则产生点估计。可以通过考虑如何将 Gibbs 采样应用于高斯混合模型来说明这两个算法之间的共性，在该模型中，各分量的方差已知但均值未知（前一章中讨论过；图 6.5 为方便起见重新绘制了该图形模型）。在这种情况下，我们希望计算在给定观测数据 x 的前提下，关于观测数据分配到聚类 z、以及分量均值 μ⁽¹⁾ 和 μ⁽²⁾ 与混合比例 π 的后验分布（这里为了简化，假设我们的高斯分布仅为一维）。从 z、μ 和 π 的某些初始值出发，Gibbs 采样器会循环遍历这些变量，每次从其余所有变量给定的条件下抽取每个变量的值。这意味着该算法有效地在“采样观测数据分配至各分量”和“采样参数值”之间交替进行。

该算法的这两部分类似于 EM 算法中的 E 步和 M 步，尽管在 Gibbs 采样中，采样同时用于处理隐变量和参数，而 EM 算法则使用一个过程（期望）处理隐变量，另一个过程（最大化）处理参数。

在此混合模型中应用 Gibbs 采样非常简单。在每次迭代中，我们循环遍历所有的 zᵢ，并从其条件分布中抽取每个 zᵢ 的值。

其中 z₋ᵢ 表示 (z₁, …, zᵢ₋₁, zᵢ₊₁, …, zₙ)，I(·) 是指示函数，当其参数为真时取值 1，否则取值 0；我们利用了该模型中所包含的条件独立性假设。接着，我们在给定这些 z 值的条件下，从 μ⁽¹⁾、μ⁽²⁾ 和 π 的后验分布中抽取样本。由于在给定 z 的条件下等价于观测到 z，这变得极为简单：我们只需计算这些参数的后验分布，假设每个观测值被分配到 z 中所指示的分量。μ⁽ʲ⁾ 的后验分布即是在所有满足 zᵢ = j 的 xᵢ 条件下的后验分布（如我们之前分析所示，该分布是高斯分布），而 π 的后验分布是 Beta(n₁ + 1, n₂ + 1)，前提是假设 π 上服从均匀先验。从这两种分布中生成样本都很容易。Gibbs 采样器随后在这两个步骤之间来回迭代：先采样分配，再采样参数，然后再次采样分配，依此类推。图 6.6 展示了一个在二维高斯混合模型中使用该采样器的例子。

在推导 Gibbs 采样器时，我们需要计算每个变量在给定其余所有变量条件下的条件分布。图形模型在此背景下极为有用。请记住，图形模型所编码的一个重要内容是变量集合中存在的依赖关系模式。特别是，当我们对除目标变量外的所有其他变量进行条件化时，我们知道该变量将独立于图中既非其父节点、也非其子节点、亦非其共同子节点的父节点的任何其他变量。这一组变量被称为该变量的马尔可夫毯（Markov blanket），因为它屏蔽了该变量与其他依赖关系的联系。因此，当我们计算 Gibbs 采样中所使用的条件分布时，我们只需考虑出现在目标变量的马尔可夫毯中的变量。

你可以通过使用图 6.5 中所示的图形模型为每个变量生成马尔可夫毯，来验证这在本讨论所涉及的高斯混合模型中确实成立。

6.6.5 一个例子：主题模型的 Gibbs 采样

使用第 5 章介绍的主题模型需要解决一个困难的概率推理问题——即找出最能解释一组文档内容的主题。作为提醒，每个词 wᵢ 都关联着一个主题 zᵢ，每个主题是一个以参数 φ 表示的、在词汇上的多项分布，而文档内主题的概率分布则是一个以参数 θ 表示的多项分布。我们对 φ 和 θ 均定义了一个狄利克雷先验（即 beta 分布的多变量推广），其中 φ 从参数为 β 的对称狄利克雷分布生成，θ 从参数为 α 的对称狄利克雷分布生成。

我们可以完全以无监督的方式，利用贝叶斯推理从一组文档中提取出 t 个主题。由于狄利克雷先验与多项分布是共轭的，我们可以像第 3 章模型选择示例中所做的那样，通过对 φ 和 θ 积分，计算联合分布 P(w, z)。然后，我们可以根据贝叶斯规则，提出关于给定 w 时 z 的后验分布的问题：

由于分母中的求和项不可解析计算（包含 tⁿ 项），我们被迫使用 MCMC 来评估该后验分布。在这种情况下，我们使用 Gibbs 采样来研究词到主题的分配 z 的后验分布。

Gibbs 采样算法包括：首先为每个词随机选择一个初始主题分配（例如，对每个词均匀随机选择一个主题），然后从条件分布 P(zᵢ | z₋ᵢ, w) 中对每个词 zᵢ 的分配进行采样。因此，该算法的每一次迭代都是对词到主题分配的一次概率重排。这一过程在图 6.7 中进行了说明，该图展示了将该算法（仅使用两个主题）应用于 TASA 语料库一小部分的结果。该部分包含 30 篇使用单词 MONEY 的文档、30 篇使用单词 OIL 的文档，以及 30 篇使用单词 RIVER 的文档。词汇表被限制为 18 个词，条目表示这 731 个词例（即实例）在 90 篇文档中出现的频率。语料库中的每个词例 wᵢ 在每次采样步骤迭代中都被分配一个主题 zᵢ。在图中，我们聚焦于三个词的词例：MONEY、BANK 和 STREAM。每个词例最初被随机分配一个主题，而每次 MCMC 迭代都会产生一组新的词例到主题的分配。经过几次迭代后，主题分配开始反映出 MONEY 和 STREAM 的不同用法模式，这些词的词例最终会落入不同的主题，而 BANK 的多重含义也会被体现出来。

该特定 Gibbs 采样算法背后的细节见 Griffiths 和 Steyvers (2004)，其中该算法被用于分析大型科学文献数据库中出现的主题。算法中使用的 zᵢ 的条件分布可通过类似于我们在硬币投掷问题中推导后验预测分布的论证得出：

其中，z₋ᵢ 表示除 zᵢ 外所有 zₖ 的分配；n⁽ʷⁱ⁾₋ᵢ,ᶻᵢ 是与 wᵢ 相同且被分配给主题 zᵢ 的词的数量；n⁽·⁾₋ᵢ,ᶻᵢ 是分配给主题 zᵢ 的总词数；n⁽ᵈⁱ⁾₋ᵢ,ᶻᵢ 是来自文档 dᵢ 并被分配给主题 zᵢ 的词的数量；ν 是模型词汇表中的词数；n⁽ᵈⁱ⁾₋ᵢ· 是文档 dᵢ 中的总词数；以上所有计数均不包含当前词 wᵢ 的分配。该表达式中的两项都有直观的解释，分别对应于在给定当前分配 z₋ᵢ 的前提下，词在主题内的后验预测分布，以及主题在文档内的后验预测分布。MCMC 算法的结果是一组来自 P(z|w) 的样本，反映了给定一组文档时主题分配的后验分布。单一样本可用于评估语料库中出现的主题，或揭示词到主题的分配，如图 6.7 所示。我们还可以通过在多个样本上取平均，利用条件概率

来计算诸如词间关联强度等量值。²

尽管存在其他可与此生成模型配合使用的推理算法（例如，Blei, Ng, & Jordan, 2003；Minka & Lafferty, 2002），但 Gibbs 采样器是一种极其简单（且效率合理）的方法，用于探究使用主题表示词之间语义关系所带来的后果。Griffiths 和 Steyvers (2002, 2003) 提出，主题模型可能为传统的语义表征方法提供一种替代方案，并证明它们在预测人类词语关联数据方面比潜在语义分析（LSA）（Landauer & Dumais, 1997）表现更好。图 6.8 是对该工作早期介绍的一个同时代记录。主题模型还可应用于一系列其他依赖语义关联的任务，例如语义启动和句子理解（Griffiths 等, 2007）。

6.7 变分推理

蒙特卡洛方法并非用于推理的唯一一类近似方法——另一种常用的方法是变分推理（Jordan 等，1999）。蒙特卡洛方法通过使用随机样本近似真实分布来处理计算复杂性，而变分方法则通过在某个更简单的分布族中寻找与真实分布最接近的成员来近似该分布。因此，与本章前面讨论的无偏但带有噪声的蒙特卡洛方法不同，变分方法提供了一种可能有偏但无噪声的近似。

寻找与真实分布最接近的简单分布族成员这一过程，正是变分方法名称的由来：它们基于变分法（calculus of variations），而变分法关注的是将函数（即真实分布和近似分布）映射到标量（即两个分布之间的偏差）的问题。由于近似分布比真实分布更简单，通常可以解析地计算出该近似。在无法解析求解的情况下，变分近似仍可表示为一种形式，能够借助近似算法（如随机梯度下降）进行求解。正因如此，随着深度神经网络训练技术的进步，出现了用于求解这类优化问题的新方法和软件，变分方法也因此变得越来越流行（参见第12章）。

我们将通过第5章讨论过的EM算法来介绍变分推断。复杂隐变量模型的参数估计正是推动变分方法发展的动因，也为理解其背后的形式体系提供了良好的直觉基础。关于更详细的技术性介绍，参见Blei等人（2017）。

6.7.1 变分EM算法

在第5章中，我们介绍了EM算法作为一种估计隐变量模型参数θ的方法，该方法无需对隐变量z进行求和运算。在该算法的E步（期望步）中，我们计算完整数据对数似然的期望：

在我们所介绍的例子中（例如混合模型），各个隐变量 zᵢ 在条件上彼此独立。我们能够利用这种条件独立性高效地计算期望完整对数似然。然而，在其他模型中，当隐变量高度相互依赖时，计算期望完整对数似然可能变得不可解析，因为它仍需对所有 z 的取值进行求和。变分EM算法是解决这一问题的一种策略（参见 Jordan 等人，1999）。在该算法中，E步通过用一个更易处理的分布族中最接近的函数替代 P(z|x, θ) 来完成。

定义一个关于 z 的任意分布 Q(z)，它在 p(x, z|θ) 非零的所有位置上也非零。然后，在公式 (6.35) 中，取关于 Q(z) 而非 P(z|x, θ) 的期望。这给出：

这意味着，关于 Q(z) 的完整对数似然的期望同时也是 log p(x|θ) 的一个下界，因为 KL 散度和熵均为非负值。

这一结果为我们提供了一种不同的方式来处理 EM 算法：我们可以从一个便于计算期望的分布族中选择 Q(z)。取 Q(z) = P(z|x, θ) 可能会得到更紧致的下界，因为它最小化了 KL 散度项，但也可以使用其他分布。特别是，我们可以选择其中 z 是（条件）独立的分布。这种通过定义一个分布 Q(z) 以使推断变得可处理，并随后最小化该分布与 P(z|x, θ) 之间 KL 散度的方法，正是变分方法的基础。

6.7.2 将推理视为优化

尽管变分方法最初是在 EM 算法的背景下发展起来的，但它们为贝叶斯推断问题提供了更通用的解决方案。即使我们并不试图估计隐变量模型的参数，一个接近后验分布 P(z|x) 的分布 Q(z) 也可以更一般地用作后验分布的近似。为了与本节后续内容强调近似后验分布保持一致，我们将省略对 θ 的依赖关系。

不幸的是，虽然寻找使得 D_KL[Q(z) || P(z|x)] 较小的 Q(z) 在直觉上是合理的，但这并非一种实用的策略。问题在于，按照公式 (6.38) 所定义的方式计算 KL 散度仍需要我们评估 P(z|x)。利用公式 (6.39) 给出的形式，我们有

6.7.3 选择近似分布族

通常需要数学推导来确定如何高效地从分布族中选出最接近真实后验分布的成员。幸运的是，如果先验分布属于指数族分布（包括高斯分布、指数分布、伽马分布和贝塔分布，如第3章所述），且先验与似然函数是共轭的（同样在第3章讨论过），那么后验分布也将属于指数族，此时可通过一个简单的迭代算法找到变分近似（详见 Blei 等人，2017）。

尽管变分方法通常比蒙特卡洛方法在计算上更简单，但一个关键的选择在于近似分布族的选取。必须在近似的准确性（即它能在多大程度上模拟真实分布）与可处理性（即计算难度）之间取得平衡。由于分布的大部分复杂性源于变量之间的依赖关系，通过假设变量间完全独立，可以避免大量复杂性——这就是平均场近似（mean-field approximation）（Jordan 等人，1999）。然而，完全独立的假设往往过于强烈，无法提供有用的近似，特别是对于那些预期具有高度相关性的变量而言。例如，将你对“一个孩子是否住在蓝色房子里”的信念与你对“她的妹妹是否住在蓝色房子里”的信念近似为相互独立，就不会非常准确。

一种折中方案是仅在对计算可处理性影响最大的地方假设独立性，这种方法被称为结构化平均场近似（structured mean-field approximation）。在这种情况下，变量被划分为若干集合 zℓ，而近似分布是

6.7.4 变分推断的替代方法

变分推断的一个缺点是，它所最小化的 KL 散度 D_KL[Q(z) || P(z|x)] 对 Q(z) 施加了一些强约束。如果你查阅公式 (6.38)，你会发现：当 P(z|x) 为零时，Q(z) 必须被强制设为零；否则，KL 散度将趋于无穷大。这可能导致 Q(z) 过于狭窄，从而相较于 P(z|x) 表现出过度自信（Minka, 2005）。

另一种替代方法是使用一个在相反方向上最小化 KL 散度的近似，即 D_KL[P(z|x) || Q(z)]。这最初是在假设密度滤波（ADF）算法中实现的，该算法旨在以类似于粒子滤波器的方式确定性地过滤序列数据（Boyen & Koller, 1998）。使用这一版本的 KL 散度会导致结构化平均场近似中产生不同的最优分布成员。在此定义下，对某个划分的最佳近似就是后验分布的边缘分布：

因此，后验分布的近似就是其每个划分的边缘分布的乘积。将 ADF 推广到非序列数据的情形被称为期望传播（expectation propagation）（Minka, 2001）。然而，使用这种方向上的 KL 散度来寻找近似分布更具挑战性，并且这两种方法均未与像 ELBO 这样的全局界相联系。

6.7.5 一个例子：局部贝叶斯学习

ADF 顺序处理输入数据，这意味着即使在一个本应不受数据到达顺序影响的概率模型中，数据到达的顺序仍会影响近似后验分布。这类顺序对推断的影响被称为“顺序效应”（order effects）。作为一个简单示例，考虑图 6.9 中所示的一对二元变量联合分布的更新。先验初始时是独立的，它将根据一次具有独立似然的观测和一次具有依赖似然的观测进行更新（例如，变量之间存在强相关性）。

此处，我们展示了观测顺序对最终后验分布的影响。当依赖似然的观测先出现时，ADF 后验分布与先验相同，因此在两次观测之后，近似后验分布并未体现出第一次观测的影响。重要的是，即使 ADF 后验的边缘分布（按列求和）也不再与精确后验的边缘分布匹配，这是由于该近似被迭代应用所致。然而，当独立观测先出现时，第一次后验分布对于 ADF 和精确推断是相同的；在两次观测之后，ADF 后验体现了两次观测的影响。虽然精确推断得到的最终后验分布在这两种顺序下是相同的，但 ADF 在完整联合分布和边缘分布上都表现出顺序效应。

ADF 的优势在于，一个独立的后验分布可以仅由其边缘分布表示，从而省去了表示依赖关系的需要。就本例而言，这种节省微不足道：表示边缘分布只需两个数字，而表示联合分布只需三个数字。但随着变量数量的增长，这种节省变得显著：对于 n 个二元变量，仅需 n 个数字即可表示边缘分布，而表示联合分布则需要 2ⁿ - 1 个数字。然而，仅用边缘分布表示是一个吸引人的想法，尤其对那些试图用近似贝叶斯认知模型逆向工程人类心智的研究者而言。这种近似的优点和缺点都很有吸引力：表征复杂性的降低有利于心理合理性，而引入的顺序效应则可与人类参与者身上观察到的现象相比较。

尽管在全局贝叶斯学习（GBL）中，观测顺序不影响后验分布的精确推断，但如果在更新后验分布时假设

之间存在有限通信（如局部贝叶斯学习（LBL）所假设），则会产生顺序效应。具体而言，LBL可以产生已在关联学习实验中观察到的复杂顺序效应：回溯性再评估和突出显示，这些效应结合了某些习得线索组合的首因效应与其他近因效应（Kruschke, 2006）。虽然在假设世界正在变化的前提下，近因效应可能被纳入模型，因为较新的观测值更有价值，因此应赋予更大的权重（如第5章介绍的贝叶斯-卡尔曼滤波器），但首因效应更难归因于模型本身，因此它们常被视为近似模型的证据。进一步的研究表明，这些顺序效应并不依赖于LBL中使用消息的确切形式，而是采用后验分布的ADF近似（称为因子化贝叶斯学习FBL），也能产生相同的人类般的顺序效应（Daw, Courville, & Dayan, 2008; Sanborn & Silva, 2013）。

6.8 总结

在许多方面，近似推断算法是贝叶斯建模的隐藏秘密，因为在几乎所有有趣的概率模型中，精确推断都是不可能的。理解不同的推断算法及其适用情境是研究这些模型的研究人员工具包的重要组成部分。通常，推断的便捷性是驱动模型选择的一个因素：写出一个极其复杂的生成模型很容易，但无法保证能够应用贝叶斯规则来计算该模型中的后验分布。最近的创新，例如高效概率编程语言的发展（参见第18章），可以使处理表达力强的概率模型变得更加容易。使用这些语言的好处之一是，用于指定概率模型的代码可以自动应用于该模型的不同推断算法。然而，基本了解这些算法是什么以及它们如何工作，对于理解它们何时会失败至关重要。

近似推断算法还为人类心智和大脑如何可能执行贝叶斯推断这一问题提供了答案。正如研究人员在贝叶斯模型中有多种算法可选以近似不同假设的后验概率一样，我们的心智和大脑也可能有许多可能的策略可用于近似贝叶斯推断。

我们将在本书后面重新探讨这一主题。第11章考虑蒙特卡洛方法如何解释人类概率推理的某些方面，第13章通过考虑心智和大脑如何最有效地利用此类计算资源来扩展这一想法，而第12章则将本章介绍的一些算法与经典和近期关于人工神经网络的工作联系起来。

原文链接：https://dokumen.pub/bayesian-models-of-cognition-reverse-engineering-the-mind-1nbsped-9780262049412-9780262381048-9780262381055.html