【摘 要】 统计推断是贝叶斯概率框架中最为重要的部分，也是概率机器学习的核心部分。几乎所有的概率机器学习模型，都会涉及根据观测量来获取隐变量或模型参数相关知识的问题，这就是统计推断。与频率主义不同，贝叶斯推断方法并不给出隐变量的确切值，而是保留了模型的不确定性，给出隐变量的概率分布。由于输出的不再是点，而是一个分布，导致贝叶斯统计推断的难度大大增加了。尤其是在复杂模型和大数据集中，问题更明显。本文对贝叶斯统计推断技术进行了综述，以便快速对相关领域知识有一个理解。

【原 文】 参考 Blei 的讲座自行整理。

1 简介

贝叶斯推断是统计学中的一个基本问题，也是许多机器学习方法中遇到的问题。例如：用于分类的高斯混合模型、用于主题建模的潜狄利克雷分配模型等概率图模型，都在拟合数据时需要解决贝叶斯推断问题。

同时应注意到，根据模型设置（如：假设、维度等）不同，贝叶斯推断问题有时很难解决。尤其是遇到大型问题中，精确推断方法往往需要大量计算，且变得难以处理，实践中通常会使用一些近似技术来克服此问题，并建立快速和可扩展的系统。

本文简单介绍贝叶斯推断存在的现实性问题，并以主题建模领域常用的潜狄利克雷分配模型为例进行讨论。

2 什么是贝叶斯推断问题？

统计推断指根据可观测到的事物来了解无法被观测到的事物的一类统计学方法。

  • 统计推断是某种估计过程，它根据总体（或总体样本）中的一些可观测变量（可观测的结果、效应等），得出关于总体中一些隐变量（不可观测的原因、模型参数等）的某些估计。

  • 贝叶斯统计方法中，常见的统计推断任务是对隐变量的 点估计、可信区间估计、分布估计 等。

  • 频率派统计方法中，常见的统计推断任务是对隐变量的 点估计 和 置信区间检验。

  • 统计推断是很多后续任务的基础，如：预测、积分等。

贝叶斯统计推断是从贝叶斯角度产生统计推断的过程。简而言之，贝叶斯范式是一种统计/概率范式，其中：由某个概率分布建模的先验知识，在每次记录由另一概率分布建模的新观测时，其不确定性被更新形成新的知识。

贝叶斯范式的整体思想可以通过贝叶斯定理完整表达，该定理表达了更新后的知识（“后验”）、先验知识（“先验”）和来自观测的知识（“似然”）三者之间的关系。

最常见的贝叶斯推断案例是 参数推断（因为模型的参数通常是需要我们求解的一种隐变量）：

  • 假设一个模型，其数据 x\mathbf{x} 从依赖于未知参数 θ\boldsymbol{\theta} 的某个概率分布中生成

  • 同时假设我们对参数 θ\boldsymbol{\theta} 已有了一些先验知识，并将其表示为概率分布 p(θ)p(\boldsymbol{\theta})

则当我们观测到数据 x\mathbf{x} 时，就可以使用贝叶斯定理来更新关于该参数的先验知识（如 图 1 所示），而更新后的知识被称为后验。

图 1 ： 贝叶斯定理在参数推断中的应用


3 贝叶斯推断的计算难题

贝叶斯定理说明，后验的计算需要三个条件：先验（Prior）、似然（Likelihood） 和 证据（Evidence，也称边缘似然）。

前两者通常可以很容易表达出来，因为它们通常是模型假设的一部分（在许多情况下，先验和似然都是通过建模明确知道的）。然而，第三项（可理解为一个归一化因子）需要在整个先验区间内求积分：

p(x)=∫θp(x∣θ)p(θ)dθp(\mathbf{x})=\int_{\boldsymbol{\theta}} p(\mathbf{x} \mid \boldsymbol{\theta}) p(\boldsymbol{\theta}) d \boldsymbol{\theta}

该积分在低维时也许能够被计算出来，但在高维时却会变得很难处理。而实际上，在高维情况下，根本无法精确计算后验分布，必须使用一些近似技术来估计后验分布，进而解决很多后续问题（如利用后验分布计算某些量的均值。一个量的均值是在某个概率分布下量值的加权平均，有些预测任务会依据后验分布对待预测的变量做加权平均以得出最终想要的结果）。此外，贝叶斯推断问题可能会引起其他一些计算困难，例如：当某些变量为离散型变量时的组合运算问题。

可以注意到， x\mathbf{x} 为给定的样本，也可被视为条件参数，一般化地，我们将面临在 θ\boldsymbol{\theta} 上定义一个概率分布（受归一化因子约束）的情况。由于归一化因子在问题域内为常数，因此有：

πx(θ)≡p(θ∣x)∝p(x∣θ)p(θ)≡gx(θ)\pi_{\mathbf{x}}(\boldsymbol{\theta}) \equiv p(\boldsymbol{\theta} \mid \mathbf{x}) \propto p(\mathbf{x} \mid \boldsymbol{\theta}) p(\boldsymbol{\theta}) \equiv g_{\mathbf{x}}(\boldsymbol{\theta})

上式中 πx(θ)\pi_{\mathbf{x}}(\boldsymbol{\theta}) 代表后验分布， gx(θ)\rm{g}_{\mathbf{x}}(\boldsymbol{\theta}) 代表贝叶斯定理的分子项。

迄今为止，克服归一化因子计算难题的主要近似统计推断方法是： 马尔可夫链蒙特卡洛方法（MCMC） 和 变分推断法（VI）。 前者通过科学可靠地抽取后验分布样本实现对后验分布的模拟，而后者则通过引入一个可参数化的变分分布族来近似后验分布，并且两者都巧妙地规避了计算归一化因子。

4 问题展示

贝叶斯推断问题自然而然会出现在基于概率图模型的机器学习方法中：「通过给定的一些观测，我们希望恢复模型中造成该观测的某些隐变量（通常具有一些不确定性）」。

举例来说，在 主题建模（Topic Modelling） 问题中，潜狄利克雷分配（Latent Dirichlet Allocation, LDA） 就是为语料库中文本的主题分析而定义的一个模型。

4.1 假设与定义

对于文档数为 DD 、单词数为 VV 、主题数为 TT 的一个语料库，LDA 模型 假定：

  • 语料库： 是一个由 MM 个文档构成的集合，记为 $$C={d_1,d_2,…,d_M}$$

  • 文档： 由具有 NN 个单词的词序列 $$d={w_1,w_2,…,w_N}$$ 表示，其中 wnw_n 为序列中第 nn 个单词。

  • 单词： 是用离散数值表示的最基础的语料单位。每个单词被定义为由 $${1,…,V}$$ 索引的词汇表中的一项。我们在这里使用 独热向量（One Hot Vector） 来做单词的向量化表示。


独热向量（One Hot Vector） 定义了一个长度为 VV 的向量，将该向量中仅仅用于索引某个单词的分量被置为 11，所有其他分量都被置为 00 ，这样每个独热向量值将仅代表一个单词。

也就是说，当使用上标来索引某个分量时，该词汇表中的第 vv 个单词将由一个长度为 VV 的向量 w{w} 表示，而该向量中只有第 vv 个分量 wv=1{w}^v=1 ，而其他单词 u≠vu \neq v 对应的分量 wu=0{w}^u=0。


  • 主题及其单词概率分布： 对于每个主题 tt，存在一个该主题下的单词概率分布 φt(w)\varphi_{t}(w)，代表了每个单词与该主体的关联强弱。整个统计推断过程假设 φt\varphi_{t} 的先验分布是一个参数为 β\beta 的狄利克雷分布。

  • 文档及其主题概率分布： 对于每个文档 dd，存在一个该文档下各种主题的概率分布 θd\boldsymbol{\theta}_d，代表了每个主体与该文档的关联强弱。整个统计推断过程假设其先验为参数为 α\alpha 的另一个狄利克雷分布。


在案例中，可观测变量是文档（以及其中的单词序列）和其对应的主题，不可观测的隐变量是两个分布的参数 α\alpha 和 β\beta。我们的诉求是通过对语料库的学习，能够建立一个文档（单词序列）与主题之间的概率模型，以便未来对某个给定文档的主题作出预测（通过该文档中的单词概率分布）。

4.2 生成过程及建模

假设文档 dd 中的所有单词 wd,nw_{d,n} 均来自于如下生成过程 （ n∈{1..N}n \in \{1..N\} ）：

  • 首先从文档 dd 的主题概率分布 θd\boldsymbol{\theta}_d 中采样得到一个主题 tt（即随机生成一个主题）

  • 然后从主题 tt 对应的单词概率分布 φt\varphi_{t} 中采样得到一个单词（即根据主题对应的单词概率分布，随机生成一个单词）

狄利克雷分配模型的名称来源于其采用了狄利克雷分布作为模型参数的先验，其目的是根据语料库中所有文档数据，求得一个有关于文档主题的概率模型，并为新文档预测主题。

即使不深入研究 LDA 的细节，也可以大致理解，对于语料库中的单词序列 （在一般化形式中表示为向量 w\mathbf{w}） 和与其相关联的主题（在一般化形式中表示为 z\mathbf{z} ），我们希望以贝叶斯方式根据观测到的单词序列 w\mathbf{w} 来预测其主题 z\mathbf{z} ：

p(z∣w)=p(w∣z)p(z)p(w)=p(w∣z)p(z)∫zp(w∣z)p(z)dzp(\mathbf{z} \mid \mathbf{w})=\frac{p(\mathbf{w} \mid \mathbf{z}) p(\mathbf{z})}{p(\mathbf{w})}=\frac{p(\mathbf{w} \mid \mathbf{z}) p(\mathbf{z})}{\int_{\mathbf{z}} p(\mathbf{w} \mid \mathbf{z}) p(\mathbf{z}) d \mathbf{z}}

4.3 面临的难题

上述案例中，后验分布中的归一化因子 p(w)p(\mathbf{w}) 由于维数巨大而难以处理。此外所有变量（ w\mathbf{w} 和 z\mathbf{z} ） 都是离散型的，后验计算将面临排列组合的挑战。

注：对主题建模及其底层贝叶斯推断问题感兴趣的读者可以看看这篇关于 LDA 的参考论文。

图 2： 潜狄利克雷分配的概率图模型图解

5 MCMC 方法

正如之前提到的，处理贝叶斯推断问题时面临的主要困难之一来自归一化因子，并且在高维度时几乎没有解析解。本节描述了一些 MCMC 采样方法，它们为解决该问题以及其他与贝叶斯推断相关的计算难题提供了一种解决方案。

5.1 什么是采样方法？

采样方法的思路如下：

假设存在一种可以从某个概率分布（既然是概率分布，那一定严格定义在某个归一化因子之上）中抽取样本的方法，则我们可以从该分布中抽取样本（有可能不需要计算归一化因子的情况下），并使用这些抽取到的样本来模拟基于该分布的各种计算（例如：计算均值、方差等 ），也可以进一步通过核密度估计方法在样本基础上获得近似的完整分布。

与下一节描述的变分推断方法不同，采样方法不假定所研究概率分布的具体形态，因此适用性较广，且该方法的偏差（bias）很小，而方差（variance）很大，意味着在大多数情况下，其获得的结果更准确，但所需成本比变分推断更昂贵。

此外，上述采样方法并不局限于后验分布的贝叶斯推断，也适用于任何归一化因子计算比较困难的概率分布。

图 3 面向一维概率分布的采样方法图解

5.2 MCMC 思想

采样方法的关键步骤在于如何根据概率分布抽取合适的随机样本，使不同值域内的样本出现频率与其概率成正比。在统计方法中，马尔科夫蒙特卡罗方法就是生成此类样本的一类高效率方法。

  • 采样与最优化不同：采样涉及从整个分布中抽取样本，得到的是样本序列；而最优化只是寻找整个分布中的最值位置，得到的是数据点。
  • 实现同比例抽取样本主要有两种方式：
    • 一是所有区域按照与概率同比例的抽取密度抽取样本（接受拒绝采样的主要思路）；
    • 二是所有区域采取同样的抽样密度抽取样本，但给抽取的样本赋予与概率同比例的系数（重要性采样的主要思路）

为了生成样本，MCMC 方法的基本想法是：建立一条马尔可夫链，并从该马尔可夫链的 平稳分布 中抽取我们想要的样本。在具体实施层面，可以模拟生成马尔可夫链的一个长度足以达到平稳状态的随机状态序列，然后保留其中一些生成状态作为样本。

➀ 什么是平稳分布？

当某个概率分布 {πj,j∈S}\{\pi_j,j \in S\} 满足如下等式时，则称其是齐次马氏链 XX 的一个平稳分布（ SS 代表状态空间 ）：

πj=∑i∈Sπipij,j∈S\pi_j = \sum_{i\in S} \pi_i p_{ij}, j \in S

或采用矩阵形式满足下式：

π=πP\pi = \pi P

其中平稳分布 $$\pi = {\pi_1,\pi_2, \ldots }$$ ， P=(pij)P = (p_{ij}) 是马尔可夫链 XX 的转移概率矩阵。

➁ 平稳分布的几个性质

若 π\pi 是平稳分布，则一条马尔可夫链在经过 nn 次状态转移后，依然满足：

（1）如果状态的初始分布为 $${\pi_i}$$，其中 $$\pi_i= P(X_0 = i),i \in S$$ ，则在 nn 次状态转移后的分布不变：

P(Xn=i)=πi,i∈SP(X_n = i) = \pi_i, i \in S

或写成矩阵形式满足

π=πPn\pi = \pi P^n

即马氏链的绝对分布是确定的、保持不变的。

(2) 对于任意正整数 nn 和 mm， 以及任意时刻 0≤t1<t2<…<tn0 \leq t_1 \lt t_2 \lt \ldots \lt t_n 和状态值 i1,i2,…,in∈Si_1,i_2,\ldots, i_n \in S ，有联合概率的时间不变性：

P(Xt1+m=i1,Xt2+m=i2,…,Xtn+m=in)=P(Xt1=i1,Xt2=i2,…,Xtn=in)P(X_{t_1+m} = i_1, X_{t_2+m}=i_2,\ldots,X_{t_n+m}=i_n) = P(X_{t_1} = i_1, X_{t_2}=i_2,\ldots,X_{t_n}=i_n)

即该马氏链是一个严平稳时间序列。

在所有的随机变量生成技术中，MCMC 是一种非常先进的方法，它可以从由定义在归一化因子上的复杂概率分布中获取样本。最重要的事情时： MCMC 方法 并不要求该分布已经做过归一化处理，也就是说马尔可夫链对归一化因子不敏感。

图 4： MCMC 方法的目的是从复杂概率分布中产生样本

5.3 马尔科夫链的定义

整个 MCMC 方法基于构建马尔可夫链的能力，而该马尔可夫链的平稳分布就是我们想要的样本。为此，两种最基础的 MCMC 采样算法（ Metropolis-Hasting 采样算法 和 Gibbs 采样算法）都使用了马尔可夫链的一个特殊性质：可逆性（ 等效于满足细致平衡条件 ）。

状态空间指所有可能的状态值构成的空间。

在状态空间 EE 上，从状态 α\alpha 到 β\beta 的转移概率被定义为：

k(α,β)≡p(Xn+1=β∣Xn=α)k(\alpha, \beta) \equiv p\left(X_{n+1}=\beta \mid X_{n}=\alpha\right)

对于该转移概率对应的马尔可夫链，如果存在另外一个概率分布 γ\gamma，能够使状态 α\alpha 和 β\beta 之间满足以下等式：

k(α,β)γ(α)=k(β,α)γ(β)∀α,β∈Ek(\alpha, \beta) \gamma(\alpha)=k(\beta, \alpha) \gamma(\beta) \quad \forall \alpha, \beta \in E

则我们称该马尔科夫链是可逆的（或者说，满足细致平衡条件）。而使得上述等式成立的概率分布 γ\gamma 是一个平稳分布。

平稳性证明

对于可逆的马尔可夫链，可以很容易地验证有

∫β∈Ek(β,α)γ(β)dβ=∫β∈Ek(α,β)γ(α)dβ=γ(α)\int_{\beta \in E} k(\beta, \alpha) \gamma(\beta) d \beta=\int_{\beta \in E} k(\alpha, \beta) \gamma(\alpha) d \beta=\gamma(\alpha)

上式符合对平稳分布的定义，因此 γ\gamma 是一个平稳分布。

上面的内容总结起来的核心是： 当一条马尔可夫链对应及其状态分布满足细致平衡条件时，该链是可逆的，且其对应的状态分布为平稳分布。

现在进一步考虑我们关心的后验分布问题，假设我们想要采样的后验概率分布 π\pi 可以定义到一个归一化因子：

π(.)=C×g(.)∝g(.)\pi(.)=C \times g(.) \propto g(.)

则下面的等价性也成立：

k(α,β)π(α)=k(β,α)π(β)∀α,β∈E⟺k(α,β)g(α)=k(β,α)g(β)∀α,β∈E\begin{aligned} k(\alpha, \beta) \pi(\alpha) &=k(\beta, \alpha) \pi(\beta) & & \forall \alpha, \beta \in E \\ \Longleftrightarrow k(\alpha, \beta) g(\alpha) &=k(\beta, \alpha) g(\beta) & & \forall \alpha, \beta \in E \end{aligned}

上式表明，一条转移概率为 k(.，.)k(.，.) 的马尔可夫链，当存在平稳分布 π\pi 时，可以利用最后一个等式，在不计算归一化因子的情况下，获得平稳分布 π\pi 的样本。

5.4 Gibbs 采样转移

假设我们要定义的马尔可夫链中状态维度为 DD ，即第 nn 次迭代对应的状态值为：

Xn=(Xn,1,Xn,2,…,Xn,D)X_{n}=\left(X_{n, 1}, X_{n, 2}, \ldots, X_{n, D}\right)

则 Gibbs 采样算法 基于如下假设：即使联合概率非常复杂，但在给定其他维度时，某个单一维度上的条件分布是便于计算的。基于该思想，可定义如下转移过程，首先生成在迭代 n+1n+1 时要访问的下一个状态提议：

从 XnX_n 的 DD 个维度中随机选择一个维度 dd 。在所有其他维度保持不变的情况下，根据下式的条件概率为该维度采样一个新值：

d∼Uniform⁡({1,2,…,D})d \sim \operatorname{Uniform}(\{1,2, \ldots, D\})

X(n+1),j=Xn,j∀j≠dX(n+1),d∼πd(⋅∣Xn,¬d)\begin{align*} X_{(n+1), j} &= X_{n, j} \quad \forall j \neq d \\ X_{(n+1), d} &\sim \pi_{d}\left(\cdot \mid X_{n, \neg d}\right) \end{align*}

其中：

πd(⋅∣Xn,¬d)=π(Xn,1,…,Xn,(d−1),⋅,Xn,(d+1),…,Xn,D)∫uπ(Xn,1,…,Xn,(d−1),u,Xn,(d+1),…,Xn,D)du=g(Xn,1,…,Xn,(d−1),⋅,Xn,(d+1),…,Xn,D)∫ug(Xn,1,…,Xn,(d−1),u,Xn,(d+1),…,Xn,D)du\pi_{d}\left(\cdot \mid X_{n, \neg d}\right)=\frac{\pi\left(X_{n, 1}, \ldots, X_{n,(d-1)}, \cdot, X_{n,(d+1)}, \ldots, X_{n, D}\right)}{\int_{u} \pi\left(X_{n, 1}, \ldots, X_{n,(d-1)}, u, X_{n,(d+1)}, \ldots, X_{n, D}\right) d u}=\frac{g\left(X_{n, 1}, \ldots, X_{n,(d-1)}, \cdot, X_{n,(d+1)}, \ldots, X_{n, D}\right)}{\int_{u} g\left(X_{n, 1}, \ldots, X_{n,(d-1)}, u, X_{n,(d+1)}, \ldots, X_{n, D}\right) du}

πd(⋅∣Xn,¬d)\pi_{d}\left(\cdot \mid X_{n, \neg d}\right) 代表给定所有其他维度时维度 dd 的条件概率分布。

如果符号

α∼dβ⟺αi=βi∀i≠d\alpha \underset{d}{\sim} \beta \Longleftrightarrow \alpha_{i}=\beta_{i} \quad \forall i \neq d

即该符号代表 dd 维以外的其他维度状态值都不变，是一个逻辑符号。

在形式上，对于任一 $$\beta \neq \alpha$$ 时， 转移概率可以写为：

k(α,β)={1Dg(β)∫γ∼dαg(γ)dγ if β∼dα0 otherwise k(\alpha, \beta)= \{\begin{array}{ll} \frac{1}{D} \frac{g(\beta)}{\int_{\gamma \underset{d} \sim \alpha} g(\gamma) d \gamma} & \text { if } \beta \underset{d}\sim \alpha \\ 0 & \text { otherwise } \end{array}

因此，细致平衡条件成立：

g(α)k(α,β)=1Dg(α)g(β)∫γ∼dαg(γ)dγ=1Dg(β)g(α)∫γ∼ddβg(γ)dγ=g(β)k(β,α)g(\alpha) k(\alpha, \beta)=\frac{1}{D} \frac{g(\alpha) g(\beta)}{\int_{\gamma \underset{d}\sim \alpha} g(\gamma) d \gamma}=\frac{1}{D} \frac{g(\beta) g(\alpha)}{\int_{\gamma \underset{d} \sim{d} \beta} g(\gamma) d \gamma}=g(\beta) k(\beta, \alpha)

也就是说 Gibbs 采样方法给出的状态分布是我们预期的平稳分布。换句话讲，在经过一定次数的 Gibbs 状态转移后，马尔可夫链最终能够收敛到平稳分布，且产生的状态值都是来自该平稳分布的样本。

5.5 Metropolis-Hasting 转移

有时，即便是 Gibbs 方法 ， 其涉及的条件分布也过于复杂，难以得到。此时，可以使用 Metropolis-Hasting 方法。

为此，先定义转移概率 h(.,.)h(.,.) ，这将有助于生成转移提议。在迭代 n+1n+1 时，通过以下过程给出马尔可夫链要访问的下一个状态。

首先从 hh 中抽取一个 “提议状态” xx ，并计算出接受它的概率 rr ：

x∼h(Xn,.) and r=min⁡(1,g(x)h(x,Xn)g(Xn)h(Xn,x))x \sim h\left(X_{n}, .\right) \quad \text { and } \quad r=\min \left(1, \frac{g(\rm{x}) h\left(x, X_{n}\right)}{g\left(X_{n}\right) h\left(X_{n}, x\right)}\right)

然后根据接受概率来选择有效的状态：

Xn+1={x with probability rXn with probability 1−rX_{n+1}=\left\{\begin{array}{ll} x & \text { with probability } r \\ X_{n} & \text { with probability } 1-r \end{array}\right.

在形式上，可以将转移概率写为：

k(α,β)=h(α,β)min⁡(1,g(β)h(β,α)g(α)h(α,β))∀β≠αk(\alpha, \beta)=h(\alpha, \beta) \min \left(1, \frac{g(\beta) h(\beta, \alpha)}{g(\alpha) h(\alpha, \beta)}\right) \quad \forall \beta \neq \alpha

因此，细致平衡条件能够得到保证：

g(α)k(α,β)=g(α)h(α,β)min⁡(1,g(β)h(β,α)g(α)h(α,β))=min⁡(g(α)h(α,β),g(β)h(β,α))=g(β)h(β,α)min⁡(1,g(α)h(α,β)g(β)h(β,α))=g(β)k(β,α)\begin{aligned} g(\alpha) k(\alpha, \beta) &=g(\alpha) h(\alpha, \beta) \min \left(1, \frac{g(\beta) h(\beta, \alpha)}{g(\alpha) h(\alpha, \beta)}\right)=\min (g(\alpha) h(\alpha, \beta), g(\beta) h(\beta, \alpha)) \\ &=g(\beta) h(\beta, \alpha) \min \left(1, \frac{g(\alpha) h(\alpha, \beta)}{g(\beta) h(\beta, \alpha)}\right)=g(\beta) k(\beta, \alpha) \end{aligned}

5.6 采样过程

一旦定义了一条马尔可夫链，我们就可以模拟一个随机的状态序列（随机初始化），并保持其中的一些状态被保留，以便获得遵循目标分布且独立的样本。

首先，为拥有遵循目标分布的样本，只需要考虑距离序列起点足够远的状态，以便确保已经达到马尔可夫链的稳定状态。因此，最初的模拟状态不能作为样本使用，通常将达到稳定所需的这一阶段称为老化时间。但请注意，实践很难知道老化时间具体需要多长时间。

其次，为拥有独立的样本，不能保留老化时间之后的所有后续状态。事实上，马尔可夫链的定义意味着两个连续状态之间存在很强相关性，因此需要将彼此相距足够远的状态作为样本，此时可认为它们几乎独立。实践可通过分析自相关函数来估计被认为几乎独立的两个状态之间所需的滞后（lag）。

因此，为获得遵循目标分布的独立样本，我们需要保留所生成序列中彼此之间以滞后 LL 并且在老化时间 BB 之后的状态。因此，如果马尔可夫链的连续状态序列被表示为：

(Xn)n≥0=X0,X1,X2,…\left(X_{n}\right)_{n \geq 0}=X_{0}, X_{1}, X_{2}, \ldots

则我们只保留样本：

XB,XB+L,XB+2L,XB+3L,…X_{B}, X_{B+L}, X_{B+2 L}, X_{B+3 L}, \ldots

图 5 MCMC 采样需要同时考虑老化时间和滞后

6 变分推断方法

另一种方法是使用变分推断方法，该方法的思路是在某个参数化的分布族中寻找与目标分布最接近的那个分布作为近似分布。为找到这个最佳近似，我们实际是在该分布族的参数空间上执行一个最优化过程。

6.1 近似的方法

变分推断方法 在给定分布族中寻找某一复杂目标分布的最佳近似。更具体地说，其思想是定义一个参数化的分布族，并对参数进行优化，以获得相对于某一误差测度而言最接近目标的结果。

仍然考虑定义到归一化因子 CC 的概率分布 π\pi：

π(.)=C×g(.)∝g(.)\pi(.)=C \times g(.) \propto g(.)

如果将参数化的分布族表示为：

Fν={fω;ω∈ν},ν≡set of all possible parameters\mathcal{F}_{\nu}=\left\{f_{\omega} ; \omega \in \nu\right\}, \nu \equiv \text{set of all possible parameters}

考虑两个分布 pp 和 qq 之间的误差测度 E(p，q)E(p，q)，我们的目的是在 ν\nu 中寻找能够使该误差最小的最优参数，：

ω∗=arg⁡min⁡ω∈νE(fω,π)\omega^{*}=\underset{\omega \in \nu}{\arg \min } E\left(f_{\omega}, \pi\right)

如果能够在不显式归一化 π\pi 的情况下解决此最小化问题，则能够在不用处理棘手的计算问题同时，得到真实后验分布的合理近似 fw⋆f_{w^\star} ，并可进一步将其用于估计更多需要的统计量。实际上，变分推断方法所隐含的优化问题比 “归一化”、“组合” 等直接计算要简单得多。

与前述 MCMC 方法不同，变分方法需要假设了一个参数化的分布族模型，这意味着可能存在着较大的偏差（bias）和较低的方差（variance）。也正因为如此，变分推断方法的精度通常不如 MCMC 方法，但其产生结果的速度则要快得多，更适用于大规模、统计类的问题。

图 6 变分推断近似方法的图示

6.2 分布族的构造

由此可见，变分方法的关键要素是构造一个参数化分布族，该分布族定义了搜索最佳近似分布的参数空间。分布族的选择定义了一个模型，同时控制了方法的偏差和复杂性。如果假设了一个非常简单的模型（简单的分布族），那么我们就有较大的偏差（bias），但优化过程很简单。相反，如果假设一个相当自由的模型（复杂的分布族），则偏差会小得多，但优化会更困难。因此，必须在 “足够复杂以确保最终近似质量的分布族” 和 “足够简单以使优化过程易于处理的分布族” 之间找到适当平衡。当然应谨记，如果分布族中没有任何接近目标分布的分布，那么即使最好的近似也可能给出很差的结果。

平均场变分族是一个将所有未知随机变量都视为独立变量的概率分布族。此分布族中的所有分布都具有乘积性密度，每个乘性的独立分量都受不同因子控制（例如：假设每个分量都为高斯分布，则每个乘性分量均受均值和方差两个因子控制）。因此，平均场变分族中的分布具有如下密度公式：

f(z)=∏j=1mfj(zj)f(z)=\prod_{j=1}^{m} f_{j}\left(z_{j}\right)

其中 zz 为 mm 维随机变量。请注意，式中所有的密度函数 fjf_j 都是参数化的。例如，如果每个密度 fjf_j 都是具有均值和方差参数的高斯函数，则全局密度 ff 由来自 mm 个独立因子的 m×2m \times 2 个参数定义，变分方法将对这组参数进行优化。

图 7 变分推断中分布族的选择决定了优化过程的难度和最终近似的质量

6.3 误差测度的选择 – KL 散度

一旦定义了分布族，下一个主要问题就是：如何在该分布族中找到给定概率分布的最佳近似？显然，最佳近似取决于误差测度的值，而且该最小化问题似乎不应该对归一化因子敏感，因为我们希望更多地比较分布（相对值）而不是概率质量（绝对值）。

因此，定义对归一化因子不敏感的 KL 散度 作为误差测度。在信息论中，如果 pp 和 qq 是两个分布，则 KL 散度 定义如下：

KL(p,q)=Ez∼p[log⁡p(z)]−Ez∼p[log⁡q(z)]K L(p, q)=\mathbb{E}_{z \sim p}[\log p(z)]-\mathbb{E}_{z \sim p}[\log q(z)]

从这个定义中，可以很容易地看出：

KL(fω,Cg)=Ez∼fω[log⁡fω(z)]−Ez∼fω[log⁡(Cg(z))]=Ez∼fω[log⁡fω(z)]−Ez∼fω[log⁡g(z)]−log⁡CK L\left(f_{\omega}, C g\right)=\mathbb{E}_{z \sim f_{\omega}}\left[\log f_{\omega}(z)\right]-\mathbb{E}_{z \sim f_{\omega}}[\log (C g(z))]=\mathbb{E}_{z \sim f_{\omega}}\left[\log f_{\omega}(z)\right]-\mathbb{E}_{z \sim f_{\omega}}[\log g(z)]-\log C

对于我们的最小化问题，这意味着下面的等式成立：

ω∗=arg⁡min⁡ω∈νKL(fω,π)=arg⁡min⁡ω∈νKL(fω,C⋅g)=arg⁡min⁡ω∈νKL(fω,g)\omega^{*}=\underset{\omega \in \nu}{\arg \min } K L\left(f_{\omega}, \pi\right)=\underset{\omega \in \nu}{\arg \min } K L\left(f_{\omega}, C \cdot g\right)=\underset{\omega \in \nu}{\arg \min } K L\left(f_{\omega}, g\right)

因此，当选择 KL 散度 作为误差测度时，优化过程对归一化因子不敏感，可以在参数化分布族中搜索最佳近似值，而不必计算目标分布的归一化因子。

最后，读者应注意到 KL 散度 是交叉熵减去熵，这在信息论中有很好的解释。

注： 这里有些说法容易引起误解，需要修改。其实正是因为 KL 散度 无法避免计算对数似然的问题，才引出了替代性的目标函数 — 变分下界 ELBO 。也就是说，变分下界才是真正与归一化因子无关的误差测度，详情见下节的推导。

6.4 优化流程

一旦定义了参数化分布族和误差测度，就可以（随机或根据明确定义的策略）初始化参数并进行优化。可以使用经典的优化技术，例如 梯度下降法 或 坐标下降法，它们在实践中将会得到局部最优解。

为了更好地理解这个优化过程，让我们举个例子。回到贝叶斯推断问题的具体情况，根据后验公式，后验分布满足：

p(z∣x)∝p(x∣z)p(z)=p(x,z)p(z \mid x) \propto p(x \mid z) p(z)=p(x, z)

如果我们想要使用变分推断得到该后验结果的近似值，必须解决以下优化过程（假设定义了参数化分布族，并将 KL 散度 作为误差测度）：

ω∗=arg⁡min⁡ω∈νKL(fω(z),p(z∣x))=arg⁡min⁡ω∈νKL(fω(z),p(x,z))=arg⁡max⁡ω∈ν(−KL(fω(z),p(x,z)))=arg⁡max⁡ω∈ν(Ez∼fω[log⁡p(z)]+Ez∼fω[log⁡p(x∣z)]−Ez∼fω[log⁡fω(z)])=arg⁡max⁡ω∈ν(Ez∼fω[log⁡p(x∣z)]−KL(fω,p(z)))\begin{aligned} \omega^{*} &=\underset{\omega \in \nu}{\arg \min } K L\left(f_{\omega}(z), p(z \mid x)\right) \\ &=\underset{\omega \in \nu}{\arg \min } K L\left(f_{\omega}(z), p(x, z)\right) \\ &=\underset{\omega \in \nu}{\arg \max }\left(-K L\left(f_{\omega}(z), p(x, z)\right)\right) \\ &=\underset{\omega \in \nu}{\arg \max }\left(\mathbb{E}_{z \sim f_{\omega}}[\log p(z)]+\mathbb{E}_{z \sim f_{\omega}}[\log p(x \mid z)]-\mathbb{E}_{z \sim f_{\omega}}\left[\log f_{\omega}(z)\right]\right) \\ &=\underset{\omega \in \nu}{\arg \max }\left(\mathbb{E}_{z \sim f_{\omega}}[\log p(x \mid z)]-K L\left(f_{\omega}, p(z)\right)\right) \end{aligned}

最后一个等式帮助我们更好地理解近似值是如何被鼓励分配其质量的。其中第一项是期望的对数似然，它倾向于调整参数以使将近似能够将质量放在隐变量 zz 的值上，以更好地解释观测数据 xx。第二项是近似和先验之间的负 KL 散度，它倾向于调整参数以使近似接近先验分布。因此，此目标函数很好地表达了常说的先验/似然平衡。

图 8 变分推断法的优化过程

7 总结

本文的主要内容是：

  • 贝叶斯推断是统计学和机器学习中的一个经典问题，它依赖于贝叶斯定理，主要缺点是在大多数情况下需要进行一些非常繁重的计算。
  • 马尔可夫链蒙特卡罗（MCMC）方法的目的是从复杂概率分布中获取样本，以便用大量样本来模拟真实分布。
  • MCMC 可以用在贝叶斯推断中，以便直接从后验中生成待处理的样本，而不用处理棘手的计算问题。
  • 变分推断(VI)是一种近似分布的方法，它使用参数优化过程在指定分布族中寻找最佳近似。
  • 变分推断的优化过程对目标分布的乘性常数（如：归一化因子）不敏感，因此可以用来近似定义到归一化因子的后验。

如前所述，MCMC 和 VI 方法 具有不同性质，因此适用于不同的应用场景。

  • MCMC 方法的采样过程相当繁重，但偏差小，因此，适用于期望得到准确结果，而不考虑所需时间的场景。
  • 虽然VI 方法中分布族的选择会明显地引入偏差，但它伴随着一个合理的优化过程，非常适合于需要快速计算的超大规模推断问题。

MCMC 和 变分推断 之间的其他比较可在论文《变分推断：统计学回顾》中找到，我们也强烈推荐给只对 变分推断 感兴趣的读者。对于更多关于 MCMC 的阅读，推荐这本《一般性介绍》以及这本《面向机器学习的介绍》。有兴趣了解更多关于应用于 LDA 和 Gibbs 采样 的读者可以参考《主题建模和 Gibbs 采样》，建议结合《关于 LDA 和 Gibbs 采样的课堂讲稿》以获得严谨的推导。