一文读懂贝叶斯推理问题：MCMC方法和变分推断

全文共6415字，预计学习时长20分钟或更长

图片来源：pexels.com/@lum3n-com-44775

贝叶斯推理（Bayesian inference）是统计学中的一个重要问题，也是许多机器学习方法中经常遇到的问题。例如，用于分类的高斯混合模型或用于主题建模的潜在狄利克雷分配（Latent Dirichlet Allocation，简称LDA）模型等概率图模型都需要在拟合数据时解决这一问题。

同时，由于模型设置（假设、维度……）不同，贝叶斯推理问题有时会很难解决。在解决大型问题时，精确的方案往往需要繁重的计算，要完成这些难以处理的计算，必须采用一些近似技术，并构建快速且有可扩展性的系统。

本文将讨论两种可用于解决贝叶斯推理问题的主要方法：基于采样的马尔可夫链蒙特卡罗(Markov Chain Monte Carlo，简称MCMC)方法和基于近似的变分推理(Variational Inference，简称VI)方法。

本文第一部分将讨论贝叶斯推理问题，并介绍几个机器学习应用的经典案例，当然，这些案例中会出现贝叶斯推理问题。第二部分将全面介绍用于解决该问题的MCMC技术，并详细介绍其中的两种算法：Metropolis-Hasting算法和吉布斯采样（Gibbs Sampling）算法。最后，第三部分将介绍变分推断，并了解如何通过优化参数化数族分布得到近似解。

注意，以a(∞)为标记的小节数学专业性非常强，跳过也不会影响对本文的整体理解。还要注意，本文中的p(.)可以用来表示概率、概率密度或概率分布，具体含义取决于上下文。

贝叶斯推理问题

这一部分提出了贝叶斯推理问题，讨论了一些计算困难，并给出了LDA算法的例子。LDA算法是一种具体的主题建模机器学习技术，能够反映贝叶斯推理问题。

统计推断旨在根据可观察到的事物来了解不可观察到的事物。即，统计推断是基于一个总体或一些样本中的某些观察变量（通常是影响）得出结论的过程，例如关于总体或样本中某些潜在变量（通常是原因）的准时估计、置信区间或区间估计等。

而贝叶斯推理则是从贝叶斯的角度产生统计推断的过程。简而言之，贝叶斯范式是一种统计/概率范式，在这种范式中，每次记录新的观测数据时就会更新由概率分布建模的先验知识，观测数据的不确定性则由另一个概率分布建模。支配贝叶斯范式的整个思想嵌入在所谓的贝叶斯定理中，该定理表达了更新知识（“后验”）、已知知识（“先验”）以及来自观察的知识（“可能性”）之间的关系。

一个经典的例子是用贝叶斯推理进行参数估计。假设一个模型中数据x是根据未知参数θ的概率分布生成的，并且有关于参数θ的先验知识，可以用概率分布p(θ)来表示。那么，当观察到数据x时，我们可以使用贝叶斯定理更新关于该参数的先验知识，如下所示：

贝叶斯定理应用于给定观测数据的参数推断的说明。

计算困难

根据贝叶斯定理，后验分布的计算需要三个条件：先验分布、可能性和证据。前两个条件很容易理解，因为它们是假设模型的一部分（在许多情况下，先验分布和可能性是显而易见的）。然而，第三个条件，即归一化因子，需要如下计算：

虽然在低维中，这个积分可以较容易地计算出来，但在高维中它会变得难以处理。在上述案例中，对后验分布进行精确计算是不可行的，必须使用一些近似技术（例如平均计算）来获得后验分布。

贝叶斯推理问题还可能会产生一些其他的计算困难。例如，当某些变量是离散的时候会产生组合学问题。马尔可夫链蒙特卡罗（Markov Chain Monte Carlo，简称MCMC）和变分推理（Variational Inference，简称VI）是最常用于解决这些问题的两种方法。下文将描述这两种方法，尤其关注“归一化因子问题”，但是应该记住，这些方法也可用于与贝叶斯推理相关的其他计算困难。

为了让接下来的章节更易于理解，可以观察到，由于x应该是给定的，因此可以作为参数，那么，θ的概率分布则被定义为归一化因子

在描述MCMC和VI两个部分之前，先来看一个具体例子，了解在机器学习LDA中存在的贝叶斯推理问题。

举例

贝叶斯推理问题通常出现在需要假设概率图模型或根据给定观测值得出模型潜变量的机器学习方法中。在主题建模中，潜在狄利克雷分配(LDA)定义了一个用于描述语料库文本的模型。因此，给定大小为V的完整语料库词汇表和给定数量为T的主题，模型假设：

· 对于每个主题，在词汇表上都存在一个“主题词”的概率分布（使用Dirichlet先验假设）

· 对于每个文档，在主题上都存在一个“文档主题”的概率分布（使用另一个Dirichlet先验假设）

· 对文档中的每个单词进行采样。首先，从文档的“文档 - 主题”分布中对主题进行采样；其次，从附加到采样话题的“主题 - 单词”分布中采样一个单词。

该方法的名称来源于模型中假设的Dirichlet先验，其目的是推断观察到的语料库中的潜在主题以及每个文档的主题分解。即使不深入研究LDA方法的细节，也可以粗略地用w来表示语料库中单词的向量，用z来表示与这些单词相关的主题向量，用贝叶斯方法根据观测到的w推断出z：

由于维度过高，这里无法推断出归一化因子，同时，还存在组合问题（因为一些变量是离散的），需要使用MCMC方法或VI方法来获得近似解。对主题建模及其特定的贝叶斯推理问题感兴趣的读者可以看看下面这篇关于LDA的参考文献。

传送门：http://www.jmlr.org/papers/vo...

LDA方法的说明。

马尔可夫链蒙特卡洛（MCMC）方法

上文提到，贝叶斯推理问题中的主要困难来自于归一化因子。本节将描述MCMC采样方法，为归一化因子以及与贝叶斯推理相关的其他计算困难提供解决方案。

采样方法

采样方法如下，首先假设有一种方法（MCMC）可以从由一个因子定义的概率分布中抽取样本。然后，可以从这个分布中得到样本（仅使用未标准化的部分定义），并使用这些样本计算各种准时统计量，如均值和方差，甚至通过核密度估计来求得近似分布，从而避免处理涉及后验的棘手计算。

与下一节所述的VI方法相反，对所研究的概率分布（贝叶斯推理中的后验分布）MCMC方法无需假设模型。因此，该方法具有低偏差但高方差，这意味着大多数情况下，获得的结果比从VI方法中得到的结果花费更多时间精力，但也更准确。

总结本小节，即上述的采样过程并不局限于后验分布的贝叶斯推理，它还可以普遍用于所有由归一化因子定义的概率分布。

采样方法（MCMC）的说明。

MCMC方法的概念

在统计学中，马尔可夫链蒙特卡罗（MCMC）算法旨在从给定的概率分布中生成样本。该方法名称中的“蒙特卡罗”部分是出于取样目的，而“马尔可夫链”部分来自获取这些样本的方式。

为了得到样本，要建立一个马尔可夫链，从其平稳分布中获得样本。然后，可以从马尔可夫链中模拟随机的状态序列，该序列足够长，能够（几乎）达到稳态，再保留生成的一些状态作为样本。

在随机变量生成技术中，MCMC是一种相当高级的方法，可以从一个非常困难的概率分布中获得样本，这个概率分布可能仅由一个乘法常数定义。更出乎意料的是，可以用MCMC从一个未经标准化的分布中获得样本，这来自于定义马尔可夫链的特定方式，马尔可夫链对这些归一化因子并不敏感。

MCMC方法旨在从一个困难的概率分布中生成样本，该概率分布可以仅由一个因子定义而成。

马尔可夫链的定义

整个MCMC方法是基于马尔可夫链的建立，并从其平稳分布中取样。为此，Metropolis-Hasting和吉布斯采样算法都使用了马氏链的一个特殊性质：可逆性。

状态空间为E的马尔可夫链，转移概率由下式表示

如果存在概率分布γ，上式则是可逆的

对于这样的马氏链，可以很容易地证明有

然后，γ是一个平稳分布（对不可约马氏链来说，也是唯一一个平稳分布）。

现在假设想要采样的概率分布π仅由一个因子定义

（其中C是未知的乘法常数）。可以注意到以下等式成立

接着，是转移概率为k(.,.)的马尔可夫链被定义为验证过去的等式，如预期那样将π定义为平稳分布。因此，我们可以定义一个马尔可夫链的平稳概率分布为π，该分布不能精确计算。

Gibbs采样转换（∞）

假设待定义的Markov链是D维的，则

吉布斯采样（Gibbs Sampling）假设即使在无法得知联合概率的情况下，也可以基于其他维度计算得出某一维度的条件分布。基于此假设，Gibbs采样转换可定义为，下一阶段状态，如在n+1次迭代的状态，可由如下步骤得出。

首先，从D维X_n中随机选择一个整数d。然后，根据相应的条件概率，通过采样赋予维度d一个新数值。这一过程中，其他维度保持如下状态不变：

其中

是基于其他维度得出的第d个维度的条件分布。

通常，设

则转换概率可以表示为

并且，在唯一有意义的情况下，局部平衡按预期得到了验证

Metropolis-Hasting转换（∞）

有时候，计算Gibbs采样中的条件分布也是很复杂的。在这种情况下，可以采用Metropolis-Hasting算法。运用该算法，需要先定义一个侧向的转换概率h(.,.)，该概率将被用于建议转换。下一阶段（n+1次迭代）Markov链的状态可由如下步骤得出。首先，从h中生成“建议转换”x，并计算一个关联概率r用于接受x：

可以得到如下有效转换

通常，转换概率可以表示为

同时，局部平衡按预期得到了验证

采样过程

定义Markov链后，模拟一串随机状态序列（随机初始化数值），并对其中一些状态进行设定，如设置为服从目标分布的独立样本。

第一步，为了让样本（近似）服从目标分布，仅考虑与初始设定序列状态相差大的状态，使Markov链近似达到稳定状态（理论上来说，渐进达到稳定状态）。这样一来，初始设定状态就没样本那么有用了。这一达到平稳的阶段被称为老化时间（burn-in time）。需要注意的是，实际操作中很难知道该阶段会持续多长时间。

第二步，为了获得（近似）独立样本，不能把所有的序列连续状态都放在老化时间之后。实际上，Markov链的定义中就已经表明了两个连续状态之间有很强的联系。因此，需要把状态相差很远的样本默认为近似独立。在实际操作中，可以通过分析自相关函数来预测两个近似独立状态间所需要的滞后（仅限于数值数据）。

所以，为了得到服从目标分布的独立样本，需要从位于老化时间B之后的、彼此间滞后为L的初始序列中分离出状态。设Markov链连续状态为

则样本状态为

MCMC采样需要考虑老化时间和滞后。

变分推断（VI）

另一个可用于解决复杂推断计算问题的方法是变分推断（Variational Inference，简称VI）。VI旨在找到参数化数族的最优近似分布。为此，需要遵循一个优化过程（优化数族里的参数），该过程需要仅由一个因子定义的目标分布。

逼近法

给定一个数族，VI旨在搜寻该数族中某些复杂目标概率分布的最优近似解。具体来说，VI定义一个参数化数族分布，并通过优化参数得到具有确定误差测量的最接近目标的元素。

将归一化因子C的概率分布π定义为：

应用数学术语，设参数化数族分布为

对于两个分布p和q的误差测量E(p,q)，搜寻如下最优参数

如果想要在未明确标准化π的情况下解决该问题，那么不需要复杂的计算，f_