Part I : 从极大似然到生成模型

1.课程概览

这部分我们重点学习概率推断和隐变量建模。首先回顾极大似然估计和最大后验估计，理解如何用参数化型 p(.|θ) 来近似真实但未知的概率函数 p^*(.)，然后介绍EM算法的定义、性质和收敛性分析，再给出一些EM算法的变种，用于解决其局限，比如E步中的难以计算的积分用蒙特卡洛代替、梯度下降、随机梯度下降、小批量梯度下降法等，M步中的参数块优化、降低最大化要求等等

最后将EM算法的思想扩展到深度生成模型，介绍深度隐变量模型的基本框架，推导ELBO详细公式，推导E步和M步的最终公式，包括其中的各种细节，比如说EM可以被理解为在ELBO上进行坐标上升，推导对 \boldsymbol{\theta} 和 \boldsymbol{\phi} 的梯度，有无解析解，重参数化技巧使得梯度可以反向传播，以及最后的完整算法流程中的Cholesky分解等等，最后得到完整的 VAE 流程图

2. MLE/MAP方法回顾

基本框架设定

一个真实但未知的概率函数 p^*(.)，这个真实函数描绘了数据的生成机制，但是他是未知的，因此我们用一个参数化的概率模型 p(.|θ) 来近似这个真实的概率函数 p^*(.)，因此我们就必须学习到这个估计概率模型的未知参数 θ \in \Theta （ \Theta 是参数空间，定义了参数可能的取值范围）。

学习这个估计模型的本质就是找到一个最优的参数 θ，使得我们的估计分布 p(.|θ) 与真实分布 p^*(.) 之间的距离度量达到最小。这里的距离概念会在后续的KLD（Kullback-Leibler散度）部分进一步明确。

判别模型与生成模型的区分

当观测数据可以被分解为两个部分时，问题会呈现出不同的建模视角。具体来说，假设观测可以分为协变量（输入）y^{(1)} \in \mathbb{R}^{p_1 \times 1} 和预测变量（输出）y^{(2)} \in \mathbb{R}^{p_2 \times 1}，其中总维度满足 p = p_1 + p_2。此时数据生成模型可以采用条件模型或联合模型两种方式。

判别模型使用条件模型的形式,即

p(\boldsymbol{y}^{(2)}|\boldsymbol{y}^{(1)}, \boldsymbol{\theta}) = p(\boldsymbol{y}|\boldsymbol{\theta})

这种形式在分类和回归等任务中很常见。判别模型直接对预测目标给定输入的条件分布建模,它关心的是"给定输入,输出应该是什么"。

与此相对,生成模型使用联合模型的形式

p(\boldsymbol{y}^{(1)}, \boldsymbol{y}^{(2)}|\boldsymbol{\theta})

这在贝叶斯估计等场景中使用。生成模型对输入和输出的联合分布建模,它不仅关心输入到输出的映射,还关心数据本身是如何生成的。这种建模方式的优势在于,一旦我们有了联合分布,我们既可以进行预测(通过条件化),也可以生成新的数据样本,还可以处理缺失数据等问题。

后续内容主要聚焦生成模型视角（联合分布）。不过需要注意，从联合分布转换到条件分布在大多数情况下都很直接，所以生成模型的方法可以很容易地应用到判别模型上。

2.1 极大似然估计的定义

MLE的核心思想

寻找参数 \boldsymbol{\theta} 的一个自然方法就是使用极大似然估计。MLE的目标本质上是要找到一种最优的方式来让我们的概率分布拟合观测到的数据。用数学语言来表达，我们要找的是

\hat{\boldsymbol{\theta}}_{\text{MLE}} = \arg\max_{\boldsymbol{\theta} \in \Theta} p(\boldsymbol{y}|\boldsymbol{\theta})\tag{1}

这个公式告诉我们，MLE就是在整个参数空间 \Theta 中搜索那个使得观测数据 \boldsymbol{y} 的概率最大的参数值。由于对数函数是单调递增的，对一个函数取对数不会改变其最大值点的位置，因此在实践中我们通常选择最大化对数似然。所以MLE的等价形式可以写为

\hat{\boldsymbol{\theta}}_{\text{MLE}} = \arg\max_{\boldsymbol{\theta} \in \Theta} \log p(\boldsymbol{y}|\boldsymbol{\theta})\tag{2}

当我们有多个独立的观测样本时，情况会变得更加具体。假设我们有 n 个独立观测 \{\boldsymbol{y}_i\}_{i=1}^n，那么根据独立性假设，联合概率可以分解为各个观测概率的乘积，对数似然就变成了求和的形式：

\hat{\boldsymbol{\theta}}_{\text{MLE}} = \arg\max_{\boldsymbol{\theta} \in \Theta} \log p(\boldsymbol{y}|\boldsymbol{\theta}) = \arg\max_{\boldsymbol{\theta} \in \Theta} \sum_{i=1}^n \log p(\boldsymbol{y}_i|\boldsymbol{\theta})\tag{3}

从目标函数的本质来看，MLE是一个让已经观测到的观测数据最有可能发生的估计量。换句话说，我们已经观测到了某些数据，那么最好的模型参数应该就是那个使得这些已经发生的事情最有可能发生的参数。

可识别性的前提条件

但在实际应用MLE之前，我们必须确保数据生成过程包含了关于模型参数的充分信息，并且模型具有全局可识别性。具体来说，这个可识别性条件要求：对于参数空间中的任意两个不同参数 \boldsymbol{\theta}_1 \in \Theta 和 \boldsymbol{\theta}_2 \in \Theta，只要 \boldsymbol{\theta}_1 \neq \boldsymbol{\theta}_2，那么对于任意的观测 \boldsymbol{y} \in \mathcal{Y}，必须有

p(\boldsymbol{y}|\boldsymbol{\theta}_1) \neq p(\boldsymbol{y}|\boldsymbol{\theta}_2)

这个条件保证了不同的参数会产生不同的数据分布，从而使得我们可以通过观测数据来唯一地识别出真实的参数，换句话说，我们最后求出的这组参数就对应着我们估计的那个分布，关系是绑定的。

从MLE到MAP

在某些情况下，我们可能拥有关于参数的先验知识，或者我们想要避免对数似然函数出现奇异性（混合模型中混合分布的数量远大于可用数据）的问题，就好比 (H^TH)^{-1} 中不适定问题。在这些场景下，我们可以引入先验信息，这就引出了最大后验估计。MAP方法直接采用贝叶斯的观点，将参数也看作随机变量，其定义为

\hat{\boldsymbol{\theta}}_{\text{MAP}} = \arg\max_{\boldsymbol{\theta} \in \Omega} \log p(\boldsymbol{\theta}|\boldsymbol{y}) = \arg\max_{\boldsymbol{\theta} \in \Omega} \left(\log p(\boldsymbol{y}|\boldsymbol{\theta}) + \log p(\boldsymbol{\theta})\right) \tag{4}

这个公式来源于贝叶斯定理。我们知道后验分布 p(\boldsymbol{\theta}|\boldsymbol{y}) 正比于似然 p(\boldsymbol{y}|\boldsymbol{\theta}) 与先验 p(\boldsymbol{\theta}) 的乘积，因此取对数后就变成了两项的和。MAP相比于MLE多了一个 \log p(\boldsymbol{\theta}) 项，这个先验项可以起到正则化的作用，防止参数估计过拟合，同时也能在数据不足时提供我们对参数的先验认知。

2.2 MLE与KL散度的联系

KL散度的定义与性质

MLE的优化目标实际上可以从信息论的角度来理解，它本质上是在求解Kullback-Leibler散度的最小化问题。

KL散度是衡量两个概率分布之间差异的一个度量，对于两个分布 p(\boldsymbol{y}) 和 q(\boldsymbol{y})，它们之间的KL散度定义为

\text{KLD}(p(\boldsymbol{y})||q(\boldsymbol{y})) = \mathbb{E}_{\boldsymbol{y} \sim p(\boldsymbol{y})}\left[\log \frac{p(\boldsymbol{y})}{q(\boldsymbol{y})}\right] = \int p(\boldsymbol{y}) \log \frac{p(\boldsymbol{y})}{q(\boldsymbol{y})} d\boldsymbol{y}\tag{5}

当我们用分布 p 来采样数据时，用 p 来编码这些数据相比于用 q 来编码，平均能节省多少信息量。KL散度有几个重要的性质我们需要明确。

首先，KL散度是非对称的，也就是说

\text{KLD}(p(\boldsymbol{y})||q(\boldsymbol{y})) \neq \text{KLD}(q(\boldsymbol{y})||p(\boldsymbol{y}))

从 p 到 q 的距离和从 q 到 p 的距离是不同的。
其次，KL散度总是非负的，即 \text{KLD}(p(\boldsymbol{y})||q(\boldsymbol{y})) \geq 0。
最后，KL散度等于零当且仅当两个分布完全相同，也就是

\text{KLD}(p(\boldsymbol{y})||q(\boldsymbol{y})) = 0 \Leftrightarrow p(\boldsymbol{y}) = q(\boldsymbol{y})

MLE如何最小化KL散度

现在我们来看MLE是如何与KL散度联系起来的。假设我们有一个训练集 \{\boldsymbol{y}_i\}_{i=1}^n，其中每个样本 \boldsymbol{y}_i 都是从真实但未知的分布 p^*(\boldsymbol{y}) 中独立采样得到的。我们的目标是找到一个参数 \boldsymbol{\theta}，使得模型分布 p(\boldsymbol{y}|\boldsymbol{\theta}) 尽可能接近真实分布 p^*(\boldsymbol{y})。如果我们用KL散度来衡量这个接近程度，那么优化目标就是

\hat{\boldsymbol{\theta}}_{\text{KLD}} = \arg\min_{\boldsymbol{\theta} \in \Theta} \text{KLD}(p^*(\boldsymbol{y})||p(\boldsymbol{y}|\boldsymbol{\theta})) = \arg\min_{\boldsymbol{\theta} \in \Theta} \mathbb{E}_{\boldsymbol{y} \sim p^*(\boldsymbol{y})}\left[\log \frac{p^*(\boldsymbol{y})}{p(\boldsymbol{y}|\boldsymbol{\theta})}\right]

我们把期望展开，得到

\mathbb{E}_{\boldsymbol{y} \sim p^*(\boldsymbol{y})}[\log p^*(\boldsymbol{y})] - \mathbb{E}_{\boldsymbol{y} \sim p^*(\boldsymbol{y})}[\log p(\boldsymbol{y}|\boldsymbol{\theta})]

第一项 \mathbb{E}_{\boldsymbol{y} \sim p^*(\boldsymbol{y})}[\log p^*(\boldsymbol{y})] 是真实分布的熵，它与我们要优化的参数 \boldsymbol{\theta} 无关，因此在优化过程中可以看作常数。所以最小化KL散度等价于最大化第二项，也就是

\hat{\boldsymbol{\theta}}_{\text{KLD}} = \arg\max_{\boldsymbol{\theta} \in \Theta} \mathbb{E}_{\boldsymbol{y} \sim p^*(\boldsymbol{y})}[\log p(\boldsymbol{y}|\boldsymbol{\theta})]

但问题是我们并不知道真实分布 p^*(\boldsymbol{y})，所以无法直接计算这个期望。这时我们可以用训练样本的经验平均来近似这个期望，即

\mathbb{E}_{\boldsymbol{y} \sim p^*(\boldsymbol{y})}[\log p(\boldsymbol{y}|\boldsymbol{\theta})] \approx \frac{1}{n} \sum_{i=1}^n \log p(\boldsymbol{y}_i|\boldsymbol{\theta})

这样我们就得到了

\hat{\boldsymbol{\theta}}_{\text{KLD}} = \arg\max_{\boldsymbol{\theta} \in \Theta} \frac{1}{n} \sum_{i=1}^n \log p(\boldsymbol{y}_i|\boldsymbol{\theta}) = \hat{\boldsymbol{\theta}}_{\text{MLE}}

这个推导揭示了MLE的深层含义：MLE实际上是在最小化模型分布与真实数据分布之间的KL散度的经验估计。从这个角度看，MLE不仅仅是一个让数据最有可能出现的估计量，更是一个让模型分布尽可能接近真实数据生成分布的方法。

例2.1：多元高斯观测的MLE

问题设定

现在我们通过一个具体例子来看看MLE是如何在实践中工作的。假设我们有一组独立同分布的观测 \{\boldsymbol{y}_i\}_{i=1}^n，每个观测都服从多元高斯分布

\boldsymbol{y}_i \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})\tag{6}

这里 \boldsymbol{y}_i \in \mathbb{R}^p 是 p 维观测向量，\boldsymbol{\mu} \in \mathbb{R}^p 是均值向量，\boldsymbol{\Sigma} \in \mathbb{R}^{p \times p} 是对称正定的协方差矩阵。多元高斯分布的概率密度函数写为

p(\boldsymbol{y}_i|\boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{1}{(2\pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}} \exp\left(-\frac{1}{2}(\boldsymbol{y}_i - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1}(\boldsymbol{y}_i - \boldsymbol{\mu})\right)

构建对数似然函数

由于观测是独立的，对数似然函数就是各个观测的对数概率之和。我们写出完整的对数似然

\ell(\boldsymbol{\mu}, \boldsymbol{\Sigma}) = \sum_{i=1}^n \log p(\boldsymbol{y}_i; \boldsymbol{\mu}, \boldsymbol{\Sigma}) \propto -\frac{n}{2}\log|\boldsymbol{\Sigma}| - \frac{1}{2}\sum_{i=1}^n (\boldsymbol{y}_i - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1}(\boldsymbol{y}_i - \boldsymbol{\mu})

这里我们用了比例符号 \propto，因为那些与参数无关的常数项（如 -\frac{np}{2}\log(2\pi)）在求最大值时不影响结果，可以忽略。这个对数似然函数第一项 -\frac{n}{2}\log|\boldsymbol{\Sigma}| 来自归一化常数，它鼓励协方差矩阵的行列式变小，第二项是所有观测的马氏距离平方和，它希望数据点尽可能接近均值。

求解均值的MLE

我们首先对均值 \boldsymbol{\mu} 求偏导。由于第一项不含 \boldsymbol{\mu}，我们只需要对第二项求导。利用矩阵微分的链式法则，对于二次型 (\boldsymbol{y}_i - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1}(\boldsymbol{y}_i - \boldsymbol{\mu}) 对 \boldsymbol{\mu} 的导数是 -2\boldsymbol{\Sigma}^{-1}(\boldsymbol{y}_i - \boldsymbol{\mu})，因此

\frac{\partial \ell(\boldsymbol{\mu}, \boldsymbol{\Sigma})}{\partial \boldsymbol{\mu}} = \boldsymbol{\Sigma}^{-1} \sum_{i=1}^n (\boldsymbol{y}_i - \boldsymbol{\mu})

令这个梯度等于零，由于 \boldsymbol{\Sigma}^{-1} 是可逆的，我们可以消去它，得到 \sum_{i=1}^n (\boldsymbol{y}_i - \boldsymbol{\mu}) = 0，展开后移项就得到均值的MLE

\hat{\boldsymbol{\mu}} = \frac{1}{n}\sum_{i=1}^n \boldsymbol{y}_i\tag{7}

多元高斯分布均值的最大似然估计就是样本均值，这是我们最熟悉的统计量。

求解协方差矩阵的MLE

接下来我们要找协方差矩阵 \boldsymbol{\Sigma} 的MLE。我们需要对 \boldsymbol{\Sigma} 求导。对数似然对 \boldsymbol{\Sigma} 的导数涉及两个技巧：首先是 \log|\boldsymbol{\Sigma}| 对 \boldsymbol{\Sigma} 的导数是 \boldsymbol{\Sigma}^{-1}；其次是二次型对协方差矩阵的导数。经过推导我们得到

\frac{\partial \ell(\hat{\boldsymbol{\mu}}, \boldsymbol{\Sigma})}{\partial \boldsymbol{\Sigma}} = -\frac{n}{2}\boldsymbol{\Sigma}^{-1} + \frac{1}{2}\sum_{i=1}^n \boldsymbol{\Sigma}^{-1}(\boldsymbol{y}_i - \hat{\boldsymbol{\mu}})(\boldsymbol{y}_i - \hat{\boldsymbol{\mu}})^\top \boldsymbol{\Sigma}^{-1}

这个导数的第一项来自 \log|\boldsymbol{\Sigma}|，第二项来自二次型。令导数等于零，两边同时左乘 \boldsymbol{\Sigma} 右乘 \boldsymbol{\Sigma}，我们得到

\hat{\boldsymbol{\Sigma}} = \frac{1}{n}\sum_{i=1}^n (\boldsymbol{y}_i - \hat{\boldsymbol{\mu}})(\boldsymbol{y}_i - \hat{\boldsymbol{\mu}})^\top\tag{8}

这就是样本协方差矩阵，它刻画了数据在各个方向上的变异性。

MLE闭式解的局限性

这个例子我们能够通过解析的方式求出MLE的闭式表达式。但MLE并不总能像这样找到闭式解。在很多实际问题中，对数似然函数可能是高度非凸的，或者包含复杂的隐变量结构，这时就需要借助数值优化方法或者更高级的算法框架，比如我们接下来要学习的EM算法。

例2.2：学习固定模型的最优混合

混合模型的问题设定

现在我们来看一个更复杂但也更实际的例子，这个例子会揭示为什么我们需要EM算法。假设我们面临这样一个学习问题：需要找到若干已知模型的最优组合。具体来说，假设我们观测到数据 \boldsymbol{y} = [\boldsymbol{y}_1, \ldots, \boldsymbol{y}_n]^\top，这些数据可能是由 k 个不同的模型 (p_1, \ldots, p_k) 生成的。更进一步，我们假设每个观测都是独立同分布地从这 k 个模型的某种组合中产生的，这种组合可以用混合分布来表达

p(\boldsymbol{y}_i|\boldsymbol{\theta}) = \sum_{j=1}^k \theta_j p_j(\boldsymbol{y}_i)\tag{9}

这里的关键参数是混合系数 \boldsymbol{\theta} = [\theta_1, \ldots, \theta_k]^\top，它们必须满足两个约束条件：

首先所有混合系数之和必须等于1，即 \sum_{j=1}^k \theta_j = 1，这保证了混合分布是一个合法的概率分布
其次每个混合系数必须非负且不超过1，即 \theta_j \in [0,1] 对所有 j 成立。

我们的目标就是学习出这个最可能的混合系数组合，使得它最好地刻画了观测到的数据 \boldsymbol{y}。

混合模型MLE的困难

在这种情况下，极大似然估计的目标函数写为

\hat{\boldsymbol{\theta}} = \arg\max_{\boldsymbol{\theta}} \sum_{i=1}^n \log \sum_{j=1}^k \theta_j p_j(\boldsymbol{y}_i)\tag{10}

从公式(10)我们立刻能看出一个严重的问题：\boldsymbol{\theta} 的MLE不存在闭式解。这个困难主要源于对数内部的求和 \sum_{j=1}^k \theta_j p_j(\boldsymbol{y}_i) 将所有的混合系数 \theta_j 耦合在了一起，使得我们无法像之前那样分别对各个参数求导并令其为零来得到解析解。本质上来说，这里面藏着一个隐变量问题，就是我们不清楚数据是来自哪个分布的，这个分布就是隐变量

高斯混合模型的现实挑战

如果 p_j(\boldsymbol{y}_i) 的每个成分都是高斯分布，那么就是高斯混合模型（GMM）的情况，这时每个分量 p_j 具有未知的均值 \boldsymbol{\mu}_j 和协方差矩阵 \boldsymbol{\Sigma}_j，即

p_j(\boldsymbol{y}_i|\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j) = \mathcal{N}(\boldsymbol{y}_i; \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)

此时极大似然变成

\hat{\boldsymbol{\theta}}, \{\hat{\boldsymbol{\mu}}_j, \hat{\boldsymbol{\Sigma}}_j\}_j = \arg\max_{\boldsymbol{\theta}, \{\boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j\}_j} \sum_{i=1}^n \log \sum_{j=1}^k \theta_j \mathcal{N}(\boldsymbol{y}_i; \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)\tag{11}

我们要找到最优的参数：

\theta_j：每个高斯分布的权重；
\mu_j：第 j 个高斯分布的均值；
\Sigma_j：第 j 个高斯分布的协方差矩阵；

这个优化问题只有在非常特定和受限的场景下才可能有闭式解，比如当 k=1 时，对数内的求和就消失了，问题就退化到了例2.1的情况。

隐变量的启发

对于第 j 个高斯分量，假设我们已经知道哪些数据点属于它，也就是如果我们知道每个观测是由哪个分布生成的，问题就会大大简化。

第 j 个高斯分布的概率权重(混合系数)的估计就是

\hat{\theta}_j = \frac{n_j}{n}

n_j 是由第 j 个模型生成的观测数量

第 j 个高斯分量的均值估计变成

\hat{\mu}_j = \frac{\sum_{i=1}^{n} \boldsymbol{y}_i \cdot \mathbb{I}_{j \to i}}{\sum_{i=1}^{n} \mathbb{I}_{j \to i}} = \frac{1}{n_j} \sum_{i=1}^{n} \boldsymbol{y}_i \mathbb{I}_{j \to i}

第 j 个高斯分量的协方差矩阵估计为

\hat{\boldsymbol{\Sigma}}_j = \frac{1}{n_j} \sum_{i=1}^n (\boldsymbol{y}_i - \hat{\boldsymbol{\mu}}_j)(\boldsymbol{y}_i - \hat{\boldsymbol{\mu}}_j)^\top \mathbb{1}_{j \to i}

这里的指示函数 \mathbb{1}_{j \to i} = 1 当且仅当第 j 个模型是第 i 个观测的生成过程，否则 \mathbb{1}_{j \to i} = 0。总的观测数 n_j = \sum_i \mathbb{1}_{j \to i} 就是属于第 j 个分量的样本数。有了这些指示变量，我们实际上就是在对每个分量分别做MLE，每个分量只使用属于它的那些样本，这就回到了例2.1的简单情况。

EM算法的动机

但遗憾的是，我们观测到的只是混合后的数据，而不知道每个数据点来自哪个分量。这就是EM算法要解决的问题：EM算法试图猜测每个数据点来自哪个分量。EM算法特别适用于这样的场景：模型中存在额外的但不可观测的信息，也就是隐变量或潜在变量，如果这些变量能够被观测到，学习任务就会被大大简化。在混合模型中，这个隐变量就是指示每个观测属于哪个分量的标签。EM算法通过迭代地估计这些隐变量和模型参数，最终收敛到一个局部最优解。

2.3 通过随机梯度下降求解MLE

梯度下降法的基本思路

当MLE不存在闭式解时，我们需要求助于数值优化方法。回顾一下在独立性假设下，MLE的目标是

对于不可处理的MLE问题，也就是无法得到闭式解的情况，我们可以利用链式法则和自动微分工具来计算目标函数的梯度，即计算 l(\boldsymbol{\theta}|\boldsymbol{y}) = \sum_{i=1}^n \log p(\boldsymbol{y}_i|\boldsymbol{\theta}) 对 \boldsymbol{\theta} 的梯度。如果我们对每个样本都计算梯度，然后求和，这就得到了所谓的批量梯度下降法。但是当样本数 n 很大时，每次迭代都要遍历所有样本来计算梯度，这在计算上会非常昂贵。

随机梯度下降的优势

为了克服批量梯度下降的计算负担，我们可以使用随机梯度下降法。SGD的核心思想是用一个随机子集（小批量）来近似真实的梯度。具体的更新规则写为

\boldsymbol{\theta}^{(m+1)} \leftarrow \boldsymbol{\theta}^{(m)} + \eta_m \nabla_{\boldsymbol{\theta}} \tilde{l}(\hat{\boldsymbol{\theta}}|\varepsilon)\tag{13}

这里 \eta_m 是第 m 步的学习率，它控制着每次更新的步长大小。关键的部分是 \tilde{l}(\boldsymbol{\theta}|\varepsilon)，这是真实对数似然 l(\boldsymbol{\theta}|\boldsymbol{y}) 的一个无偏估计。所谓无偏估计意味着

E_{\varepsilon \sim p(\varepsilon)}\left\{\tilde{l}(\boldsymbol{\theta}|\varepsilon)\right\} = l(\boldsymbol{\theta}|\boldsymbol{y})\tag{14}

也就是说，虽然每次迭代我们使用的是一个随机的、可能带有噪声的梯度估计，但这个估计在期望意义下是准确的。

证明: 随机梯度下降的估计是真实似然的无偏估计

完整批量梯度：

\nabla l(\boldsymbol{\theta}|\boldsymbol{y}) = \sum_{i=1}^{n} \nabla \log p(\boldsymbol{y}_i|\boldsymbol{\theta})

平均梯度（真实梯度）：

\frac{1}{n}\sum_{i=1}^{n} \nabla \log p(\boldsymbol{y}_i|\boldsymbol{\theta})

随机梯度（随机选一个样本 \varepsilon）

\tilde{l}(\boldsymbol{\theta}|\varepsilon) = \nabla \log p(\boldsymbol{y}_\varepsilon|\boldsymbol{\theta})

假设我们等概率随机选择样本，即每个样本被选中的概率都是 \frac{1}{n}。

E_{\varepsilon \sim p(\varepsilon)}\left[\tilde{l}(\boldsymbol{\theta}|\varepsilon)\right] = E_{\varepsilon}\left[\nabla \log p(\boldsymbol{y}_\varepsilon|\boldsymbol{\theta})\right]

展开期望（对所有可能的样本求加权平均）：

= \sum_{i=1}^{n} P(\varepsilon = i) \cdot \nabla \log p(\boldsymbol{y}_i|\boldsymbol{\theta})

由于等概率选择，P(\varepsilon = i) = \frac{1}{n}：

= \sum_{i=1}^{n} \frac{1}{n} \nabla \log p(\boldsymbol{y}_i|\boldsymbol{\theta}) = \frac{1}{n}\sum_{i=1}^{n} \nabla \log p(\boldsymbol{y}_i|\boldsymbol{\theta})

这正是真实的平均梯度

小批量梯度估计的构造

随机变量 \varepsilon 可以被设计为一个随机均匀地选择小批量数据的机制。具体来说，我们在每次迭代时随机选择一个大小为 n_{\text{mini}} 的小批量 \mathcal{M}_{\text{mini}}，这个小批量包含了 n 个样本中的 n_{\text{mini}} 个。利用这个小批量，我们可以构造出梯度的无偏估计

\frac{1}{n_{\text{mini}}} \sum_{i, \boldsymbol{y}_i \in \mathcal{M}_{\text{mini}}} \nabla_{\boldsymbol{\theta}} l(\boldsymbol{\theta}|\boldsymbol{y}_i)

这个小批量梯度是完整批量梯度的无偏估计

\frac{1}{n} \sum_{i=1}^n \nabla_{\boldsymbol{\theta}} l(\boldsymbol{\theta}|\boldsymbol{y}_i)

这是因为每个样本被选入小批量的概率是均等的，所以对小批量梯度求期望就会得到完整梯度。小批量的大小 n_{\text{mini}} 是一个需要权衡的超参数：太小会导致梯度估计方差很大，收敛不稳定；太大则会增加每次迭代的计算成本，失去了SGD的计算优势。

证明: 小批量梯度是完整真实批量 \frac{1}{n} \sum_{i=1}^{n} \nabla l(\boldsymbol{\theta}|\boldsymbol{y}_i) 的无偏估计

小批量梯度的构造：随机选择 n_{\text{mini}} 个样本组成小批量 \mathcal{M}_{\text{mini}}

小批量梯度估计：

\tilde{l}_{\text{mini}}(\boldsymbol{\theta}) = \frac{1}{n_{\text{mini}}} \sum_{i, \boldsymbol{y}_i \in \mathcal{M}_{\text{mini}}} \nabla l(\boldsymbol{\theta}|\boldsymbol{y}_i)

假设小批量是随机均匀抽取的，那么对于任意样本 \boldsymbol{y}_i：

P(\boldsymbol{y}_i \in \mathcal{M}_{\text{mini}}) = \frac{n_{\text{mini}}}{n}

对小批量梯度求期望：

E_{\mathcal{M}_{\text{mini}}}\left[\tilde{l}_{\text{mini}}(\boldsymbol{\theta})\right] = E\left[\frac{1}{n_{\text{mini}}} \sum_{i, \boldsymbol{y}_i \in \mathcal{M}_{\text{mini}}} \nabla l(\boldsymbol{\theta}|\boldsymbol{y}_i)\right]

利用期望的线性性：

= \frac{1}{n_{\text{mini}}} \sum_{i=1}^{n} E[\mathbb{I}(\boldsymbol{y}_i \in \mathcal{M}_{\text{mini}})] \cdot \nabla l(\boldsymbol{\theta}|\boldsymbol{y}_i)

其中指示函数的期望就是概率：

E[\mathbb{I}(\boldsymbol{y}_i \in \mathcal{M}_{\text{mini}})] = P(\boldsymbol{y}_i \in \mathcal{M}_{\text{mini}}) = \frac{n_{\text{mini}}}{n}

代入得：

= \frac{1}{n_{\text{mini}}} \sum_{i=1}^{n} \frac{n_{\text{mini}}}{n} \nabla l(\boldsymbol{\theta}|\boldsymbol{y}_i) = \frac{1}{n} \sum_{i=1}^{n} \nabla l(\boldsymbol{\theta}|\boldsymbol{y}_i)

证毕

SGD算法流程

随机梯度下降的完整算法流程可以总结如下。在初始化阶段，我们需要选择一个初始参数估计 \boldsymbol{\theta}^{(0)} 以及一系列学习率 \{\eta_m\}。学习率的选择对算法的收敛性能有重要影响，通常我们会使用随迭代次数递减的学习率序列。

然后我们重复以下步骤直到算法收敛：首先在第 m 次迭代时，随机选择一个大小为 n_{\text{mini}} 的小批量 \mathcal{M}_{\text{mini}}^{(m)}；接着使用这个小批量来更新参数

\boldsymbol{\theta}^{(m+1)} \leftarrow \boldsymbol{\theta}^{(m)} + \frac{\eta_m}{n_{\text{mini}}} \sum_{i, \boldsymbol{y}_i \in \mathcal{M}_{\text{mini}}^{(m)}} \nabla_{\boldsymbol{\theta}} l(\boldsymbol{\theta}|\boldsymbol{y}_i)

这个更新规则将当前参数沿着小批量梯度的方向移动一小步。通过不断重复这个过程，参数会逐渐向对数似然的局部最大值靠近。虽然SGD引入了随机性，使得每一步的更新都可能偏离真实的最优方向,但正是这种随机性帮助算法以较小的计算代价在高维参数空间中探索，并最终找到好的解。在实践中，SGD及其各种变体（如带动量的SGD、Adam等）已经成为训练大规模机器学习模型的标准工具。

3 EM算法

3.1 EM算法的定义与性质

3.1.1 EM算法的基本定义

完全数据的概念

EM算法引入了一个关键的概念，那就是完全数据 \boldsymbol{x} ，这个完全数据的参数化概率密度记为 p(\boldsymbol{x}|\boldsymbol{\theta})。而我们实际上并不能观测到这个完全数据 \boldsymbol{x}，我们观测到的只是 \boldsymbol{y}

EM算法通过引入一个我们无法直接观测但能够大大简化问题的隐变量结构，我们可以迭代地逼近原本难以处理的优化问题。

EM算法的两步迭代结构

EM算法的核心由两个交替进行的步骤组成，分别是期望步骤和最大化步骤，这也是算法名称的由来。在开始迭代之前，我们需要进行初始化：当 m=0 时，我们需要给出参数的初始估计 \boldsymbol{\theta}^{(0)}。这个初始值可以通过先验知识来设定，或者使用某个次优但已有的算法来获得。初始化对算法的收敛速度和最终结果都有影响。

E步骤的深层含义

期望步骤是EM算法的第一步。在这一步中，我们基于上一次迭代得到的参数估计 \boldsymbol{\theta}^{(m)}，构造出完全数据在观测数据和当前参数下的条件概率密度 p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})。有了这个后验分布之后，我们计算一个关键的量 Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)})，它的定义是

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \int \log p(\boldsymbol{x}|\boldsymbol{\theta}) p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) d\boldsymbol{x}

这个积分也可以写成期望的形式

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}|\boldsymbol{\theta})\}\tag{15}

这个 Q 函数是完全数据对数似然 \log p(\boldsymbol{x}|\boldsymbol{\theta}) 关于后验分布 p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 的期望。换句话说，由于我们不知道完全数据 \boldsymbol{x} 的真实值，我们就对所有可能的完全数据取一个加权平均，权重就是在当前参数估计下每个完全数据的后验概率。

M步骤的优化目标

在最大化步骤中，我们要找到使得 Q 函数最大化的新参数，即第 m+1 次的参数估计定义为

\boldsymbol{\theta}^{(m+1)} = \arg\max_{\boldsymbol{\theta} \in \Theta} Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)})\tag{16}

这个最大化问题通常比直接最大化原始的对数似然 \log p(\boldsymbol{y}|\boldsymbol{\theta}) 要简单得多。这是因为我们精心设计的完全数据结构使得 \log p(\boldsymbol{x}|\boldsymbol{\theta}) 具有更好的数学形式。M步的目标不是直接最大化观测数据的似然，而是最大化完全数据对数似然的条件期望，这个替代目标函数往往有更好的优化性质。

迭代收敛的过程

完成一次E步和M步后，我们增加迭代计数 m，然后重复上述的E步和M步过程，直到满足某个停止准则。这个停止准则可以是参数的变化量小于某个阈值，或者对数似然的增量足够小，或者达到预设的最大迭代次数。

EM算法的收敛性质

EM估计算法只会单调不降，也就是说，每次迭代后的似然值绝不会比之前差。这个性质我们后面会详细讨论其数学原理。

EM算法通常会找到似然函数 p(\boldsymbol{y}|\boldsymbol{\theta}) 的一个峰值，但如果似然函数存在多个局部最大值，EM不一定能找到全局最大值，算法最终收敛到哪个局部最优解很大程度上取决于初始值的选择。因此在实践中，一个常见的策略是从多个不同的随机初始猜测出发运行EM算法，然后在所有收敛结果中选择具有最大似然值的那个作为 \boldsymbol{\theta} 的最终估计。这种多次随机初始化的策略可以提高找到全局最优解或接近全局最优解的概率。

3.1.2 完全数据的性质要求

完全数据的可优化性

在选择完全数据时，指定的完全数据 X 必须使得对于给定的 \boldsymbol{x}，最大化 p(\boldsymbol{x}|\boldsymbol{\theta}) 相对容易。这个要求直接关系到M步能否有效执行。如果我们选择的完全数据结构使得 p(\boldsymbol{x}|\boldsymbol{\theta}) 的优化仍然非常困难，那么引入完全数据就失去了意义。理想情况下，完全数据应该被设计成使得 Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) 的最大化能够得到闭式解，或者至少比原始问题容易很多。

马尔可夫关系的必要性

除了可优化性之外，完全数据 X 还必须满足马尔可夫关系，图形化地表示为

\boldsymbol{\theta} \to X \to Y\tag{17}

表示条件独立性结构：给定完全数据 X，观测数据 Y 与参数 \boldsymbol{\theta} 条件独立。用概率语言来表达，这意味着

p(\boldsymbol{y}|\boldsymbol{x}, \boldsymbol{\theta}) = p(\boldsymbol{y}|\boldsymbol{x})\tag{18}

这个等式告诉我们，一旦我们知道了完全数据 \boldsymbol{x}，观测数据 \boldsymbol{y} 的分布就不再依赖于参数 \boldsymbol{\theta}。换言之，参数 \boldsymbol{\theta} 对观测数据的影响完全是通过完全数据 \boldsymbol{x} 传递的。

确定性函数的特例

为了更好地理解这个马尔可夫关系，我们可以看一个简单但常见的例子。当完全数据与观测数据通过某个确定性函数相关联时，即存在关系 Y = F(X)，其中 F(\cdot) 是一个确定性函数，那么上述公式所表达的马尔可夫关系自然成立。这是因为如果 Y 完全由 X 通过确定性映射决定，那么知道了 X 之后，Y 就完全确定了，自然不需要再依赖参数 \boldsymbol{\theta}。许多实际的EM应用都属于这种情况，比如在混合模型中，观测数据就是完全数据（包含观测值和隐藏的分量标签）通过去掉标签这个确定性操作得到的。

3.1.3 独立同分布样本情况下的EM形式

i.i.d.假设下的问题简化

在许多实际应用中，我们观测到的数据是独立同分布的。在这种情况下，我们自然希望完全数据也具有相应的独立同分布结构，使得第 i 个观测 \boldsymbol{y}_i 仅仅是完全数据中对应部分 \boldsymbol{x}_i 的函数。这种设定会带来一个重要的简化，使得EM算法的计算可以分解为对各个样本的独立处理。

引理3.1：i.i.d.情况下Q函数的分解

为了形式化这个简化，我们给出下面的引理。假设完全数据的概率密度可以分解为

p(\boldsymbol{x}|\boldsymbol{\theta}) = \prod_{i=1}^n p(\boldsymbol{x}_i|\boldsymbol{\theta})

这里完全数据 \boldsymbol{x} = [\boldsymbol{x}_1^\top, \ldots, \boldsymbol{x}_n^\top]^\top 由 n 个独立的部分组成，每一部分 \boldsymbol{x}_i 都服从相同的参数化分布。同时，马尔可夫关系对每个样本都成立，即

\boldsymbol{\theta} \to X_i \to Y_i, \quad \forall i = 1, \ldots, n

这个条件用概率语言表达就是

p(\boldsymbol{y}_i|\boldsymbol{x}, \boldsymbol{y}_1, \ldots, \boldsymbol{y}_{i-1}, \boldsymbol{y}_{i+1}, \ldots, \boldsymbol{y}_n, \boldsymbol{\theta}) = p(\boldsymbol{y}_i|\boldsymbol{x}_i, \boldsymbol{\theta})\tag{19}

这个等式表明，给定完全数据的第 i 个分量 \boldsymbol{x}_i 后，第 i 个观测 \boldsymbol{y}_i 与其他所有观测、其他所有完全数据分量以及参数都条件独立。

在这些假设下，引理告诉我们Q函数可以分解为各个样本贡献的和

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \sum_{i=1}^n Q_i(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)})\tag{20}

其中每个样本的贡献定义为

Q_i(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = E_{\boldsymbol{x}_i \sim p(\boldsymbol{x}_i|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}_i|\boldsymbol{\theta})\}, \quad \forall i = 1, \ldots, n

这个结果的重要性在于：原本我们需要对整个完全数据向量 \boldsymbol{x} 计算后验分布 p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 并对其求期望，现在我们只需要对每个 \boldsymbol{x}_i 分别计算其在给定对应观测 \boldsymbol{y}_i 下的后验分布 p(\boldsymbol{x}_i|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)})，然后求和即可。这大大降低了计算复杂度，尤其是当样本数很大时。

证明过程：建立样本的相互独立性

现在我们来证明这个引理。证明分为两个部分。首先我们要证明完全数据和观测数据的联合分布具有样本间的相互独立性。记所有观测为 \boldsymbol{y} = [\boldsymbol{y}_1^\top, \ldots, \boldsymbol{y}_n^\top]^\top，我们从联合分布的链式法则开始

p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta}) = p(\boldsymbol{y}_1|\boldsymbol{y}_2, \ldots, \boldsymbol{y}_n, \boldsymbol{x}, \boldsymbol{\theta}) \ldots p(\boldsymbol{y}_n|\boldsymbol{x}, \boldsymbol{\theta}) p(\boldsymbol{x}|\boldsymbol{\theta})

这是概率的链式法则，我们将联合概率分解为一系列条件概率的乘积。接下来我们应用条件独立性，每个 \boldsymbol{y}_i 给定 \boldsymbol{x}_i 后就与其他所有变量独立，因此

p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta}) = p(\boldsymbol{y}_1|\boldsymbol{x}_1, \boldsymbol{\theta}) \ldots p(\boldsymbol{y}_n|\boldsymbol{x}_n, \boldsymbol{\theta}) p(\boldsymbol{x}|\boldsymbol{\theta})

现在我们利用完全数据的独立性假设，将 p(\boldsymbol{x}|\boldsymbol{\theta}) 展开为各个分量的乘积

重新组织这些项，我们得到

p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta}) = \prod_{i=1}^n p(\boldsymbol{y}_i|\boldsymbol{x}_i, \boldsymbol{\theta}) p(\boldsymbol{x}_i|\boldsymbol{\theta})

根据概率的乘法法则，p(\boldsymbol{y}_i|\boldsymbol{x}_i, \boldsymbol{\theta}) p(\boldsymbol{x}_i|\boldsymbol{\theta}) = p(\boldsymbol{x}_i, \boldsymbol{y}_i|\boldsymbol{\theta})，因此

p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta}) = \prod_{i=1}^n p(\boldsymbol{x}_i, \boldsymbol{y}_i|\boldsymbol{\theta})\tag{21}

它表明在我们的假设下，完全数据和观测数据的联合分布可以分解为各个样本对 (\boldsymbol{x}_i, \boldsymbol{y}_i) 的联合分布的乘积。

证明过程：推导后验分布的分解

证明的第二部分是要说明后验分布 p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}) 实际上只依赖于对应的观测 \boldsymbol{y}_i。我们从贝叶斯定理开始，后验分布等于联合分布除以边缘分布

p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}) = \frac{p(\boldsymbol{x}_i, \boldsymbol{y}|\boldsymbol{\theta})}{p(\boldsymbol{y}|\boldsymbol{\theta})}

分子是 \boldsymbol{x}_i 和所有观测 \boldsymbol{y} 的联合分布，我们需要将其他完全数据变量 \boldsymbol{x}_j, j \neq i 积分掉。这个积分在 n-1 维空间 \mathcal{X}^{n-1} 上进行

p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}) = \frac{\int_{\mathcal{X}^{n-1}} p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta}) d\boldsymbol{x}_1 \ldots d\boldsymbol{x}_{j-1}d\boldsymbol{x}_{j} d\boldsymbol{x}_{j+1} \ldots d\boldsymbol{x}_n}{\int_{\mathcal{X}^n} p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta}) d\boldsymbol{x}}

分母是对所有完全数据变量积分，得到观测数据的边缘分布。现在我们利用刚才证明的公式，将联合分布展开为乘积形式

p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}) = \frac{\int_{\mathcal{X}^{n-1}} \prod_{j=1}^n p(\boldsymbol{x}_j, \boldsymbol{y}_j|\boldsymbol{\theta}) d\boldsymbol{x}_1 \ldots d\boldsymbol{x}_{j-1} d\boldsymbol{x}_{j+1} \ldots d\boldsymbol{x}_n}{\int_{\mathcal{X}^n} \prod_{j=1}^n p(\boldsymbol{x}_j, \boldsymbol{y}_j|\boldsymbol{\theta}) d\boldsymbol{x}_1 \ldots d\boldsymbol{x}_n}

在分子中，p(\boldsymbol{x}_i, \boldsymbol{y}_i|\boldsymbol{\theta}) 不涉及其他 \boldsymbol{x}_j，所以可以提到积分外面，剩下的项是对 j \neq i 的乘积

p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}) = \frac{p(\boldsymbol{x}_i, \boldsymbol{y}_i|\boldsymbol{\theta}) \prod_{j=1, j \neq i}^n \int_{\mathcal{X}} p(\boldsymbol{x}_j, \boldsymbol{y}_j|\boldsymbol{\theta}) d\boldsymbol{x}_j}{\prod_{j=1}^n \int_{\mathcal{X}} p(\boldsymbol{x}_j, \boldsymbol{y}_j|\boldsymbol{\theta}) d\boldsymbol{x}_j}

这里我们利用了积分的可交换性，将多重积分分解为单个积分的乘积。每个积分 \int_{\mathcal{X}} p(\boldsymbol{x}_j, \boldsymbol{y}_j|\boldsymbol{\theta}) d\boldsymbol{x}_j 就是观测 \boldsymbol{y}_j 的边缘概率 p(\boldsymbol{y}_j|\boldsymbol{\theta})，因此

p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}) = \frac{p(\boldsymbol{x}_i, \boldsymbol{y}_i|\boldsymbol{\theta}) \prod_{j=1, j \neq i}^n p(\boldsymbol{y}_j|\boldsymbol{\theta})}{\prod_{j=1}^n p(\boldsymbol{y}_j|\boldsymbol{\theta})}

现在分子分母都是乘积形式，我们可以约掉那些 j \neq i 的项

p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}) = \frac{p(\boldsymbol{x}_i, \boldsymbol{y}_i|\boldsymbol{\theta})}{p(\boldsymbol{y}_i|\boldsymbol{\theta})}

根据贝叶斯定理，这正是

p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}) = p(\boldsymbol{x}_i|\boldsymbol{y}_i, \boldsymbol{\theta})\tag{22}

表明给定所有观测 \boldsymbol{y} 后 \boldsymbol{x}_i 的后验分布，实际上只依赖于对应的单个观测 \boldsymbol{y}_i。这是一个非常强的结果：尽管我们有 n 个观测，但第 i 个隐变量的后验推断并不需要用到其他观测的信息。

完成Q函数分解的证明

现在我们可以利用完全数据 X 的独立性和刚才证明的公式(22)来推导Q函数的分解。从Q函数的定义出发

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}|\boldsymbol{\theta})\}

利用完全数据的独立性，我们有 \log p(\boldsymbol{x}|\boldsymbol{\theta}) = \log \prod_{i=1}^n p(\boldsymbol{x}_i|\boldsymbol{\theta})，对数可以将乘积变为求和

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \left\{\log \prod_{i=1}^n p(\boldsymbol{x}_i|\boldsymbol{\theta})\right\}

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \left\{\sum_{i=1}^n \log p(\boldsymbol{x}_i|\boldsymbol{\theta})\right\}

期望算子对求和是线性的，所以可以将期望移到求和内部

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \sum_{i=1}^n E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}_i|\boldsymbol{\theta})\}

现在关键的一步来了：虽然期望是关于整个完全数据向量 \boldsymbol{x} 的分布 p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})，但被积函数 \log p(\boldsymbol{x}_i|\boldsymbol{\theta}) 只依赖于 \boldsymbol{x}_i。因此我们可以将这个期望简化为只关于 \boldsymbol{x}_i 的边缘分布的期望

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \sum_{i=1}^n E_{\boldsymbol{x}_i \sim p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}_i|\boldsymbol{\theta})\}

最后利用公式(22)，后验分布 p(\boldsymbol{x}_i|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) = p(\boldsymbol{x}_i|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)})，我们得到

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \sum_{i=1}^n E_{\boldsymbol{x}_i \sim p(\boldsymbol{x}_i|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}_i|\boldsymbol{\theta})\}

这正是我们要证明的公式(20)，其中每一项就是 Q_i(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)})。这个结果的实际意义是：在i.i.d.假设下，E步可以对每个样本独立进行，我们只需要计算 n 个单样本的后验分布 p(\boldsymbol{x}_i|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)})，而不需要处理维度为 n 倍的联合后验分布。这大大简化了EM算法的实现。

3.2 EM算法的收敛性分析

EM算法的单调性保证

EM算法有一个非常重要的理论性质，那就是单调性。这个性质告诉我们，改进Q函数至少不会让对数似然 l(\boldsymbol{\theta}) 变得更糟。接下来我们会通过定理3.2来严格阐述和证明这个性质。

定理3.2：EM算法的收敛性

定理的前提条件

为了陈述这个定理,我们首先需要明确几个前提条件。设随机变量 X 和 Y 具有关于参数 \boldsymbol{\theta} \in \Theta 的参数化概率密度。我们要求 X 的支撑集不依赖于参数 \boldsymbol{\theta}，也就是说，无论参数取什么值，完全数据 X 的可能取值范围都是相同的。此外，马尔可夫关系 \boldsymbol{\theta} \to X \to Y 必须成立，即

p(\boldsymbol{y}|\boldsymbol{x}, \boldsymbol{\theta}) = p(\boldsymbol{y}|\boldsymbol{x})\tag{23}

这个等式对所有 \boldsymbol{\theta} \in \Theta、\boldsymbol{x} \in \mathcal{X} 和 \boldsymbol{y} \in \mathcal{Y} 都成立。

定理的结论

在这些前提下，对于参数空间中的任意参数 \boldsymbol{\theta} \in \Theta 和任意观测 \boldsymbol{y} \in \mathcal{Y}，只要支撑集 \mathcal{X}(\boldsymbol{y}) \neq \emptyset 并且对数似然满足 l(\boldsymbol{\theta}) \geq l(\boldsymbol{\theta}^{(m)})，如果M步使得 Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) \geq Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)})，那么就有

l(\boldsymbol{\theta}) \geq Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) - Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)}) + l(\boldsymbol{\theta}^{(m)})\tag{24}

这个不等式的深刻含义我们会在证明过程中逐步揭示。

证明的整体思路

证明的基本策略是利用似然函数的定义，结合马尔可夫关系和全概率公式。我们从对数似然的定义出发

l(\boldsymbol{\theta}) = \log p(\boldsymbol{y}|\boldsymbol{\theta}) = \log \int_{\mathcal{X}(\boldsymbol{y})} p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta}) d\boldsymbol{x}

这里我们将观测数据的边缘似然表示为对完全数据的积分。注意积分域是 \mathcal{X}(\boldsymbol{y})，这是给定观测 \boldsymbol{y} 后完全数据 \boldsymbol{x} 的可能取值范围。

引入辅助分布进行重写

现在我们做一个技巧性的变换，引入上一次迭代的后验分布 p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 作为辅助分布。我们可以将积分改写为

l(\boldsymbol{\theta}) = \log \int_{\mathcal{X}(\boldsymbol{y})} \frac{p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta})}{p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) d\boldsymbol{x}

实际上我们只是乘以和除以同一个量 p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})，等式依然成立。但这样做的好处是我们现在可以把这个积分看作是一个期望

l(\boldsymbol{\theta}) = \log E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \left\{\frac{p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta})}{p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})}\right\}

这个期望是关于后验分布 p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 的，被积函数是似然比 \frac{p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta})}{p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})}。

应用Jensen不等式

接下来是证明的关键步骤。由于对数函数是凹函数，我们可以应用Jensen不等式。Jensen不等式告诉我们，对于凹函数，函数的期望小于等于期望的函数，即对于凹函数 f 和随机变量 X，有 E[f(X)] \leq f(E[X])。但对于我们这里的情况，由于对数在外面，我们有

l(\boldsymbol{\theta}) \geq E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \left\{\log \frac{p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta})}{p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})}\right\}

这个不等式将对数从期望外面移到了里面，代价是等号变成了不等号。Jensen不等式在这里的应用是整个证明的核心，它将似然函数与Q函数联系了起来。

利用贝叶斯定理和马尔可夫性质

现在我们来简化右边的期望。首先处理分子，利用贝叶斯定理的链式法则，我们有 p(\boldsymbol{x}, \boldsymbol{y}|\boldsymbol{\theta}) = p(\boldsymbol{x}|\boldsymbol{\theta})p(\boldsymbol{y}|\boldsymbol{x})。这里关键的一点是，根据马尔可夫关系公式(23)，p(\boldsymbol{y}|\boldsymbol{x}, \boldsymbol{\theta}) = p(\boldsymbol{y}|\boldsymbol{x})，所以 p(\boldsymbol{y}|\boldsymbol{x}) 不依赖于参数。因此

l(\boldsymbol{\theta}) \geq E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \left\{\log \frac{p(\boldsymbol{x}|\boldsymbol{\theta})p(\boldsymbol{y}|\boldsymbol{x})}{p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})}\right\}

接下来我们处理分母。对于分母，我们也可以用贝叶斯定理来展开。根据贝叶斯定理

p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) = \frac{p(\boldsymbol{y}|\boldsymbol{x})p(\boldsymbol{x}|\boldsymbol{\theta}^{(m)})}{p(\boldsymbol{y}|\boldsymbol{\theta}^{(m)})}

这里再次用到了马尔可夫关系，p(\boldsymbol{y}|\boldsymbol{x}, \boldsymbol{\theta}^{(m)}) = p(\boldsymbol{y}|\boldsymbol{x})。

代入后得到

化简后，p(\boldsymbol{y}|\boldsymbol{x}) 约掉了，我们得到

l(\boldsymbol{\theta}) \geq E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \left\{\log \frac{p(\boldsymbol{x}|\boldsymbol{\theta})p(\boldsymbol{y}|\boldsymbol{\theta}^{(m)})}{p(\boldsymbol{x}|\boldsymbol{\theta}^{(m)})}\right\}

分解对数并展开期望

现在我们利用对数的性质，将分数的对数分解为两个对数之差

l(\boldsymbol{\theta}) \geq E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \left\{\log p(\boldsymbol{x}|\boldsymbol{\theta}) + \log p(\boldsymbol{y}|\boldsymbol{\theta}^{(m)}) - \log p(\boldsymbol{x}|\boldsymbol{\theta}^{(m)})\right\}

由于期望算子是线性的，我们可以将这个期望拆分为三项

l(\boldsymbol{\theta}) \geq E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}|\boldsymbol{\theta})\} - E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}|\boldsymbol{\theta}^{(m)})\} + \log p(\boldsymbol{y}|\boldsymbol{\theta}^{(m)})

注意最后一项 \log p(\boldsymbol{y}|\boldsymbol{\theta}^{(m)}) 不依赖于 \boldsymbol{x}，所以可以直接从期望中提出来。现在观察前两项，第一项正是Q函数的定义

E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}|\boldsymbol{\theta})\} = Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)})

第二项是将 \boldsymbol{\theta} 替换为 \boldsymbol{\theta}^{(m)}

E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}|\boldsymbol{\theta}^{(m)})\} = Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)})

而第三项就是上一次迭代的对数似然

\log p(\boldsymbol{y}|\boldsymbol{\theta}^{(m)}) = l(\boldsymbol{\theta}^{(m)})

将这三项代入，我们得到

l(\boldsymbol{\theta}) \geq Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) - Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)}) + l(\boldsymbol{\theta}^{(m)})

这正是公式(24)，它建立了对数似然 l(\boldsymbol{\theta}) 与Q函数之间的下界关系。

推导单调性结论

现在我们来看这个不等式的实际含义。在M步中，我们选择 \boldsymbol{\theta} 来最大化 Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)})，这意味着我们找到的新参数 \boldsymbol{\theta}^{(m+1)} 满足

Q(\boldsymbol{\theta}^{(m+1)}|\boldsymbol{\theta}^{(m)}) \geq Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)})

这是因为 \boldsymbol{\theta}^{(m+1)} 是最大化Q函数的点，所以Q函数在这一点的值至少不会小于在 \boldsymbol{\theta}^{(m)} 处的值。将这个结果代入公式(24)，我们得到

l(\boldsymbol{\theta}^{(m+1)}) \geq Q(\boldsymbol{\theta}^{(m+1)}|\boldsymbol{\theta}^{(m)}) - Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)}) + l(\boldsymbol{\theta}^{(m)})\tag{25}

由于 Q(\boldsymbol{\theta}^{(m+1)}|\boldsymbol{\theta}^{(m)}) - Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)}) \geq 0，我们立即得到

l(\boldsymbol{\theta}^{(m+1)}) \geq l(\boldsymbol{\theta}^{(m)})\tag{25}

这正是公式(25)所表达的单调性：每次EM迭代后，对数似然不会下降。这就完成了证明。这个结果告诉我们，EM算法保证每次迭代都不会让似然值变差。

注释3.3：单调性的深层含义

单调性与收敛性的区别

EM算法的单调性保证了似然值序列 \{l(\boldsymbol{\theta}^{(m)})\} 是单调不降的，也就是说，算法产生的猜测在似然意义下不会变糟。然而，单调性本身并不能保证参数序列 \{\boldsymbol{\theta}^{(m)}\} 一定会收敛。似然值不降并不意味着参数一定会稳定到某个固定点，理论上参数可能会在保持似然值不变的情况下振荡。

有界性对收敛性的作用

但是如果对数似然函数 l(\boldsymbol{\theta}) 在参数空间 \Theta 上有上界，那么情况就不同了。在这种情况下，单调不降的似然值序列 \{l(\boldsymbol{\theta}^{(m)})\} 必然会收敛到某个极限值。虽然这保证了似然值序列的收敛，但仍然不直接保证参数序列 \{\boldsymbol{\theta}^{(m)}\} 的收敛。不过在实践中，如果似然值已经稳定，参数通常也会接近某个稳定的解。

注释3.4：在EM算法中引入先验

为何需要引入先验

正如我们之前提到的，在某些情况下引入先验知识是有益的，甚至是必要的。当我们拥有关于参数的先验信息时，或者为了避免对数似然函数的奇异性问题时，我们可以在EM框架中引入先验信息。这通过使用MAP估计来实现

这里 \log p(\boldsymbol{\theta}) 是参数的对数先验概率。加入这一项相当于对参数施加了正则化，可以防止过拟合，同时也能在数据不足时利用我们的先验知识。

MAP-EM算法的形式

将MAP的思想融入EM算法。在E步，我们仍然给定当前参数估计 \boldsymbol{\theta}^{(m)}，计算Q函数

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = E_{\boldsymbol{x} \sim p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{x}|\boldsymbol{\theta})\}\tag{27}

在M步，最大化步骤需要修改为

\boldsymbol{\theta}^{(m+1)} = \arg\min_{\boldsymbol{\theta} \in \Theta} \left(Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) + \log p(\boldsymbol{\theta})\right)\tag{28}

这里我们在原来的Q函数基础上加上了对数先验项 \log p(\boldsymbol{\theta})。这个修改后的目标函数同时考虑了数据拟合（通过Q函数）和先验约束（通过先验项）。

MAP-EM的单调性

MAP-EM算法的单调性可以用与定理3.2完全相同的方式来证明。唯一的区别是我们优化的目标从 l(\boldsymbol{\theta}) 变成了 l(\boldsymbol{\theta}) + \log p(\boldsymbol{\theta})，也就是从最大化似然变成了最大化后验。证明的每一步推导都基本类似，最后单调性保证：MAP-EM每次迭代都会使得后验概率不降。

3.3 通过EM算法学习隐变量模型

隐变量模型的基本框架

沿用相同的方法论，EM算法可以被用来估计包含隐变量的有向模型的参数。这些隐变量用 \boldsymbol{z} 来表示，在我们实际可获得的数据中并不能直接观测到他们。当我们对观测变量 \boldsymbol{y} 进行无条件建模时，相应的有向图模型描述的是一个联合概率分布 p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta})，这个分布同时涉及观测数据和隐藏变量。为了得到仅关于观测变量的分布，我们需要对隐变量进行边缘化，即将隐变量积分掉

p(\boldsymbol{y}|\boldsymbol{\theta}) = \int p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta}) d\boldsymbol{z}\tag{29}

这个积分得到的 p(\boldsymbol{y}|\boldsymbol{\theta}) 在MLE估计的语境下被称为边缘似然，而在生成模型的语境下被称为模型证据。当我们把它看作参数 \boldsymbol{\theta} 的函数时，它就是我们要最大化的目标。

隐变量模型中的完全数据定义

在这种情况下，完全数据自然地被定义为观测数据和隐变量的组合，即 \boldsymbol{x} = \{\boldsymbol{y}, \boldsymbol{z}\}。有了这个定义，我们就可以用EM算法来学习参数 \boldsymbol{\theta}。E步的形式需要根据这个完全数据的定义来展开。我们从Q函数的一般定义出发

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \int \log p(\boldsymbol{x}|\boldsymbol{\theta}) p(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) d\boldsymbol{x}

将完全数据 \boldsymbol{x} = \{\boldsymbol{y}, \boldsymbol{z}\} 代入，我们需要对 \boldsymbol{x} 的积分实际上就是对 \boldsymbol{z} 的积分，因为 \boldsymbol{y} 是已知的观测。因此

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \int \log p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta}) p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) d\boldsymbol{x}

这里 p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 看起来有些冗余，因为已经给定了 \boldsymbol{y}，这实际上就是 p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})。所以我们可以改写为

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \int \log p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta}) p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) d\boldsymbol{z}

用期望的形式来表达

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = E_{\boldsymbol{z} \sim p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta})\}

最后，由于 \boldsymbol{z} 是唯一的随机变量，我们也可以写成

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = E_{\boldsymbol{z} \sim p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta})\}\tag{30}

这就是公式(30)，它告诉我们在隐变量模型中，E步需要计算的是联合对数概率 \log p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta}) 关于隐变量后验分布 p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 的期望。

深度隐变量模型的表达能力

在深度隐变量模型的语境下，涉及的分布 p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta}) = p(\boldsymbol{z}|\boldsymbol{\theta})p(\boldsymbol{y}|\boldsymbol{z}, \boldsymbol{\theta}) 通常由神经网络参数化。即使我们假设相对简单的分布形式（比如参数化的高斯分布），边缘分布 p(\boldsymbol{y}|\boldsymbol{\theta}) 仍然非常复杂。这种表达能力来自于隐变量的引入和神经网络的非线性变换。通过对隐变量积分（或求和），简单的条件分布可以组合成复杂的边缘分布。这种表达性使得深度隐变量模型特别适合用来近似那些复杂的真实数据分布 p^*(\boldsymbol{y})。我们可以用相对简单的构建模块（简单的先验和似然），通过隐变量的层次结构和神经网络的非线性映射，来捕捉数据中复杂的依赖关系和结构。

3.4 学习固定模型的最优混合

问题的具体设定

现在我们通过一个例子来说明EM算法的应用。

假设我们观测到数据 \boldsymbol{y} = [\boldsymbol{y}_1, \ldots, \boldsymbol{y}_n]^\top，这些数据可能是由 k 个模型 (p_1, \ldots, p_k) 生成的。举个例子，在高斯混合的场景中，p_i 表示一个已知均值和协方差矩阵的高斯概率密度函数。我们进一步假设观测是独立同分布的，来自这 k 个模型的某种组合，具体的混合分布形式为

p(\boldsymbol{y}_i|\boldsymbol{\theta}) = \sum_{j=1}^k \theta_j p_j(\boldsymbol{y}_i)\tag{31}

这里的约束条件是 \sum_{j=1}^k \theta_j = 1 并且 \theta_j \in [0,1] 对所有 j 成立。我们的目标是学习出这个最可能的混合系数组合 \boldsymbol{\theta} = [\theta_1, \ldots, \theta_k]^\top，使得它最好地刻画了观测数据 \boldsymbol{y}。

完全数据的选择与直觉

从之前例2.2的讨论我们知道，如果我们知道每个观测是由哪个模型生成的，问题会大大简化。因此，很自然地，我们选择EM的完全数据为观测数据加上一组隐藏数据 \boldsymbol{z} = [z_1, \ldots, z_n]^\top，其中 z_i 表示第 i 个观测是由 k 个模型中的哪一个生成的。具体来说，z_i \in \{1, \ldots, k\}，对于 i \in \{1, \ldots, n\}。这样一来，给定了隐变量 z_i = j 的情况下，第 i 个观测的条件概率就非常简单

p(\boldsymbol{y}_i, z_i = j|\boldsymbol{\theta}) = \theta_j p_j(\boldsymbol{y}_i)\tag{32}

这个公式表明，如果我们知道第 i 个观测来自第 j 个分量，那么它的联合概率就是混合系数 \theta_j 乘以第 j 个模型的概率密度 p_j(\boldsymbol{y}_i)。

E步：计算后验概率

根据公式(32)，我们可以很容易地推导出隐变量的后验概率。利用贝叶斯定理

P(z_i = j|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)}) = \frac{p(z_i = j, \boldsymbol{y}_i|\boldsymbol{\theta}^{(m)})}{p(\boldsymbol{y}_i|\boldsymbol{\theta}^{(m)})}\tag{33}

分子是联合概率，根据公式(32)就是 \theta_j^{(m)} p_j(\boldsymbol{y}_i)。分母是边缘概率，需要对所有可能的 z_i 求和，即 p(\boldsymbol{y}_i|\boldsymbol{\theta}^{(m)}) = \sum_{l=1}^k \theta_l^{(m)} p_l(\boldsymbol{y}_i)。因此

P(z_i = j|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)}) = \frac{\theta_j^{(m)} p_j(\boldsymbol{y}_i)}{\sum_{l=1}^k \theta_l^{(m)} p_l(\boldsymbol{y}_i)}\tag{33}

这就是公式(33)，它给出了在当前参数估计下，第 i 个观测属于第 j 个分量的后验概率。这个后验概率也被称为"责任度"，它量化了第 j 个分量对生成第 i 个观测的责任程度。

M步：推导更新公式

接下来是M步，我们需要最大化Q函数。根据隐变量模型的EM形式，M步变成

\boldsymbol{\theta}^{(m)} = \arg\max_{\boldsymbol{\theta} \in \Theta} E_{\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}} \{\log p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta})\}

这是使用缺失数据形式的EM公式。由于数据是独立同分布的，我们可以利用i.i.d.情况下的EM形式，将期望分解为对各个样本的求和

\boldsymbol{\theta}^{(m)} = \arg\max_{\boldsymbol{\theta} \in \Theta} \sum_{i=1}^n \sum_{j=1}^k p(z_i = j|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)}) \{\log p(\boldsymbol{y}_i, z_i = j|\boldsymbol{\theta})\}

现在我们展开 \log p(\boldsymbol{y}_i, z_i = j|\boldsymbol{\theta})。根据公式(31)中的混合模型定义和公式(32)

\boldsymbol{\theta}^{(m)} = \arg\max_{\boldsymbol{\theta} \in \Theta} \sum_{i=1}^n \sum_{j=1}^k p(z_i = j|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)}) \{\log \theta_{z_i=j} + \log p_{z_i=j}(\boldsymbol{y}_i)\}

\log p_{z_i=j}(\boldsymbol{y}_i) 不依赖于我们要优化的参数 \boldsymbol{\theta}（因为模型 p_j 是已知固定的），所以在优化过程中可以当作常数忽略。因此优化问题简化为

\boldsymbol{\theta}^{(m)} = \arg\max_{\boldsymbol{\theta} \in \Theta} \sum_{i=1}^n \sum_{j=1}^k p(z_i = j|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)}) \log \theta_j

为了简化记号，我们定义一个辅助量 \alpha_j^{(m)}，它是所有样本对第 j 个分量的后验概率之和

\alpha_j^{(m)} = \sum_{i=1}^n p(z_i = j|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)})

将公式(33)代入，我们得到

\alpha_j^{(m)} = \sum_{i=1}^n \frac{\theta_j^{(m)} p_j(\boldsymbol{y}_i)}{\sum_{l=1}^k \theta_l^{(m)} p_l(\boldsymbol{y}_i)}

有了这个定义，优化问题就变成了

\boldsymbol{\theta}^{(m)} = \arg\max_{\boldsymbol{\theta} \in \Theta} \sum_{j=1}^k \alpha_j^{(m)} \log \theta_j

使用拉格朗日乘数法处理约束

现在我们面临一个带约束的优化问题。我们需要最大化 \sum_{j=1}^k \alpha_j^{(m)} \log \theta_j，同时满足约束 \sum_j \theta_j = 1 和 \theta_j \in [0,1] 对所有 j 成立。求解这个约束优化问题可以用拉格朗日乘数法来处理。不过在这里，我们暂时忽略不等式约束 \theta_j \in [0,1]，只考虑等式约束 \sum_j \theta_j = 1。构造拉格朗日函数

\mathcal{L} = \sum_{l=1}^k \alpha_l^{(m)} \log \theta_l - \lambda \left(\sum_{l=1}^k \theta_l - 1\right)

对 \theta_j 求偏导并令其为零

\frac{\partial}{\partial \theta_j} \left(\sum_{l=1}^k \alpha_l^{(m)} \log \theta_l - \lambda \left(\sum_{l=1}^k \theta_l - 1\right)\right) = \frac{\alpha_j^{(m)}}{\theta_j} - \lambda = 0

这给出 \theta_j^* = \frac{\alpha_j^{(m)}}{\lambda}。现在我们需要确定拉格朗日乘数 \lambda。将这个解代入约束条件 \sum_j \theta_j = 1，我们得到

\sum_j \frac{\alpha_j^{(m)}}{\lambda} = 1

因此 \lambda = \sum_{l=1}^k \alpha_l^{(m)}。这样最终的解就是

\theta_j^* = \frac{\alpha_j^{(m)}}{\sum_{l=1}^k \alpha_l^{(m)}}

公式(33)中定义的优化问题实际上是一个凸优化问题，这个解确实是全局最优的(不再验证解的最优性了)。

EM算法的完整流程

综合E步和M步，我们得到学习固定模型最优混合的完整EM算法。在初始化阶段，我们选择一个初始参数估计 \boldsymbol{\theta}^{(0)}。

在E步，我们计算后验概率

\alpha_j^{(m)} = p(z_i = j|\boldsymbol{y}_i, \boldsymbol{\theta}^{(m)}) = \sum_{i=1}^n \frac{\theta_j^{(m)} p_j(\boldsymbol{y}_i)}{\sum_{l=1}^k \theta_l^{(m)} p_l(\boldsymbol{y}_i)}, \quad j = 1, \ldots, k

这一步对每个分量 j 和每个样本 i 计算责任度，然后求和得到 \alpha_j^{(m)}。

在M步，我们更新参数

\theta_j^{(m+1)} = \frac{\alpha_j^{(m)}}{\sum_{l=1}^k \alpha_l^{(m)}}, \quad j = 1, \ldots, k

新的混合系数正比于该分量对所有样本的总责任度。

最后是收敛性检查。我们计算新参数下的对数似然 l(\boldsymbol{\theta}^{(m+1)})，如果似然的改进量 |l(\boldsymbol{\theta}^{(m+1)}) - l(\boldsymbol{\theta}^{(m)})| 大于某个预设的阈值 \delta，我们就返回E步继续迭代；否则算法终止。这个停止准则保证了我们在似然值基本不再改善时停止迭代，避免不必要的计算。整个算法通过E步和M步的交替进行，单调地提升似然值，最终收敛到一个局部最优解。

4 超越标准EM算法

标准EM算法的优势与挑战

EM算法专门用于处理不完整数据或隐变量数据时的极大似然估计。它的迭代过程在估计缺失信息和优化模型参数之间交替进行，将一个困难的优化问题分解成了两个相对简单的步骤

然而，尽管标准EM算法有这些优势，它在实践中仍然面临一些重要的局限性。首先，EM算法可能收敛得很慢，特别是在似然表面的平坦区域附近。当算法接近某个最优解时，如果该区域的似然函数梯度很小，参数的更新步长就会变得非常微小，导致需要大量迭代才能达到收敛。

其次，EM算法对初始参数值往往很敏感，不同的初始值可能导致收敛到不同的局部最优解。如果初始值选择不当，算法可能会陷入一个很差的局部最优，得到远离全局最优的解。

第三个也是最关键的问题是，在许多实际应用中，特别是那些涉及复杂隐变量结构或大规模数据集的问题中，期望步骤或最大化步骤在计算上变得不可处理。E步可能因为需要对高维隐变量空间进行积分而变得不可行，M步可能因为优化目标函数过于复杂而无法求解。

这些局限性让人开发了EM算法的各种变体和扩展。这些改进方法的设计目标是提高收敛速度、降低计算成本，或者处理精确更新不可能实现的情况。在本节中，我们将看到如何在保持EM核心框架的基础结构的同时，对其进行调整和扩展来应对这些挑战。这些变体虽然在具体实现上有所不同，但都通过引入和处理隐变量来简化原本困难的优化问题。

4.1 处理不可处理的E步

随机方法的基本思想

当E步在计算上不可处理时（复杂的概率密度函数或非线性设置），一个自然的选择是使用EM算法的随机变体。核心思想是：与其计算完整的期望，不如使用随机样本通过蒙特卡罗方法来近似这个期望。

随机EM算法

随机EM算法用基于样本的近似来替代完整的E步。具体来说，SEM不是计算期望 E_{\boldsymbol{x} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})}[\log f(\boldsymbol{x}|\boldsymbol{\theta})]，而是从隐变量的条件分布中抽取一个或几个随机样本，然后直接使用这些样本来近似期望。在最简单的形式中，SEM在每次迭代只使用一个样本。假设我们从后验分布 f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 中采样得到 \boldsymbol{x}^{(m)}，那么Q函数的近似就变成

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) \approx \log f(\boldsymbol{x}^{(m)}|\boldsymbol{\theta})

其中 \boldsymbol{x}^{(m)} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})

这个近似虽然是有偏的（单个样本不能准确代表整个分布的期望），但随着迭代的进行，在一定条件下算法仍然可以收敛。SEM特别适用于处理大规模或复杂数据集，在这些情况下计算精确期望是不现实的。通过牺牲每一步的精确性，我们换取了计算上的可行性。

蒙特卡罗EM算法

蒙特卡罗EM算法通过使用蒙特卡罗采样来近似E步，在精度上比SEM有所改进。MCEM不是只用一个样本，而是使用多个样本来估计期望。假设我们从后验分布中抽取 K 个独立样本 \boldsymbol{x}^{(m),1}, \boldsymbol{x}^{(m),2}, \ldots, \boldsymbol{x}^{(m),K}，那么Q函数可以近似为

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) \approx \frac{1}{K} \sum_{c=1}^K \log f(\boldsymbol{x}^{(m),c}|\boldsymbol{\theta})

其中每个 \boldsymbol{x}^{(m),c} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})

根据大数定律，当样本数 K 足够大时，这个样本平均会收敛到真实的期望。MCEM对期望的估计更加准确，因此在M步中使用的目标函数更接近真实的Q函数。这在解析解不可用但采样是可行的情况下特别有用。MCEM的代价是每次迭代需要生成多个样本，计算量比SEM大，但得到的估计质量更高。实践中，样本数 K 可以根据计算资源和精度要求进行调整，甚至可以在迭代过程中动态调整样本数。

随机近似EM算法

随机近似EM算法在随机更新的基础上进一步引入了平滑机制，通过增量式地更新期望来提高稳定性和收敛性。SAEM的核心思想是不在每一步完全替换Q函数的估计，而是将当前估计与历史估计进行加权组合。具体来说，SAEM使用一个递归更新公式

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m+1)}) = Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) + \gamma_k \left[\log f(\boldsymbol{x}^{(m)}|\boldsymbol{\theta}) - Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)})\right]

这里 \boldsymbol{x}^{(m)} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 是从当前后验分布采样得到的样本，\gamma_k 是一个介于0和1之间的步长参数，满足 0 < \gamma_k < 1。这个更新规则可以理解为一个移动平均：新的Q函数估计是旧估计和当前样本贡献的加权平均。步长参数 \gamma_k 通常随着迭代次数递减，这样早期迭代中算法对新信息的响应更快，而后期迭代中估计更加稳定。这种平滑机制使得SAEM即使在每次迭代只使用少量样本的情况下也能获得较好的稳定性和收敛性。SAEM特别适合处理那些采样成本较高但需要稳定收敛的问题。

EM变体的比较总结

标准EM算法计算精确的条件期望，其E步描述为：计算完全数据对数似然的精确条件期望，对应的解析方程是

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) = \mathbb{E}_{\boldsymbol{x} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})}[\log f(\boldsymbol{x}|\boldsymbol{\theta})]

随机EM则使用来自后验分布的单个样本来近似E步，其近似形式为

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) \approx \log f(\boldsymbol{x}^{(m)}|\boldsymbol{\theta})

其中 \boldsymbol{x}^{(m)} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})。

蒙特卡罗EM使用多个样本来估计E步中的期望，其近似为

Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) \approx \frac{1}{K} \sum_{c=1}^K \log f(\boldsymbol{x}^{(m),c}|\boldsymbol{\theta})

其中每个 \boldsymbol{x}^{(m),c} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})。

随机近似EM则使用递归更新来结合当前和过去的估计，其更新公式为

其中 \boldsymbol{x}^{(m)} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})，步长满足 0 < \gamma_k < 1。

4.2 处理不可处理的M步

MAP框架下的M步复杂性

在MAP估计的框架下，EM算法的E步保持不变，期望仍然只是关于隐变量的。但M步会发生变化，因为我们现在要最大化的不仅是Q函数，还要加上先验的对数概率。具体来说，MAP情况下的M步变成

\boldsymbol{\theta}_{\text{MAP}}^{(k+1)} = \arg\max_{\boldsymbol{\theta} \in \Omega_{\theta}} Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(m)}) + \log f(\boldsymbol{\theta})

这里 f(\boldsymbol{\theta}) 表示参数的先验分布。由于先验的存在，即使原本的Q函数优化是可处理的，加上先验项后M步也可能变得难以计算。先验可能引入额外的约束或复杂的惩罚项，使得联合优化问题没有闭式解，或者求解的计算成本过高。因此我们需要一些变体算法来应对M步不可处理的情况。

EM变体的核心思想

在经典的EM算法中，M步涉及对所有参数进行完全的联合最大化，这在复杂模型中可能在计算上非常昂贵甚至不可行。为了克服这个困难，提出了几种策略。

期望条件最大化算法

期望条件最大化算法通过将M步分解为一系列更简单的条件最大化来解决复杂性问题。ECM不是一次性对所有参数进行优化，而是将参数分成若干个块，然后依次对每个参数块进行条件最大化，同时保持其他参数块固定。假设我们将参数向量 \boldsymbol{\theta} 分成 B 个块 \{\boldsymbol{\theta}_b\}_{b=1}^B，那么ECM在每次迭代中循环执行以下步骤：对于每个参数块 b，更新

\boldsymbol{\theta}_b^{(m+1)} = \arg\max_{\boldsymbol{\theta}_b \in \Theta_b} \mathbb{E}_{\boldsymbol{x} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}_1^{(m)}, \boldsymbol{\theta}_{\bar{b}})}[\log f(\boldsymbol{x}|\boldsymbol{\theta}_b, \boldsymbol{\theta}_{\bar{b}})]

这里 \boldsymbol{\theta}_{\bar{b}} 表示除了第 b 块之外的所有其他参数，在优化第 b 块时这些参数保持固定。这种分块优化策略将一个高维的联合优化问题分解为多个低维的条件优化问题，每个子问题通常更容易求解。

举例:

要估计2分量高斯混合模型的参数：

\boldsymbol{\theta} = \{\underbrace{\theta_1, \theta_2}_{\text{块1: 混合系数}}, \underbrace{\boldsymbol{\mu}_1, \boldsymbol{\mu}_2}_{\text{块2: 均值}}, \underbrace{\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2}_{\text{块3: 协方差}}\}

分成3个块，依次优化：

第1步：固定 \boldsymbol{\mu}, \boldsymbol{\Sigma}，只优化混合系数

\{\theta_1^{(m+1)}, \theta_2^{(m+1)}\} = \arg\max_{\theta_1, \theta_2} Q(\theta_1, \theta_2, \boldsymbol{\mu}^{(m)}, \boldsymbol{\Sigma}^{(m)}|\boldsymbol{\theta}^{(m)})

得到：

\theta_j^{(m+1)} = \frac{1}{n}\sum_{i=1}^{n} \gamma_{ij}^{(m)}

第2步：固定 \theta^{(m+1)}, \boldsymbol{\Sigma}^{(m)}，只优化均值

\{\boldsymbol{\mu}_1^{(m+1)}, \boldsymbol{\mu}_2^{(m+1)}\} = \arg\max_{\boldsymbol{\mu}_1, \boldsymbol{\mu}_2} Q(\theta^{(m+1)}, \boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \boldsymbol{\Sigma}^{(m)}|\boldsymbol{\theta}^{(m)})

得到：

\boldsymbol{\mu}_j^{(m+1)} = \frac{\sum_{i=1}^{n} \gamma_{ij}^{(m)} \boldsymbol{y}_i}{\sum_{i=1}^{n} \gamma_{ij}^{(m)}}

第3步：固定 \theta^{(m+1)}, \boldsymbol{\mu}^{(m+1)}，只优化协方差

\{\boldsymbol{\Sigma}_1^{(m+1)}, \boldsymbol{\Sigma}_2^{(m+1)}\} = \arg\max_{\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2} Q(\theta^{(m+1)}, \boldsymbol{\mu}^{(m+1)}, \boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2|\boldsymbol{\theta}^{(m)})

期望条件最大化二选一算法

期望条件最大化二选一算法在ECM的基础上进一步改进，通过在某些步骤中直接最大化观测数据的对数似然来加快收敛。ECME的思想是：对于某些参数（记为 \boldsymbol{\theta}_{\text{OL}}，OL代表observed likelihood），直接最大化观测数据的似然 \log f(\boldsymbol{y}|\boldsymbol{\theta}_{\text{OL}}, \boldsymbol{\theta}_{\overline{\text{OL}}})，当然这取决于参数的性质条件，比如说有些参数是对角矩阵，那么就有简单的闭式解，那么就可以直接最大化似然来求解

具体来说，ECME的更新分为两类步骤。对于那些使用观测似然的参数

\boldsymbol{\theta}_{\text{OL}}^{(m+1)} = \arg\max_{\boldsymbol{\theta}_{\text{OL}} \in \Theta_{\text{OL}}} \log f(\boldsymbol{y}|\boldsymbol{\theta}_{\text{OL}}, \boldsymbol{\theta}_{\overline{\text{OL}}})

而对于其他参数，仍然使用标准的Q函数最大化

\boldsymbol{\theta}_{\overline{\text{OL}}}^{(m+1)} = \arg\max_{\boldsymbol{\theta}_{\overline{\text{OL}}} \in \Theta_{\overline{\text{OL}}}} \mathbb{E}_{\boldsymbol{x} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}_{\text{OL}}^{(m)}, \boldsymbol{\theta}_{\overline{\text{OL}}}^{(m)})}[\log f(\boldsymbol{x}|\boldsymbol{\theta}_{\text{OL}}, \boldsymbol{\theta}_{\overline{\text{OL}}})]

这种混合策略结合了两种优化方式的优势：对于那些直接优化观测似然更有效的参数使用观测似然，对于其他参数使用完全数据似然。这往往能够在保持单调性的同时获得更快的收敛速度。

广义EM算法

广义EM算法通过完全放松最大化的要求来提供最大的灵活性。GEM不要求M步必须找到Q函数的最大值，只要求新的参数能够使Q函数相比上一次迭代有所增加即可。也就是说，我们只需要找到 \boldsymbol{\theta}^{(m+1)} 使得

Q(\boldsymbol{\theta}^{(m+1)}|\boldsymbol{\theta}^{(m)}) \geq Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)})

这个条件比完全最大化弱得多，因此给了我们更多的选择空间。我们可以使用部分更新、基于梯度的方法，或者任何能够改进Q函数的优化策略。虽然GEM的每一步改进可能不如完全最大化那么大，但它保证了单调性，并且在很多实际问题中仍然能够收敛到好的解。

M步变体的系统比较

为了更系统地理解这些变体，我们可以从它们的M步公式来对比。

标准EM算法的M步描述是：最大化关于所有参数的完全数据对数似然的期望，其解析方程为

\boldsymbol{\theta}^{(m+1)} = \arg\max_{\boldsymbol{\theta} \in \Theta} \mathbb{E}_{\boldsymbol{x} \sim f(\boldsymbol{x}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})}[\log f(\boldsymbol{x}|\boldsymbol{\theta})]

这是对所有参数的联合优化。

ECM的M步则被替换为对参数子集的条件最大化序列。对于每个参数块 b，更新公式为

这里关键是每次只优化一个参数块，同时其他参数块固定。

ECME的M步分为两类。对于某些参数 \boldsymbol{\theta}_{\text{OL}}，直接使用观测数据对数似然

而对于其他参数 \boldsymbol{\theta}_{\overline{\text{OL}}}，仍然使用期望形式

这种混合方式在实践中往往能带来计算效率和收敛速度的显著提升。

GEM的M步要求最为宽松，只需要找到满足条件的 \boldsymbol{\theta}^{(m+1)} 使得

Q(\boldsymbol{\theta}^{(m+1)}|\boldsymbol{\theta}^{(m)}) \geq Q(\boldsymbol{\theta}^{(m)}|\boldsymbol{\theta}^{(m)})

这个条件不要求必须找到最大值，只要求函数值有所改进。这为各种近似优化方法打开了大门，使得算法在复杂情况下仍然可行。

4.3 变分EM算法

blabla

5 变分自编码器

EM算法的局限性与变分推断的崛起

EM算法长期以来一直是估计带有隐变量的概率模型参数的核心工具。对于许多结构化问题，经典的EM算法及其随机变体都能有效工作。然而在现代机器学习应用中，这些方法遇到了显著的可扩展性问题。具体来说，即使是EM的随机版本，它们依赖采样方法来近似期望，在大规模数据集或深度生成模型上也难以高效扩展。当数据量巨大或模型复杂度很高时，采样本身就成为了计算瓶颈。

为了克服这些挑战，变分推断提供了一个强大的替代方案。变分推断的核心思想是用一个可处理分布族中的近似分布来替代精确但难以计算的后验分布。这个近似分布通过优化技术（比如证据下界ELBO）来获得。变分自编码器正是将这种变分推断思想与神经网络结合起来，使得我们能够高效地学习复杂的生成模型。

VAE的本质可以这样理解：EM算法依赖于迭代地计算或采样不可处理的后验分布，而VAE则通过学习到的编码器网络来进行推断，这使得它们特别适合处理高维和大规模问题，在这些场景中经典EM算法会失效。

5.1 深度隐变量模型

隐变量模型的结构回顾

正如我们之前解释的，隐变量（用 \boldsymbol{z} 表示）可以被纳入到模型的结构中，但它们在我们实际可获得的数据中并不能直接观测到。当我们对观测变量 \boldsymbol{y} 进行无条件建模时，相应的有向图模型描述的是一个联合概率分布 p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta})，这个分布同时涉及观测数据和隐藏变量。为了得到仅关于观测变量的分布，我们需要对隐变量进行边缘化，通过积分将隐变量消去

p^*(\boldsymbol{y}) \approx p(\boldsymbol{y}|\boldsymbol{\theta}) = \int p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta}) d\boldsymbol{z}\tag{34}

深度隐变量模型的动机

当真实的数据分布 p^*(\boldsymbol{y}) 很复杂时，直接对 p(\boldsymbol{y}|\boldsymbol{\theta}) 建模往往会导致EM算法变得不可处理。这是因为复杂的边缘分布通常没有好的解析形式，后验分布 p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}) 也无法精确计算或有效采样。在这种情况下，一个常见的策略是假设深度隐变量模型分布。深度隐变量模型的关键特点是用神经网络来参数化联合分布 p(\boldsymbol{y}, \boldsymbol{z}|\boldsymbol{\theta}) = p(\boldsymbol{z}|\boldsymbol{\theta})p(\boldsymbol{y}|\boldsymbol{z}, \boldsymbol{\theta}) 中涉及的各个条件分布。

这种参数化方式带来了巨大的表达能力：即使我们为先验 p(\boldsymbol{z}|\boldsymbol{\theta}) 和似然 p(\boldsymbol{y}|\boldsymbol{z}, \boldsymbol{\theta}) 选择相对简单的形式（比如高斯分布），通过神经网络的非线性变换和隐变量的边缘化，得到的边缘分布 p(\boldsymbol{y}|\boldsymbol{\theta}) 仍然可以非常复杂。

变分推断框架的引入

为了在深度隐变量模型中获得可处理的推断和学习过程，我们引入一个参数化的推断模型 q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) 来近似真实的后验分布 p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta})。这个近似后验被称为编码器，而 \boldsymbol{\phi} 代表变分参数，包括神经网络的权重和偏置。在实践中，我们通常让编码器输出后验分布的参数。例如，如果我们假设后验是高斯分布，编码器网络可以输出均值和协方差

\boldsymbol{\mu}, \boldsymbol{\Sigma} = \text{ENN}(\boldsymbol{y}, \boldsymbol{\phi})

q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = \mathcal{N}(\boldsymbol{z}; \boldsymbol{\mu}, \boldsymbol{\Sigma})

这里ENN表示编码器神经网络。给定观测数据 \boldsymbol{y}，编码器网络输出一个高斯分布的参数，这个分布近似了隐变量的后验。这种参数化方式的优势在于，一旦网络训练完成，我们可以通过一次前向传播快速获得任意观测的后验近似，而不需要像EM那样每次都进行迭代计算或采样。

5.2 将MLE与ELBO联系起来

ELBO的定义与重要性

在深度隐变量模型的语境下，证据下界通常被选作优化目标。ELBO的定义为

\text{(ELBO)} \quad \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) = E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \left\{\log \frac{p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta})}{q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})}\right\}

我们可以将这个期望展开，利用对数的性质将分数分解为两项

\mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) = E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \{\log p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta}) - \log q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})\}\tag{35}

这就是公式(35)。ELBO的名称"证据下界"来自于它与对数似然（证据）的关系，这个关系我们接下来会详细推导。

从对数似然到ELBO的推导

现在我们来建立ELBO与对数似然之间的精确关系。从对数似然的定义出发

\log p(\boldsymbol{y}|\boldsymbol{\theta}) = E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \{\log p(\boldsymbol{y}|\boldsymbol{\theta})\}

这个等式成立是因为对数似然不依赖于 \boldsymbol{z}，所以对 q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) 求期望只是乘以1。接下来我们对 \log p(\boldsymbol{y}|\boldsymbol{\theta}) 进行巧妙的变换。利用条件概率的定义，我们可以写出 p(\boldsymbol{y}|\boldsymbol{\theta}) = \frac{p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta})}{p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta})}，因此

现在我们在分子和分母同时乘以和除以 q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})，这不改变等式的值

\log p(\boldsymbol{y}|\boldsymbol{\theta}) = E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \left\{\log \frac{p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta}) q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})}{q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta})}\right\}

利用对数的性质 \log \frac{ab}{cd} = \log \frac{a}{c} + \log \frac{b}{d}，我们可以将这个对数拆分为两项

\log p(\boldsymbol{y}|\boldsymbol{\theta}) = E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \left\{\log \frac{p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta})}{q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})}\right\} + E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \left\{\log \frac{q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})}{p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta})}\right\}\tag{36}

这就是公式(36)的形式。观察这个等式的结构，我们发现第一项正是ELBO的定义

E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \left\{\log \frac{p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta})}{q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})}\right\} = \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y})

而第二项则是KL散度的定义

E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \left\{\log \frac{q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})}{p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta})}\right\} = \text{KLD}(q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})||p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}))

这个KL散度衡量的是近似后验 q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) 与真实后验 p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}) 之间的差异。

ELBO作为下界的证明

由于KL散度总是非负的，即 \text{KLD}(q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})||p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta})) \geq 0，我们可以从公式(36)直接得出

\mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) = \log p(\boldsymbol{y}|\boldsymbol{\theta}) - \text{KLD}(q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})||p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}))\tag{37}

\mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) \leq \log p(\boldsymbol{y}|\boldsymbol{\theta})

这就是公式(37)，它明确表明ELBO确实是对数似然的一个下界。KL散度项 \text{KLD}(q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})||p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta})) 可以被理解为ELBO与边缘对数似然之间的差距，它量化了这个下界的紧致程度。

最大化ELBO的双重目标

理解了ELBO与对数似然的关系后，我们可以看到最大化ELBO同时达成了两个重要目标。

首先，最大化ELBO会最小化KL散度 \text{KLD}(q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})||p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}))，这使得近似后验 q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) 尽可能接近真实后验 p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta})。这意味着我们的变分推断质量得到了提升，编码器学会了更准确地推断隐变量。

其次，最大化ELBO同时也在最大化边缘似然 p(\boldsymbol{y}|\boldsymbol{\theta})，这正是我们在生成模型中的最终目标。当KL散度趋近于零，也就是 q 完美地近似了 p 时，ELBO就等于对数似然，下界变得紧致。

因此，通过联合优化生成参数 \boldsymbol{\theta} 和变分参数 \boldsymbol{\phi} 来最大化ELBO，我们既改进了推断网络（编码器），也改进了生成模型（解码器），最终得到一个更好的生成模型。

5.3 将EM算法与ELBO联系起来

EM算法作为ELBO的坐标上升

EM算法可以被理解为在ELBO \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) 上进行坐标上升的过程。为了看清这一点，我们从公式(37)出发，考虑在第 m 次迭代时，首先固定 \boldsymbol{\theta}^{(m)} 来最大化ELBO。具体来说，我们要找到最优的 q 分布

q^{(m)}(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = \arg\max_{q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}^{(m)}}(\boldsymbol{y}) = \arg\min_{q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \text{KLD}(q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})||p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}))

这个最大化问题等价于最小化KL散度，因为从公式(37)我们知道

而第一项在固定 \boldsymbol{\theta}^{(m)} 时是常数。KL散度的最小值为零，当且仅当两个分布相等时达到。因此最优的 q 分布就是真实的后验分布

q^{(m)}(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})

这正对应于EM算法的E步：给定当前参数估计，计算隐变量的后验分布。

M步作为对生成参数的优化

接下来的第二步，我们用刚才得到的 q^{(m)}(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 来优化参数 \boldsymbol{\theta}。利用公式(35)中的ELBO表达式，我们要解决的优化问题是

\boldsymbol{\theta}^{(m+1)} = \arg\max_{\boldsymbol{\theta} \in \Theta} E_{\boldsymbol{z} \sim p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})} \{\log p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta}) - \log p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)})\}

注意这个表达式中的第二项 \log p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 不依赖于我们要优化的变量 \boldsymbol{\theta}，所以在优化过程中可以忽略。因此优化问题简化为

\boldsymbol{\theta}^{(m+1)} = \arg\max_{\boldsymbol{\theta} \in \Theta} E_{\boldsymbol{z} \sim p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(t)})} \{\log p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta})\}

这正是标准EM算法M步的形式。

EM与变分推断的本质联系

因此我们可以看到，当我们选择特定的 q^{(m)}(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = p(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\theta}^{(m)}) 时，就得到了公式(15)中给出的标准EM算法，其中完全数据被定义为 \boldsymbol{x} = \{\boldsymbol{y}, \boldsymbol{z}\}。

这揭示了EM算法与变分推断之间的深刻联系：EM可以被看作是一种特殊的变分推断，其中近似后验就是精确后验。而当精确后验不可计算时，我们就需要用参数化的分布族来近似，这就是VAE的思想。

5.4 通过随机梯度下降最大化ELBO

联合优化的必要性

接下来我们讨论如何使用随机梯度下降来优化ELBO，这需要同时对 \boldsymbol{\theta} 和 \boldsymbol{\phi} 进行联合优化。在ELBO不可处理或其梯度难以计算的情况下（常见），这种方法特别有用。基本思路是从 \boldsymbol{\theta} 和 \boldsymbol{\phi} 的随机初始值开始，然后随机地优化它们的值直到收敛。

i.i.d.样本情况下的ELBO分解

在独立同分布样本的情况下，总的ELBO是各个样本ELBO的平均，即

\mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) = \sum_{i=1}^n \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}_i)\tag{38}

这个分解使得我们可以使用小批量梯度下降：每次迭代只使用一个小批量的样本来估计梯度，而不需要遍历整个数据集。单个样本的ELBO及其梯度通常是不可处理的，因此我们需要使用SGD作为 \nabla_{\boldsymbol{\theta}, \boldsymbol{\phi}} \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) 的无偏估计器。

5.4.1 推导对 \boldsymbol{\theta} 的梯度

对于生成参数 \boldsymbol{\theta} 的梯度，我们可以使用简单的蒙特卡罗估计器来获得ELBO的无偏梯度。具体来说，梯度可以写为

\nabla_{\boldsymbol{\theta}} \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) = E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \{\nabla_{\boldsymbol{\theta}}(\log p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta}) - \log q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}))\}

注意到 \log q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) 不依赖于 \boldsymbol{\theta}，所以这一项的梯度为零，我们得到

\nabla_{\boldsymbol{\theta}} \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) = E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \{\nabla_{\boldsymbol{\theta}} \log p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta})\} \approx \nabla_{\boldsymbol{\theta}} \log p(\tilde{\boldsymbol{z}}, \boldsymbol{y}|\boldsymbol{\theta})

其中 \tilde{\boldsymbol{z}} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) 是从变分后验中采样的一个样本。这个单样本蒙特卡罗估计是无偏的，并且我们只需要从 q 中采样，然后计算对数似然对 \boldsymbol{\theta} 的梯度。

5.4.2 通过重参数化技巧推导对 \boldsymbol{\phi} 的梯度

计算关于 \boldsymbol{\phi} 的梯度则更加复杂，因为期望本身依赖于 \boldsymbol{\phi}。具体来说，我们有

\nabla_{\boldsymbol{\phi}} E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \{\log q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})\} \neq E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})} \{\nabla_{\boldsymbol{\phi}} \log q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})\}

问题在于梯度算子不能直接穿过期望算子，因为分布本身依赖于 \boldsymbol{\phi}。但有一个被称为重参数化技巧的巧妙方法可以克服这个问题，它可以被理解为一种变量替换。我们将 \boldsymbol{z} 与另一个随机变量 \boldsymbol{e} 联系起来，通过一个可微且可逆的函数

\boldsymbol{z} = \boldsymbol{r}(\boldsymbol{e}, \boldsymbol{y}, \boldsymbol{\phi})

这里 \boldsymbol{r}(\cdot) 是可微且可逆的，而 \boldsymbol{e} \sim p(\boldsymbol{e}) 是一个不依赖于 \boldsymbol{y} 的随机变量。重参数化的关键是将随机性从依赖 \boldsymbol{\phi} 的 \boldsymbol{z} 转移到不依赖 \boldsymbol{\phi} 的 \boldsymbol{e}，而让 \boldsymbol{\phi} 的依赖性体现在确定性的变换函数 \boldsymbol{r} 中。

有了这个重参数化，我们可以重写关于 \boldsymbol{\phi} 的梯度

\nabla_{\boldsymbol{\phi}} \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) = \nabla_{\boldsymbol{\phi}} E_{\boldsymbol{e} \sim p(\boldsymbol{e})} \{\log p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta}) - \log q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi})\}

现在期望是关于不依赖 \boldsymbol{\phi} 的分布 p(\boldsymbol{e})，所以梯度可以移到期望里面

\nabla_{\boldsymbol{\phi}} \mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}) \approx \nabla_{\boldsymbol{\phi}} \log p(\tilde{\boldsymbol{z}}, \boldsymbol{y}|\boldsymbol{\theta}) - \log q(\tilde{\boldsymbol{z}}|\boldsymbol{y}, \boldsymbol{\phi})

其中 \tilde{\boldsymbol{z}} = \boldsymbol{r}(\tilde{\boldsymbol{e}}, \boldsymbol{y}, \boldsymbol{\phi})，然后 \tilde{\boldsymbol{e}} \sim p(\boldsymbol{e})。这样我们就得到了一个可以通过采样和自动微分计算的梯度估计。

例5.1：高斯推断模型的ELBO无偏估计计算

模型设定

接下来我们通过一个具体例子来说明如何计算ELBO的无偏估计。我们考虑以下设定。

首先是高斯推断模型。我们假设

q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = \mathcal{N}(\bar{\boldsymbol{\mu}}(\boldsymbol{\phi}, \boldsymbol{x}), \bar{\boldsymbol{\Sigma}}(\boldsymbol{\phi}, \boldsymbol{x}))

这个分布的目标是近似真实的后验 p(\boldsymbol{z}|\boldsymbol{y})。其中均值和协方差矩阵将通过一个参数为 \boldsymbol{\phi} 的编码器神经网络来学习，输入是观测 \boldsymbol{y}，记为 \text{ENN}(\boldsymbol{\phi}, \boldsymbol{y})。

其次是生成模型。生成过程定义为

p(\boldsymbol{z}, \boldsymbol{y}|\boldsymbol{\theta}) = p(\boldsymbol{y}|\boldsymbol{z}, \boldsymbol{\theta})p(\boldsymbol{z})

其中条件似然 p(\boldsymbol{y}|\boldsymbol{z}, \boldsymbol{\theta}) = \mathcal{N}(\boldsymbol{\mu}(\boldsymbol{\theta}, \boldsymbol{z}), \boldsymbol{\Sigma}(\boldsymbol{\theta}, \boldsymbol{y}))。

这里均值和协方差也将通过一个参数为 \boldsymbol{\theta} 的解码器神经网络来学习，输入是隐变量 \boldsymbol{z}，记为 \text{DNN}(\boldsymbol{\phi}, \boldsymbol{y})。先验分布我们选择标准正态分布 p(\boldsymbol{z}) = \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})。

重参数化的具体形式

从高斯推断模型 q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = \mathcal{N}(\bar{\boldsymbol{\mu}}(\boldsymbol{\phi}, \boldsymbol{y}), \bar{\boldsymbol{\Sigma}}(\boldsymbol{\phi}, \boldsymbol{y})) 出发，我们可以写出重参数化的形式。引入一个辅助噪声变量

\boldsymbol{e} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})\tag{39}

然后通过线性变换生成 \boldsymbol{z}

\boldsymbol{z} = \bar{\boldsymbol{\mu}}(\boldsymbol{\phi}, \boldsymbol{y}) + \bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y})\boldsymbol{e}\tag{39}

这个参数化对应于协方差矩阵的Cholesky分解，即 \bar{\boldsymbol{\Sigma}}(\boldsymbol{\phi}, \boldsymbol{y}) = \bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y})\bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y})^\top。为了对协方差矩阵强制使用Cholesky形式，我们可以让编码器网络输出

\bar{\boldsymbol{\mu}}(\boldsymbol{\phi}, \boldsymbol{y}), \bar{\boldsymbol{\sigma}}(\boldsymbol{\phi}, \boldsymbol{y}), \bar{\boldsymbol{L}}'(\boldsymbol{\phi}, \boldsymbol{y}) \leftarrow \text{ENN}(\boldsymbol{\phi}, \boldsymbol{y})

然后构造Cholesky因子

\bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y}) \leftarrow \boldsymbol{L}_{\text{mask}} \odot \bar{\boldsymbol{L}}'(\boldsymbol{\phi}, \boldsymbol{y}) + \text{diag}(\bar{\boldsymbol{\sigma}}(\boldsymbol{\phi}, \boldsymbol{y}))

这里 \boldsymbol{L}_{\text{mask}} 是一个在对角线及以上元素为零的掩码矩阵，\odot 表示逐元素乘法。这样构造保证了 \bar{\boldsymbol{L}} 是下三角矩阵，满足Cholesky分解的要求。有了这个参数化，我们就可以从公式(39)中采样 \boldsymbol{z}。

ELBO各项的计算

从公式(35)和(38)出发，单个样本的ELBO可以写为

\mathcal{E}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}_i) = E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi})} \{\log p(\boldsymbol{y}_i|\boldsymbol{z}, \boldsymbol{\theta})\} + \log p(\boldsymbol{z}) - E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi})} \{\log q(\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi})\}\tag{40}

这就是公式(40)。ELBO包含三项：重构项 E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi})} \{\log p(\boldsymbol{y}_i|\boldsymbol{z}, \boldsymbol{\theta})\}，先验项 \log p(\boldsymbol{z})，以及熵项 -E_{\boldsymbol{z} \sim q(\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi})} \{\log q(\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi})\}。

公式(40)中的ELBO包含需要计算的项 q(\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi})。由于我们使用了重参数化，\boldsymbol{e} 和 \boldsymbol{z} 之间通过可逆函数 \boldsymbol{r}(\cdot) 联系，根据随机变量变换的公式，两个密度之间的关系为

q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = \log p(\boldsymbol{e}) - \log \left|\det \frac{\partial \boldsymbol{z}}{\partial \boldsymbol{e}}\right|

由于 \boldsymbol{z} = \bar{\boldsymbol{\mu}}(\boldsymbol{\phi}, \boldsymbol{y}) + \bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y})\boldsymbol{e}，Jacobian矩阵就是 \bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y})，其行列式等于对角元素的乘积。因此

q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = \log \mathcal{N}(\boldsymbol{0}, \boldsymbol{I}) - \sum_j \log |[\bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y})]_{jj}|

由于标准正态分布的对数密度和下三角矩阵对角元素的对数，我们得到

q(\boldsymbol{z}|\boldsymbol{y}, \boldsymbol{\phi}) = -\frac{1}{2} \sum_j \left([\boldsymbol{e}]_j^2 + \log(2\pi) + \log([\bar{\boldsymbol{\sigma}}]_j)\right)\tag{41}

这就是公式(41)。

完整算法流程

将所有这些步骤总结起来，我们得到了计算高斯推断模型下ELBO无偏估计的完整算法。算法首先初始化生成参数 \boldsymbol{\theta} 和变分参数 \boldsymbol{\phi}，然后对每个样本 \boldsymbol{y}_i 执行以下步骤。

第一步是通过编码器网络获得变分后验的参数

\bar{\boldsymbol{\mu}}(\boldsymbol{\phi}, \boldsymbol{y}_i), \bar{\boldsymbol{\sigma}}(\boldsymbol{\phi}, \boldsymbol{y}_i), \bar{\boldsymbol{L}}'(\boldsymbol{\phi}, \boldsymbol{y}_i) \leftarrow \text{ENN}(\boldsymbol{\phi}, \boldsymbol{y}_i)

编码器网络输出三个量：后验分布的均值向量、对角方差元素，以及Cholesky因子的初步形式。

接下来我们需要构造完整的Cholesky因子。这个步骤确保协方差矩阵具有正确的下三角结构

\bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y}_i) \leftarrow \boldsymbol{L}_{\text{mask}} \odot \bar{\boldsymbol{L}}'(\boldsymbol{\phi}, \boldsymbol{y}_i) + \text{diag}(\bar{\boldsymbol{\sigma}}(\boldsymbol{\phi}, \boldsymbol{y}_i))

这里 \boldsymbol{L}_{\text{mask}} 是掩码矩阵，\odot 表示逐元素乘法，\text{diag}(\cdot) 构造对角矩阵。

有了后验分布的参数后，我们开始采样过程。首先从标准正态分布采样噪声变量

\boldsymbol{e} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})

然后通过重参数化技巧生成隐变量样本

\boldsymbol{z} \leftarrow \bar{\boldsymbol{\mu}}(\boldsymbol{\phi}, \boldsymbol{y}_i) + \bar{\boldsymbol{L}}(\boldsymbol{\phi}, \boldsymbol{y}_i)\boldsymbol{e}

这个变换将不依赖参数的随机性（通过 \boldsymbol{e}）与依赖参数的确定性变换（通过均值和Cholesky因子）分离开来，使得梯度计算成为可能。

得到隐变量样本 \boldsymbol{z} 后，我们将其输入解码器网络以获得条件似然的参数

\boldsymbol{\mu}(\boldsymbol{\theta}, \boldsymbol{z}_i), \boldsymbol{\Sigma}(\boldsymbol{\theta}, \boldsymbol{y}_i) \leftarrow \text{DNN}(\boldsymbol{\theta}, \boldsymbol{z})

解码器网络输出生成分布 p(\boldsymbol{y}|\boldsymbol{z}, \boldsymbol{\theta}) 的均值和协方差参数。

现在我们可以计算ELBO的三个组成部分。首先是先验项的对数概率

\tilde{\mathcal{E}}_{\boldsymbol{z}} = \log \mathcal{N}(\boldsymbol{z}; \boldsymbol{0}, \boldsymbol{I})

这一项衡量采样得到的隐变量在先验分布下的概率。

第二项是变分后验熵的近似，通过重参数化我们可以将其表示为

\tilde{\mathcal{E}}_{\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi}} = -\frac{1}{2} \sum_j \left([\boldsymbol{e}]_j^2 + \log(2\pi) + \log([\bar{\boldsymbol{\sigma}}]_j)\right)

这一项通过噪声变量 \boldsymbol{e} 和编码器输出的方差参数来计算变分分布的对数概率。

第三项是条件似然项的近似

\tilde{\mathcal{E}}_{\boldsymbol{y}_i|\boldsymbol{z}, \boldsymbol{\theta}} = \log \mathcal{N}(\boldsymbol{z}; \boldsymbol{\mu}(\boldsymbol{\theta}, \boldsymbol{z}_i), \boldsymbol{\Sigma}(\boldsymbol{\theta}, \boldsymbol{y}_i))

这一项衡量在给定隐变量 \boldsymbol{z} 的情况下，观测数据 \boldsymbol{y}_i 的重构质量。

最后，我们将这三项组合起来得到单个样本ELBO的无偏蒙特卡罗估计

\tilde{\mathcal{E}}_{\boldsymbol{\phi}, \boldsymbol{\theta}}(\boldsymbol{y}_i) = \tilde{\mathcal{E}}_{\boldsymbol{y}_i|\boldsymbol{z}, \boldsymbol{\theta}} + \tilde{\mathcal{E}}_{\boldsymbol{z}} - \tilde{\mathcal{E}}_{\boldsymbol{z}|\boldsymbol{y}_i, \boldsymbol{\phi}}

这个组合遵循ELBO的定义：重构项加先验项减去变分后验的熵项。

算法的实际应用

这个算法给出了ELBO的一个无偏的单样本蒙特卡罗估计，可以直接用于随机梯度下降优化。通过现代深度学习框架提供的自动微分工具，我们可以高效地计算这个估计对生成参数 \boldsymbol{\theta} 和变分参数 \boldsymbol{\phi} 的梯度。得到梯度后，我们就可以使用标准的梯度下降或其变体（如Adam优化器）来更新参数。这个完整的流程——从数据到编码器、通过重参数化采样、到解码器、计算ELBO、反向传播梯度、更新参数——构成了变分自编码器训练的核心算法。通过在整个数据集上迭代这个过程，编码器学会了如何有效地推断隐变量，而解码器学会了如何从隐变量生成逼真的数据，两个网络协同工作最终得到一个强大的生成模型。

Part II : 贝叶斯计算与采样模拟方法

这部分内容将涵盖贝叶斯估计的理论回顾、各类模拟算法、采样方法以及实际应用。

首先回顾贝叶斯估计的基础框架和数值积分方法。贝叶斯方法的核心在于通过贝叶斯定理将先验信息与观测数据结合，得到参数的后验分布。然而，由于后验分布往往没有解析形式，需要通过数值积分或采样方法来近似计算后验均值、方差等统计量。这部分会讨论为什么需要这些数值方法，以及它们在贝叶斯推断中扮演的角色。

然后针对采样部分，我们希望能够从分布中直接采样，但实际应用中，很多分布并不能直接采样，因此我们介绍几类主要模拟抽样方法:

累积分布函数（CDF）反演方法基于概率积分变换。其核心是采用概率分布函数的逆变换来将均分分布的采样变成目标分布的采样，其缺点在于逆函数可能难以计算

重要性采样方法不直接从目标分布采样，而是从一个容易采样的提议分布采样，然后通过重要性权重来修正。在数学理论上，我们引入一个新的分布（提议分布），对其进行积分期望的变换，对新分布进行新函数的计算。

接受-拒绝采样方法则是从提议分布采样后，根据接受概率决定是否保留该样本，最终得到的样本严格服从目标分布。

最后我们会介绍Gibbs采样器与Metropolis-Hastings算法，其本质是马尔可夫链蒙特卡洛（MCMC）方法，适用于目标分布是多维的情况（且难以采样）。

Gibbs采样器适用于虽然联合分布复杂，但条件分布相对简单的情况，其核心是固定其他维度的值，只从当前维度的条件分布中采样，然后不断循环这个过程，生成的样本会形成一个马尔可夫链，其平稳分布正是我们想要的联合分布。该方法的优势在于如果条件分布是常见的标准分布，我们可以直接采样而不需要接受-拒绝步骤，因此效率较高

Metropolis-Hastings算法是更通用的MCMC方法。它通过构造一个转移核，从当前状态提议新状态，然后根据接受概率决定是否转移到新状态。接受概率的设计保证了马尔可夫链的平稳分布就是目标分布。其优势在于它不要求条件分布易于采样，只需要能够计算目标分布的值（甚至不需要归一化常数）。

除此之外，在学习这些算法的过程中我们还会回顾蒙特卡洛、马尔可夫链等细节内容，比如随机采样等无偏性，MC的遍历性等，并使用随机游走作为相关示例

逆问题的贝叶斯方法回顾

问题的基本构成

一个典型的逆问题通常包含三个基本要素。

第一个要素是待求的未知量，这通常是我们最感兴趣的目标，例如一张清晰的图像，我们用变量 x 来表示。

第二个要素是问题的物理模型，它描述了未知量 x 是如何产生观测数据 y 的，这个模型本身可能包含一些已知的参数。

第三个要素是先验知识，也就是关于未知量 x 本身性质的预先了解。比如，我们可能知道自然图像通常是平滑的，而不是完全随机的噪声。先验知识往往通过正则项来体现，它能帮助我们从众多可能的解中找到一个更符合实际的解。

正向模型的数学描述

我们首先从一个基础的线性模型来描述观测过程。这个模型是正向模型，它描述了从未知信号 x 生成观测数据 y 的完整过程。其核心结构可以表示为一个流程:未知的原始信号 x 经过观测算子 H 的作用，然后与噪声 b 相加，最终得到我们能够直接观测到的数据 y。

在这个流程中，我们首先有一个未知的原始信号 x，例如一张高分辨率的清晰图像。这个信号经过系统 H （观测算子）的作用，它将原始信号从一个高维空间 \mathbb{R}^N 映射到一个观测空间 \mathbb{R}^M，它可以是线性的，也可以是非线性的。在这个过程中，不可避免地会引入噪声，我们用变量 b 来表示。噪声与系统作用后的信号 Hx 相加，最终得到了我们能够直接观测到的数据 y。

因此，整个正向模型的数学方程为:

y = Hx + b

y \in \mathbb{R}^M 是我们得到的观测数据
x \in \mathbb{R}^N 是我们希望恢复的未知信号
H 是描述物理过程的观测算子
b \in \mathbb{R}^M 是加性噪声

所谓直接问题或正向问题，指的是在已知 x、H 和 b 的情况下，去计算 y 的过程，也就是模拟观测数据的生成。而我们的核心任务——逆问题，则是在已知观测数据 y 和算子 H 的情况下，反向推断出最有可能的原始信号 x。

统计推断的基本思想

为了从数学上解决这个逆问题，我们引入统计推断的思想。核心目标是根据已有的观测数据 y，构建一个有效的估计器，用它来计算出对未知信号 x 的一个估计值，我们记作 \hat{x}。这个估计器本质上是一个函数，我们将它记为 f，它以观测数据 y 作为输入，输出一个估计值 \hat{x}:

\hat{x} = f(y)

这个过程可以看作是一个从观测空间 \mathbb{R}^M 到解空间 \mathbb{R}^N 的映射。一个好的估计器 f 是基于某个最优化准则来构建的，比如说最小化某个误差（估计值 \hat{x} 与真实值 x 之间的差距），或者是最大化某个概率（在给定观测 y 的条件下某个 x 出现的后验概率）。

这种最优化问题的求解，有时可以直接得到一个明确的解析解，但更多时候需要依赖于迭代的算法来逐步逼近最优解。

由于噪声 b 的存在以及观测过程可能造成的信息损失，例如 M < N，我们通过估计器得到的解 \hat{x} 无法与真实信号 x （永远无法直接观测到）完全一致，即 \hat{x} \neq x。我们只能让 \hat{x} 尽可能地接近 x。

概率建模的引入

为了处理问题中的不确定性，我们进行概率建模。逆问题中的不确定性主要来源于噪声 b。因此想法是，不再将 b 视为一个固定的未知数，而是将其建模为一个随机变量。我们可以假设 b 服从某个特定的概率分布 p(b)。一个最常见且应用广泛的假设是，噪声 b 服从高斯分布，其均值为 \mu，协方差矩阵为 \Sigma:

b \sim \mathcal{N}(\mu, \Sigma)

一旦我们为噪声 b 建立了概率模型，我们就可以推导出似然函数，它描述了在给定一个具体的 x 的条件下，观测到数据 y 的概率 p(y|x)。根据我们的模型 y = Hx + b，如果 x 是一个已知的、确定的值，那么 y 就是随机变量 b 经过一个线性变换后得到的新随机变量，且如果 b \sim \mathcal{N}(\mu， \Sigma)，那么 y 在 x 给定的条件下，也服从一个高斯分布:

p(y|x) \sim \mathcal{N}(Hx + \mu, \Sigma)

也就是说，对噪声 b 分布的假设 p(b)，直接决定了似然函数 p(y|x) 的具体形式，因此贝叶斯思想的关键一步就是将噪声 b 和观测 y 都看作是随机变量。在更完整的贝叶斯框架下，我们甚至将待求解的 x 也视为一个随机变量。

概率密度函数与似然函数的区别

p(y|x) 和 L(x|y) 在数值上是相等的，但物理意义不同。

p(y|x) 是一个关于 y 的概率密度函数。对于一个固定的 x，它描述了 y 的概率分布，因此它对 y 的积分等于1，即 \int p(y|x) dy = 1。
L(x|y) 被看作是关于未知参数 x 的函数。在观测到具体的数据 y 之后，我们通过改变 x 的值来考出现这个观测数据 y 的可能性有多大。

最大似然估计的实例

假设我们有一组观测数据 (\gamma_1， \gamma_2， \dots， \gamma_n)，它们是独立同分布的，并且我们知道它们都来自于一个伽马分布，但这个分布的参数是未知的。伽马分布的概率密度函数形式为:

p(\gamma | \alpha, \beta) = \frac{1}{\Gamma(\alpha)\beta^\alpha} \gamma^{\alpha-1} e^{-\frac{\gamma}{\beta}}, \quad \text{其中 } \alpha>0, \beta>0

由于所有观测是独立同分布的，那么观测到整个数据集 (\gamma_1， \dots， \gamma_n) 的联合概率，也就是似然函数 L(\alpha， \beta | \gamma_1， \dots， \gamma_n)，就是每个观测点概率的乘积:

L(\alpha, \beta | \gamma_1, \dots, \gamma_n) = p(\gamma_1, \dots, \gamma_n | \alpha, \beta) = \prod_{i=1}^{n} p(\gamma_i | \alpha, \beta)

我们应该选择一组参数 (\hat{\alpha}， \hat{\beta})，使得我们当前观测到的这组数据出现的可能性（似然函数的值），达到最大。

(\hat{\alpha}, \hat{\beta}) = \underset{\alpha, \beta}{\text{argmax}} \ L(\alpha, \beta | \gamma_1, \dots, \gamma_n)

数取对数，将连乘变为连加，因此对数似然函数为:

(\hat{\alpha}, \hat{\beta}) = \underset{\alpha, \beta}{\text{argmax}} \ \log L(\alpha, \beta | \gamma_1, \dots, \gamma_n)

我们需要对它求导并令导数为零。对于参数 \beta，这个方程通常有显式解，也就是一个封闭的数学表达式。但对于参数 \alpha，往往没有显式解，必须依赖数值优化算法来迭代求解。

不适定问题

一个数学问题如果被称为适定的，它必须同时满足以下三个条件:解必须存在，解必须是唯一的，解必须是稳定的。

稳定性意味着当输入数据，即我们的观测值 y，发生微小扰动时，输出的解 x 也只发生微小的变化。

如果以上任何一个条件不被满足，该问题就被称为不适定的，绝大多数有实际意义的逆问题都是不适定的。

我们还是通过线性模型 y = Hx + b 来具体说明。假设噪声 b 服从均值为0、协方差矩阵为单位矩阵 I 的高斯分布，即 b \sim \mathcal{N}(0， I)。在这种情况下，我们之前推导的似然函数 p(y|x) 是一个高斯分布。其对数形式为:

\log L(x|y) \propto -\frac{1}{2} \|y - Hx\|^2

最大化这个对数似然函数，等价于最小化欧几里得距离的平方 \|y - Hx\|^2。这个问题被称为最小二乘问题。它的解 \hat{x} 满足:

\hat{x} = \underset{x}{\text{argmin}} \ \|y - Hx\|^2

通过对目标函数关于 x 求梯度并令其为零，我们可以得到所谓的正规方程:

(H^T H) \hat{x} = H^T y

如果矩阵 H^T H 是可逆的，我们可以直接得到解:

\hat{x} = (H^T H)^{-1} H^T y

但是这里面有很大的问题。首先，如果 H^T H 不可逆，例如当观测数据的维度 M 小于未知信号的维度 N 时，解就不是唯一的，问题是不适定的。更严重的是，即使 H^T H 在理论上是可逆的，它也可能具有非常小的特征值，这会导致其逆矩阵 (H^T H)^{-1} 具有非常大的特征值。

我们可以将 y = Hx + b 代入解的表达式中来看看这对噪声的影响:

\hat{x} = (H^T H)^{-1} H^T (Hx + b) = x + (H^T H)^{-1} H^T b

我们的估计值 \hat{x} 等于真实信号 x 加上一个被严重放大的噪声项 (H^T H)^{-1} H^T b。这种噪声放大的现象是造成解不稳定性的根本原因 —— 微小的观测噪声 b 会被矩阵 (H^T H)^{-1} 放大到足以完全淹没真实信号 x 的程度。

为了克服这个问题，我们需要引入正则化。其核心思想是在最小二乘的目标函数中加入一个惩罚项，这个惩罚项代表了我们对解 x 的先验知识。一个最常见的正则化方法是Tikhonov正则化:

\hat{x} = \underset{x}{\text{argmin}} \ \left( \|y - Hx\|^2 + \lambda \|x\|^2 \right)

这里的 \lambda \|x\|^2 就是正则项，它会惩罚那些范数过大的解。参数 \lambda > 0 是正则化参数，它控制着数据拟合项和正则项之间的平衡，这个正则化后的问题有稳定的解析解:

\hat{x} = (H^T H + \lambda I)^{-1} H^T y

通过在 H^T H 的对角线上加上一个正数 \lambda，我们保证了矩阵 (H^T H + \lambda I) 始终是可逆且良态的，从而有效地抑制了噪声的爆炸式增长。

贝叶斯方法的核心思想

正则化是一种有效但带有经验性的方法，而贝叶斯方法的核心是:将问题中所有的未知量都建模为随机变量，并用概率分布来描述它们的不确定性。在贝叶斯框架下，我们不仅认为观测 y 和噪声 b 是随机的，我们连希望求解的未知信号 x 本身也看作是一个随机变量。我们最关心的量是后验概率分布 p(x|y)。这个分布描述了在观测到数据 y 之后，我们对于未知量 x 的全部认知。

这个后验分布可以通过贝叶斯定理得到。根据概率的乘法法则，联合概率 p(x,y) 有两种分解方式:

p(x,y) = p(y|x)p(x) = p(x|y)p(y)

对上式进行简单的移项，我们就得到了贝叶斯定理的核心形式:

p(x|y) = \frac{p(y|x)p(x)}{p(y)}

p(x|y)：后验概率 (Posterior)。这是我们推断的目标，代表了结合数据和先验知识后对 x 的认识。
p(y|x)：似然函数 (Likelihood)。它与我们之前讨论的似然函数完全一样，描述了在给定 x 的情况下，观测到 y 的概率。
p(x)：先验概率 (Prior)。它代表了在进行任何观测之前，我们对 x 已有的知识或经验。例如，我们可以用一个高斯先验 p(x) = \mathcal{N}(0, C) 来表明我们相信 x 的分量不会太大，本质上与正则化中的惩罚项是一样的。
p(y)：证据 (Evidence) 或边际似然。它等于联合概率对 x 的积分：p(y) = \int p(y|x)p(x) dx。对于一个给定的观测 y，它的值是固定的，因此它起到归一化常数的作用，确保后验概率 p(x|y) 对 x 积分后等于1。

在很多情况下，我们只关心哪个 x 使得后验概率最大，而不需要计算后验概率的具体数值。在这种情况下，我们可以忽略分母 p(y)，得到一个更简洁的比例关系:

p(x|y) \propto p(y|x)p(x)

这个关系式阐述了贝叶斯推断的本质:后验知识是似然和先验的结合。一个特别理想的情况是，当我们为似然 p(y|x) 和先验 p(x) 都选择高斯分布时，例如 p(y|x) \sim \mathcal{N}(Hx， \Sigma) 和 p(x) \sim \mathcal{N}(0， C)，它们相乘得到的后验分布 p(x|y) 也会是一个高斯分布。

贝叶斯决策理论

贝叶斯推断的结果是完整的后验分布 p(x|y)，它包含了关于 x 的所有不确定性信息。但我们需要给出一个确定的估计值 \hat{x}，这就需要贝叶斯决策过程。

为了从一个概率分布中选择一个最优的点，我们需要一个评判标准——损失函数 L(x， \hat{x})。它量化了当真实值为 x 而我们的估计值为 \hat{x} 时所付出的代价或损失。损失函数一个非负函数，即 L(x， \hat{x}) \in \mathbb{R}^+。

我们要找到的估计器 \hat{x} = f(y) 必须在所有可能情况下的平均损失最小。这个平均损失被称为贝叶斯风险 RB，其定义为损失函数在 x 和 y 的联合分布下的期望值:

RB = \iint L(x, f(y)) p(x, y) dy dx = \iint L(x, f(y)) p(y|x)p(x) dy dx= \iint L(x, f(y)) p(x|y)p(y) dy dx

直接最小化这个关于函数 f 的双重积分很复杂，我们可以把这个双重积分改写成嵌套积分的形式。利用最后一个等式中的p(x|y)p(y)，我们可以写成：

RB = \int \left[\int L(x, f(y))p(x|y)dx\right]p(y)dy

外层积分是对y积分，内层积分是对x积分。先看内层积分，表示当观测值固定到某个 y 时，损失函数对所有可能的真实值 x 的后验分布期望，外层积分则是把所有可能的观测值 y 再求一次期望，因此，如果我们想要最小化整个贝叶斯风险RB，也就是选择最优的估计器函数f(y)，我们实际上可以逐点地进行优化。

具体来说，对于每一个可能的观测值 y，我们只需要最小化内层的积分：

\hat{x}(y) = \underset{\hat{x}}{\text{argmin}} \int L(x, \hat{x}) p(x|y) dx

这个积分是对后验分布 p(x|y) 求期望，代表如果我们对每个 y 都找到了使内层积分最小的 \hat{x}(y)，那么外层积分（对 y 的加权平均）自然也会达到最小值，即优化问题从对所有可能的 (x,y) 优化简化为对每个固定的 y，对所有可能的 x 优化。

因此后续的目标我们只需要看内层积分即可，现在我们选择一个非常常用且直观的损失函数——二次损失:

L(x, \hat{x}) = \|x - \hat{x}\|^2=\|x - f(y)\|^2

在这种损失函数下，可以证明，能使后验期望损失最小的那个估计值 \hat{x}，恰好就是后验分布 p(x|y) 的均值:

\hat{x}_{MMSE} =\underset{\hat{x}}{\text{argmin}} \ RB = \int x p(x|y) dx = \mathbb{E}_{p(x|y)}[x]

下面给出证明

给定二次损失函数 L(x， \hat{x}) = \|x - \hat{x}\|^2，我们要最小化后验期望损失:

RB(\hat{x}) = \int L(x, \hat{x})p(x|y)dx = \int \|x - \hat{x}\|^2 p(x|y)dx

第一步是展开二次项:

\|x - \hat{x}\|^2 = x^Tx - 2x^T\hat{x} + \hat{x}^T\hat{x}

第二步是代入积分:

RB(\hat{x}) = \int x^Tx \, p(x|y)dx - 2\hat{x}^T\int x \, p(x|y)dx + \hat{x}^T\hat{x}\int p(x|y)dx

由于 \int p(x|y)dx = 1，即概率密度积分为1，简化为:

RB(\hat{x}) = \int x^Tx \, p(x|y)dx - 2\hat{x}^T\int x \, p(x|y)dx + \hat{x}^T\hat{x}

第三步是对 \hat{x} 求导并令其为零:

\frac{\partial RB}{\partial \hat{x}} = -2\int x \, p(x|y)dx + 2\hat{x} = 0

第四步解得最优估计值:

\hat{x}_{MMSE} = \int x \, p(x|y)dx = \mathbb{E}_{p(x|y)}[x]

这正是后验分布 p(x|y) 的均值。为了验证这确实是最小值点，我们对 \hat{x} 求二阶导数:

\frac{\partial^2 RB}{\partial \hat{x}^2} = 2I > 0

其中 I 是单位矩阵，所以这确实是最小值点。这个估计器被称为最小均方误差估计器。它告诉我们，在二次损失的准则下，最好的点估计就是对后验分布求均值。这是贝叶斯决策理论中一个极其重要的结果。

0-1损失函数

然而这并非唯一的选择，不同的损失函数反映了对误差的不同考量，从而引出不同的最优估计策略。例如，我们可以采用0-1损失函数，该准则规定，只有当估计值 \hat{x}(y) 与真实值 x 完全一致时，损失才为零，否则损失就是一个固定的正值，通常设为1。其数学形式为:

L(x, \hat{x}(y)) = \begin{cases} 0 & \text{if } x^* = f(y) \\ 1 & \text{if } x^* \neq f(y) \end{cases}

在这种要么对要么错的评价体系下，为了最小化平均损失，我们的最佳策略是选择那个最有可能成为真实值的点。这恰好对应了后验概率分布 p(x|y) 中概率密度最大的点，即众数。通过这种方式得到的估计器被称为最大后验概率估计器:

\hat{x}_{MAP}(y) = \underset{x}{\text{argmax}} \ p(x|y)

与需要计算积分的后验均值相比，MAP估计将问题转化为一个优化问题，在实际计算中通常更为便捷。

绝对值损失函数

另一种常见的选择是绝对值损失函数，它用估计误差的绝对值来衡量损失:

L(x, \hat{x}(y)) = |x - \hat{x}(y)|

当采用此损失函数时，可以证明，能够最小化后验期望损失的估计值是后验分布 p(x|y) 的中位数。

总结来说，三种最核心的贝叶斯估计器分别对应了三种不同的损失准则:均方误差损失导向后验均值，0-1损失导向后验众数，绝对值损失则导向后验中位数。

积分问题

无论是为了获得MMSE估计而去计算后验均值，还是为了评估估计结果的不确定性而去计算后验方差，我们都面临着一个共同的计算瓶颈：求解积分。在大多数有实际价值的问题中，后验分布 p(x|y) 的形式都相当复杂，并且未知量 x 的维度可能非常高，这使得上述积分无法通过解析方法求得精确解。因此我们必须转向能够近似计算这些积分的数值方法，其中最强大和最通用的工具就是蒙特卡洛方法。

蒙特卡洛方法

蒙特卡洛方法是一大类依赖于随机抽样来获得数值结果的计算算法，适合解决高维、复杂的积分问题。

假设我们的目标是计算积分 I = \int h(x) p(x) dx，这等价于计算函数 h(x) 在概率分布 p(x) 下的期望值 \mathbb{E}_p[h(x)]。

蒙特卡洛方法的基本流程非常直观：首先我们设法从目标分布 p(x) 中抽取大量独立的样本，记为 \{x_1， x_2， \dots， x_N\}，然后将这些样本代入函数 h(x) 中得到一系列的函数值 \{h(x_1)， h(x_2)， \dots， h(x_N)\}，最后计算这些函数值的算术平均值，以此作为对真实积分值 I 的估计，记作 \hat{I}_N：

\hat{I}_N = \frac{1}{N} \sum_{i=1}^{N} h(x_i)

根据概率论中的大数定律，只要样本数量 N 足够大，这个由样本均值构成的估计量 \hat{I}_N 就会依概率收敛到真实的期望值 I。

蒙特卡洛估计的统计性质

一个有效的估计方法需要具备优良的统计性质。蒙特卡洛估计器 \hat{I}_N 的表现可以通过其无偏性、方差和渐近分布来刻画。

首先该估计器是无偏的，这意味着从期望意义上讲我们的估计是准确的。我们可以通过计算其期望值来验证这一点。根据期望的线性性质：

\mathbb{E}[\hat{I}_N] = \mathbb{E}\left[\frac{1}{N} \sum_{i=1}^{N} h(x_i)\right] = \frac{1}{N} \sum_{i=1}^{N} \mathbb{E}[h(x_i)]

由于每个样本 x_i 都来自于同一个分布 p(x)，所以 \mathbb{E}[h(x_i)] 的值都等于我们想要估计的积分 I。因此：

\mathbb{E}[\hat{I}_N] = \frac{1}{N} (NI) = I

其次我们可以通过分析其方差来评估估计的稳定性。估计量的方差越小，单次估计的结果就越可能接近真实值。由于我们采集的样本是独立同分布的，所以 h(x_i) 这些随机变量也是相互独立的。对于独立随机变量，和的方差等于方差的和，因此：

\text{Var}(\hat{I}_N) = \text{Var}\left(\frac{1}{N} \sum_{i=1}^{N} h(x_i)\right) = \frac{1}{N^2} \sum_{i=1}^{N} \text{Var}(h(x_i))

最终我们得到一个非常关键的结果：

\text{Var}(\hat{I}_N) = \frac{\text{Var}(h(x))}{N}

这个公式表明，蒙特卡洛估计的方差与单个样本的方差 \text{Var}(h(x)) 成正比，与样本数量 N 成反比。

最后根据中心极限定理，当样本数量 N 足够大时，蒙特卡洛估计器 \hat{I}_N 的概率分布会趋向于一个正态分布：

\hat{I}_N \sim \mathcal{N}\left(I, \text{Var}(\hat{I}_N)\right)

这个渐近正态性使我们能够为估计结果构建置信区间，从而对估计的不确定性进行量化。

方差缩减技术

从估计量的方差公式可以看出：

\text{Var}(\hat{I}_N) = \frac{\text{Var}(h(x))}{N}

蒙特卡洛方法的收敛效率直接受制于被积函数 h(x) 本身的方差。如果这个方差过大，我们就需要海量的样本才能达到满意的精度。因此研究如何降低估计量方差的方差缩减技术，是蒙特卡洛方法中的一个核心课题。

Rao-Blackwell定理为我们提供了强大的方差缩减思路，其理论基础是全方差公式。对于任意两个随机变量 Y 和 Z，变量 Y 的总方差可以被分解为两部分：

\text{Var}(Y) = \mathbb{E}[\text{Var}(Y|Z)] + \text{Var}(\mathbb{E}[Y|Z])

由于方差项永远是非负的，即 \mathbb{E}[\text{Var}(Y|Z)] \ge 0，我们可以得到不等式：

\text{Var}(Y) \ge \text{Var}(\mathbb{E}[Y|Z])

我们先详细证明全方差公式:

\text{Var}(Y) = \mathbb{E}[\text{Var}(Y|Z)] + \text{Var}(\mathbb{E}[Y|Z])

步骤一

在步骤一中，我们首先从条件方差的定义出发。对于给定的 Z，Y 的条件方差被定义为：

\text{Var}(Y|Z) = \mathbb{E}[Y^2|Z] - (\mathbb{E}[Y|Z])^2

现在对这个等式的两边同时取关于 Z 的期望：

\mathbb{E}[\text{Var}(Y|Z)] = \mathbb{E}\left[ \mathbb{E}[Y^2|Z] - (\mathbb{E}[Y|Z])^2 \right]

利用期望的线性性质，我们可以将上式拆分为两项：

\mathbb{E}[\text{Var}(Y|Z)] = \mathbb{E}[\mathbb{E}[Y^2|Z]] - \mathbb{E}[(\mathbb{E}[Y|Z])^2]

这里我们需要引入迭代期望定律，也常被称为塔式法则。该定律指出，对一个条件期望再次求期望会得到原始的无条件期望：

\mathbb{E}[\mathbb{E}[X|Z]] = \mathbb{E}[X]

将此定律应用于上式的第一项：

\mathbb{E}[\mathbb{E}[Y^2|Z]] = \mathbb{E}[Y^2]

因此步骤一的最终表达式为：

\mathbb{E}[\text{Var}(Y|Z)] = \mathbb{E}[Y^2] - \mathbb{E}[(\mathbb{E}[Y|Z])^2]

步骤二

在步骤二中，我们将 \mathbb{E}[Y|Z] 看作一个整体，它是一个新的随机变量，因为它是 Z 的函数。根据方差的标准定义 \text{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2，我们可以得到：

\text{Var}(\mathbb{E}[Y|Z]) = \mathbb{E}[(\mathbb{E}[Y|Z])^2] - (\mathbb{E}[\mathbb{E}[Y|Z]])^2

再次应用迭代期望定律于上式的第二项：

\mathbb{E}[\mathbb{E}[Y|Z]] = \mathbb{E}[Y]

代入后，步骤二的最终表达式为：

\text{Var}(\mathbb{E}[Y|Z]) = \mathbb{E}[(\mathbb{E}[Y|Z])^2] - (\mathbb{E}[Y])^2

将步骤一和步骤二的结果相加：

\mathbb{E}[\text{Var}(Y|Z)] + \text{Var}(\mathbb{E}[Y|Z]) = \left( \mathbb{E}[Y^2] - \mathbb{E}[(\mathbb{E}[Y|Z])^2] \right) + \left( \mathbb{E}[(\mathbb{E}[Y|Z])^2] - (\mathbb{E}[Y])^2 \right)

我们可以清晰地看到，中间项 \mathbb{E}[(\mathbb{E}[Y|Z])^2] 一正一负，正好相互抵消。于是等式右边简化为：

\mathbb{E}[Y^2] - (\mathbb{E}[Y])^2

而这个表达式正是随机变量 Y 的方差 \text{Var}(Y) 的定义。因此我们证明了全方差公式：

\text{Var}(Y) = \mathbb{E}[\text{Var}(Y|Z)] + \text{Var}(\mathbb{E}[Y|Z])

这个不等式是Rao-Blackwell方法的核心。假设我们想要估计的量是 \mathbb{E}[Y]，标准的蒙特卡洛方法是直接对 Y 进行采样，其估计器的方差与 \text{Var}(Y) 成正比。然而如果我们能找到另一个相关的随机变量 Z，并且可以解析地计算出条件期望 \mathbb{E}[Y|Z]，我们就可以构造一个改进的估计器。这个新估计器首先对 Z 进行采样，然后计算出 \mathbb{E}[Y|Z] 的值，最后对这些值求平均。

由于 \mathbb{E}[\mathbb{E}[Y|Z]] = \mathbb{E}[Y]，这个新估计器仍然是无偏的。但它的方差是 \text{Var}(\mathbb{E}[Y|Z])，根据不等式，这个方差必然小于或等于原估计器的方差。其本质思想是，通过将问题中的一部分不确定性用解析积分，即求条件期望的方式消除掉，我们减少了最终估计中的随机波动，用同样多的样本就能获得更精确的估计结果。

方差缩减的进阶技术

控制变量法

在前面我们学习了Rao-Blackwell定理提供的方差缩减思路，现在我们来看另一个同样强大的方差缩减技术——控制变量法。这个方法的核心思想是利用一个我们已知期望值的辅助随机变量来校正我们真正想要估计的量，从而降低估计的方差。

假设我们的目标是估计某个函数的期望 \mathbb{E}[g(x)]，其中 x 服从某个分布 p(x)。在标准的蒙特卡洛方法中，我们会从分布 p(x) 中抽取样本，然后计算 g(x) 的样本均值。但现在我们假设存在另一个函数 f(x)，它与 g(x) 相关，并且关键是我们已经知道它的期望值 \mathbb{E}[f(x)] = \mu。这个已知期望的函数 f(x) 就被称为控制变量。

有了这个控制变量，我们可以构造一个新的随机变量 W：

W = g(x) + \alpha(f(x) - \mu)

这里 \alpha 是一个待定的实数参数。这个构造背后有深刻的意义，让我们先验证这个新随机变量 W 的期望值。根据期望的线性性质：

\mathbb{E}[W] = \mathbb{E}[g(x)] + \alpha(\mathbb{E}[f(x)] - \mu)

由于我们已知 \mathbb{E}[f(x)] = \mu，因此上式中的第二项 \alpha(\mathbb{E}[f(x)] - \mu) 恰好等于零，所以：

\mathbb{E}[W] = \mathbb{E}[g(x)]

这个结果非常重要：无论参数 \alpha 取什么值，新构造的随机变量 W 的期望始终等于我们想要估计的原始期望 \mathbb{E}[g(x)]。换句话说，用 W 的样本均值来估计 \mathbb{E}[g(x)] 仍然是无偏的。

无偏性保证了估计的正确性，但我们真正关心的是：这个新的估计器能否比原来的估计器具有更小的方差？

我们计算 W 的方差，根据方差的性质，对于随机变量 X 和 Y 以及常数 a， b，有：

\text{Var}(aX + bY) = a^2\text{Var}(X) + b^2\text{Var}(Y) + 2ab\text{Cov}(X,Y)

将这个公式应用到 W = g(x) + \alpha(f(x) - \mu) 上。注意到 f(x) - \mu 本身就是一个随机变量，因为 \mu 是常数，且 \text{Var}(f(x) - \mu) = \text{Var}(f(x))，我们得到：

\text{Var}(W) = \text{Var}(g(x)) + \alpha^2 \text{Var}(f(x)) + 2\alpha \text{Cov}(g(x), f(x))

方差 \text{Var}(W) 是关于参数 \alpha 的二次函数。如果我们能够选择一个合适的 \alpha 值，就有可能让 \text{Var}(W) 小于原始的 \text{Var}(g(x))。

现在的问题是：我们应该如何选择参数 \alpha 才能使方差 \text{Var}(W) 尽可能小？这是一个标准的优化问题。这是一个关于 \alpha 的二次函数，并且二次项系数 \text{Var}(f(x)) 是正的，因此这个抛物线开口向上，必然存在唯一的最小值点。要找到这个最小值点，我们对 \alpha 求导并令导数为零：

\frac{d \text{Var}(W)}{d\alpha} = 2\alpha \text{Var}(f(x)) + 2\text{Cov}(g(x), f(x)) = 0

解这个方程，得到最优的参数值：

\alpha = -\frac{\text{Cov}(g(x), f(x))}{\text{Var}(f(x))}

这个公式给出了控制变量法中参数 \alpha 的最优取值。当我们把这个最优的 \alpha 代入方差表达式时，可以验证得到的最小方差必然满足：

\text{Var}(W) \leq \text{Var}(g(x))

这个不等式确认了控制变量法的有效性：通过引入一个已知期望的辅助变量 f(x) 并选择合适的参数 \alpha，我们构造的新估计器 W 的方差总是不大于原始估计器的方差。更进一步，只有当 g(x) 和 f(x) 完全不相关，即协方差为零时，等号才成立。只要两者存在一定的相关性，我们就能获得方差的实际降低。

重要性采样

到目前为止，我们讨论的蒙特卡洛方法都有一个基本前提：我们需要能够从目标分布 p(x) 中直接抽取样本。但在很多实际问题中，从目标分布 p(x) 中抽样可能非常困难，甚至是不可行的。此时重要性采样提供了一个解决方案：我们可以从一个更容易采样的分布中抽取样本，然后通过适当的权重调整来修正这种采样分布偏离带来的影响。更重要的是，即使我们能够从 p(x) 中采样，按照 p(x) 进行采样也不一定是最优或最必要的策略。根据具体问题的特点，选择一个合适的替代分布往往能够显著提高估计效率。

回到我们的核心目标：估计函数 h(x) 在分布 p(x) 下的期望值 \mathbb{E}_p[h(x)]。根据期望的定义，这个量可以写成积分形式：

\mathbb{E}_p[h(x)] = \int h(x) p(x) dx

现在重要性采样的核心技巧是引入一个新的分布 g(x)，我们称它为提议分布或重要性分布。我们在积分中同时乘以和除以 g(x)，这个操作不会改变积分的值（只要 g(x) \neq 0），但会改变积分的形式：

\mathbb{E}_p[h(x)] = \int h(x) \frac{p(x)}{g(x)} g(x) dx

现在这个积分可以被理解为函数 h(x) \frac{p(x)}{g(x)} 在分布 g(x) 下的期望。换句话说：

\mathbb{E}_p[h(x)] = \mathbb{E}_g\left[h(x) \frac{p(x)}{g(x)}\right]

这个等式告诉我们：我们可以从更容易采样的分布 g(x) 中抽取样本 x_1， x_2， \dots， x_N，然后用这些样本来估计原本需要从 p(x) 采样才能得到的期望值。对于期望，我们依旧使用平均值来进行估算，估计器形式为：

\hat{I}_N = \frac{1}{N} \sum_{i=1}^{N} \frac{p(x_i)}{g(x_i)} h(x_i)= \frac{1}{N} \sum_{i=1}^{N} W_i h(x_i)

其中 W_i = \frac{p(x_i)}{g(x_i)} 被称为重要性权重。

这个权重修正了由于我们从错误的分布 g(x) 而不是正确的分布 p(x) 进行采样所带来的偏差。当某个样本点 x_i 在目标分布 p(x) 下比在提议分布 g(x) 下更可能出现时，即 p(x_i) > g(x_i)，权重 W_i > 1 会增加这个样本的贡献；反之则减少其贡献。

重要性采样有一个至关重要的前提条件：提议分布 g(x) 的支撑集必须包含目标分布 p(x) 的支撑集，用数学语言表达为：

\text{supp}(g) \supseteq \text{supp}(p)

这里分布的支撑集指的是该分布的概率密度为正值的所有点的集合。这个条件的物理含义是：对于任何在目标分布 p(x) 下可能出现的点，即 p(x) > 0 的点，提议分布 g(x) 在该点处也必须有正的概率密度，即 g(x) > 0。如果违反这个条件，就会存在某些点 x 满足 p(x) > 0 但 g(x) = 0，这会导致重要性权重 W = \frac{p(x)}{g(x)} 变为无穷大，整个估计过程就会失败。

在重要性采样中，样本 \{x_m\}_{m=1}^N 可以通过归一化的权重来处理。具体来说，即使我们只知道 p(x) 和 g(x) 的未归一化形式，比如只知道它们各自乘以一个未知常数，我们仍然可以使用归一化的重要性权重来进行估计。这在贝叶斯推断中特别有用，因为后验分布的归一化常数往往是未知的。

提议分布选择的关键性

重要性采样的效果在很大程度上取决于提议分布 g 的选择，一个糟糕的 g 不仅不能提高效率，反而可能导致估计器方差爆炸，使得方法完全失效。

那么什么样的 g 是最优的呢？理论上，最优的提议分布选择能够最小化重要性采样估计器的方差。然而这里存在一个悖论：要确定最优的 g，我们需要知道期望值 \mathbb{E}_g[h(x)] = I，但这个 I 恰恰是我们试图通过重要性采样来估计的未知量。

随机抽样的实际实现

从理论到实践

到目前为止，我们讨论的所有蒙特卡洛方法都依赖于一个基本前提：我们能够从目标分布 p(x) 中生成样本。理论上我们需要样本序列 \{x_1， x_2， \dots， x_N\}，使得当样本数量 N 趋于无穷时，样本均值收敛到真实期望：

\frac{1}{N} \sum_{n=1}^{N} h(x_n) \rightarrow \mathbb{E}_p[h(x)]

这个收敛的前提是样本 x_n 必须服从分布 p(x)，即 x_n \sim p(x)。但是在实际的计算机实现中，我们如何才能真正从分布中抽样呢？这涉及到随机数生成。

均匀分布的模拟

所有随机抽样的基础是均匀分布。在区间 [0，1] 上的均匀分布 \mathcal{U}_{[0，1]} 是所有其他随机分布的构建基础。通过适当的变换方法，如逆变换采样、拒绝采样等，我们可以从均匀分布出发，生成服从任意其他分布的样本。因此模拟均匀分布是实现所有随机算法的第一步。

在实际计算中，我们使用的是伪随机数生成器。这些生成器并不是真正意义上的随机，而是通过确定性的算法产生一个看起来像随机的数列。具体来说，伪随机数生成器的工作方式是：从一个初始值 x_0 开始，这个初始值被称为种子，然后通过反复应用某个函数 f 来生成一个序列：

x_0, \quad f(x_0), \quad f(f(x_0)), \quad \dots

也就是说，第 n 个生成的数是 x_n = f^{(n)}(x_0)，即对种子连续应用函数 f 共 n 次的结果。

伪随机数生成器的核心目标是：生成的数列 \{x_n\} 的经验分布应该尽可能接近均匀分布 \mathcal{U}_{[0，1]}。换句话说，虽然这些数是通过确定性算法生成的，但它们在统计意义上应该表现得像是从真正的均匀分布中独立抽取的。

最常用且标准的均匀随机数生成器是线性同余生成器。这类生成器的递推公式非常简洁：

x_{i+1} = (ax_i + c) \mod m

这个公式中包含三个整数参数：a 是乘数，c 是增量，m 是模数。从一个初始种子 x_0 开始，每次通过这个线性变换加上取模运算来生成下一个数。取模运算保证了所有生成的数都落在 \{0， 1， 2， \dots， m-1\} 这个有限集合中。通过将这些整数除以 m，我们得到近似均匀分布在 [0，1] 区间上的数。

线性同余生成器的优点是计算速度极快，实现简单。但它的局限性在于生成的序列必然是周期性的，因为状态空间是有限的，并且如果参数选择不当，会表现出明显的统计相关性。因此参数 a、c、m 的选择需要经过精心设计，以保证生成序列的统计质量。

从均匀分布到任意分布的变换

一旦我们能够生成均匀分布 \mathcal{U}_{[0，1]} 的样本，如何将这些均匀分布的样本变换成我们需要的任意目标分布的样本？

逆变换采样法基于累积分布函数的逆变换。假设我们有一个随机变量 X 服从概率密度函数 f(x)，其对应的累积分布函数为 F(x)：

F(x) = \mathbb{P}(X \leq x) = \int_{-\infty}^{x} f(t) dt

逆变换采样法的核心思想是利用累积分布函数的广义逆函数，记作 Q(u)，它被定义为：

Q(u) := \inf\{x: F(x) \geq u\}

对于给定的概率值 u \in [0，1]，Q(u) 是使得累积分布函数 F(x) 达到或超过 u 的最小的 x 值。这个定义即使在 F 不连续或不严格单调的情况下也是有效的，这就是为什么它被称为广义逆。

如果累积分布函数 F 是连续且严格单调递增的，那么它存在通常意义上的逆函数 F^{-1}，并且这个逆函数也是连续的且构成双射。在这种理想情况下，广义逆 Q(u) 就退化为标准的逆函数 F^{-1}(u)。

逆变换采样法

逆变换采样法的定理表明，如果我们有一个随机变量 X 服从概率密度函数 f(x)，其累积分布函数为 F，对应的广义逆函数为 Q，同时我们有一个均匀分布的随机变量 U \sim \mathcal{U}(0，1)，那么通过广义逆变换得到的随机变量 Q(U) 就服从目标分布 g，即 Q(U) \sim g。

对于一个 [0，1] 区间上的均匀随机数 U，我们在累积分布函数曲线上找到纵坐标为 U 的点，然后该点对应的横坐标 x = F^{-1}(U) 就是服从目标分布的样本。

在连续情况下证明这个定理

假设 U \sim \mathcal{U}(0，1)，我们定义随机变量 Y = F^{-1}(U)。我们的目标是证明 Y 服从目标分布 g，也就是说 Y 的累积分布函数应该等于 F(a)（概率密度函数为 g 的分布）。

对于任意实数 a，我们计算 Y 的累积分布函数：

F_Y(a) = \mathbb{P}(Y \leq a) = \mathbb{P}(F^{-1}(U) \leq a)

这里的关键步骤是利用累积分布函数 F 的性质。由于在连续情况下 F 是严格单调递增的，我们可以对不等式两边同时应用 F，且不等号方向不变。因此：

\mathbb{P}(F^{-1}(U) \leq a) = \mathbb{P}(F(F^{-1}(U)) \leq F(a))

由于 F 和 F^{-1} 互为逆函数，我们有 F(F^{-1}(U)) = U，所以上式简化为：

\mathbb{P}(U \leq F(a))

由于 U 服从 [0，1] 上的均匀分布，对于任何 0 \leq u \leq 1，都有 \mathbb{P}(U \leq u) = u。因此：

\mathbb{P}(U \leq F(a)) = F(a)

这就证明了 Y 的累积分布函数 F_Y(a) 等于目标分布的累积分布函数 F(a)，从而 Y \sim g。证明完毕。

这个证明的精妙之处在于它巧妙地利用了均匀分布的特殊性质：均匀分布的累积分布函数就是恒等函数，即 \mathbb{P}(U \leq u) = u。正是这个性质使得逆变换法能够完美地将均匀分布转换为任意目标分布。

指数分布的实例

让我们通过一个具体的例子来理解逆变换采样法。考虑参数为 \lambda 的指数分布，其概率密度函数为：

g(x) = \lambda \exp(-\lambda x), \quad x \geq 0

指数分布用于描述随机事件之间的等待时间等。参数 \lambda 控制了分布的衰减速度，\lambda 越大，分布越集中在原点附近。

要应用逆变换采样法，我们首先需要推导指数分布的累积分布函数。根据定义，累积分布函数是概率密度函数从负无穷到 x 的积分。对于指数分布，由于 x < 0 时密度为零，我们只需要从 0 积分到 x：

F(x) = \int_0^x \lambda \exp(-\lambda t) dt = 1 - \exp(-\lambda x), \quad x \geq 0

累积分布函数 F(x) = 1 - \exp(-\lambda x) 是一个从 0 单调递增到 1 的函数，满足作为累积分布函数的所有要求。

下一步是计算这个累积分布函数的逆函数 F^{-1}(u)。计算逆函数的标准步骤：我们设 u = F(x)，然后解出 x 关于 u 的表达式。从 u = 1 - \exp(-\lambda x) 出发，我们移项得到 \exp(-\lambda x) = 1 - u。对两边取自然对数得到 -\lambda x = \log(1-u)，最后解出：

x=F^{-1}(u) = -\frac{1}{\lambda} \log(1-u)

在实际应用中，算法的流程非常简单。我们首先生成一个服从 [0，1] 均匀分布的随机数 u，然后将这个 u 代入逆函数公式 F^{-1}(u) = -\frac{1}{\lambda} \log(1-u)，计算得到的 x 就是服从参数为 \lambda 的指数分布的样本。

逆变换采样法的限制

逆变换采样法的前提是我们能够获得累积分布函数的逆函数 F^{-1} 或其广义逆 Q。这个要求看似简单，但在实践中却是一个严格的限制。许多常见的概率分布，如正态分布、伽马分布等，的累积分布函数并没有解析形式的逆函数表达式。对于这些分布，虽然我们可以用数值方法近似计算逆函数，但这会显著增加计算成本，使得逆变换法不再那么高效。

因此逆变换法主要适用于那些逆函数有显式表达式的分布，比如我们前面看到的指数分布。

即使对于那些有显式逆函数的分布，逆变换法在计算机实现时也容易受到数值误差的影响。这主要体现在两个方面：

首先，计算逆函数时可能涉及对数、指数等超越函数的运算，这些运算本身就存在舍入误差

其次，当累积分布函数在某些区域变化非常陡峭或非常平缓时，即使输入的均匀随机数 u 有微小扰动，输出的 x = F^{-1}(u) 也可能产生较大偏差。这种数值不稳定性在处理分布的尾部区域时尤为明显。

因此逆变换法主要适用于实值标量随机变量。

并且当我们需要生成多维随机向量 X \in \mathbb{R}^d 时，逆变换法不能直接推广。对于高维情况，一种可行的策略是通过逐次积分的方式来处理：我们可以先对联合分布的某些变量积分得到边缘分布，然后再处理条件分布，将高维采样问题分解为一系列一维采样问题。但这种分解往往依赖于变量之间的具体依赖关系，并不总是可行或高效的。这就是为什么在高维情况下，我们通常需要借助其他更复杂的采样方法，比如马尔可夫链蒙特卡洛方法。

二维采样的策略

在开始介绍变量变换方法之前，我们先看一个重要的二维采样例子。对于二维随机向量 (x_1， x_2)，其联合概率密度可以通过条件概率和边缘概率进行分解:

f_{x_1,x_2}(x_1,x_2) = f_{x_1|x_2}(x_1) f_{x_2}(x_2)

这个分解公式为二维采样提供了一个实用的策略。我们不需要直接从二维联合分布中采样，而是可以分两步进行: 第一步从边缘分布 f_{x_2}(x_2) 中采样得到 x_2，第二步在已知 x_2 的条件下从条件分布 f_{x_1|x_2}(x_1) 中采样得到 x_1。这样我们就将一个二维采样问题转化为了两个一维采样问题的序列。这种分解的关键在于能够计算或表示出边缘分布和条件分布，如果这两个分布都是我们熟悉且容易采样的形式，那么二维采样就变得可行了。

变量变换方法

基本原理

变量变换方法也是为了能从均匀分布生成其他分布样本，其核心思想是通过对已知分布的随机变量进行适当的代数变换来获得目标分布的样本。

假设我们有一个服从标准均匀分布的随机变量 U \sim \mathcal{U}(0，1)，现在想要通过变换来生成其他区间上的均匀分布。第一个变换是简单的平移: X = m + U。这个变换将 [0，1] 区间整体向右或向左平移 m 个单位，因此得到的 X 服从 [m， m+1] 区间上的均匀分布，即 X \sim \mathcal{U}(m， m+1)。

第二个变换是尺度变换: X = \delta U。这个变换将 [0，1] 区间拉伸（当 \delta > 1 ）或压缩（当 0 < \delta < 1 ） \delta 倍，得到 [0， \delta] 区间上的均匀分布，即 X \sim \mathcal{U}(0， \delta)。参数 \delta 直接控制了新区间的长度。

如果我们想要生成任意区间 [a，b] 上的均匀分布 X \sim \mathcal{U}(a，b)，可以结合平移和拉伸两种操作:

X = a + (b-a)U

这个公式的几何意义是:首先用系数 (b-a) 将 [0，1] 区间拉伸到长度为 (b-a) 的区间 [0， b-a]，然后通过加上 a 将整个区间平移到起点为 a 的位置，最终得到区间 [a，b]。

一般变换公式

假设我们有一个随机变量 X 服从概率密度函数 f_X(x)，现在通过一个函数变换 Y = \phi(X) 得到新的随机变量 Y。那么 Y 的概率密度函数 f_Y(y) 是什么?

当变换函数 \phi 是单调的，即严格单调递增或严格单调递减时，我们有一个重要的变换公式:

f_Y(y) = \left|\frac{dy}{d(\phi^{-1}(y))}\right|^{-1} f_X(\phi^{-1}(y))

首先是原分布在逆变换点处的密度 f_X(\phi^{-1}(y))，其次是变换的雅可比行列式在一维情况下就是导数的倒数 |\phi'(\phi^{-1}(y))|^{-1}，最后是绝对值符号，它确保密度函数始终为正，无论变换是递增还是递减。

证明的思路是:先写出 Y 的累积分布函数 F_Y(y) = \mathbb{P}(Y \leq y) = \mathbb{P}(\phi(X) \leq y)，然后利用 \phi 的单调性将其转化为关于 X 的概率，接着对 y 求导并应用链式法则，就能得到上述密度变换公式。

对于最常见的仿射变换，即线性变换加平移 Y = aX + b，上述一般公式可以简化为:

f_Y(y) = \frac{1}{|a|} f_X\left(\frac{y-b}{a}\right)

系数 \frac{1}{|a|} 是因为变换改变了坐标的尺度，需要用这个因子来保证概率密度的归一化。

余弦变换的例子

假设 U \sim \mathcal{U}(0，\pi) 是标准均匀分布，我们定义 Y = \cos(U)。求 Y 的概率密度函数

应用变换公式，我们需要计算余弦函数的导数和逆函数。在区间 [0，\pi] 上，余弦函数是严格单调递减的，其逆变换是 u = \arccos(y)。将这些代入变换公式:

f_Y(y) = \frac{1}{\pi} \frac{1}{\sqrt{1-y^2}}, \quad y \in [-1,1]

这个结果说明，虽然输入的 U 是均匀分布，但经过余弦变换后，Y 的分布变成了一个U形的密度函数。其形状如下: 在 y = \pm 1 附近，即余弦函数变化最缓慢的地方，密度趋于无穷;而在 y = 0 附近，即余弦函数变化最快的地方，密度最小。这个分布被称为反正弦分布或Arcsin分布，分母中的 \sqrt{1-y^2} 因子正是导致两端密度爆炸的原因。

经典的变换例子

通过变量变换方法，我们可以从均匀分布的样本出发，通过巧妙的变换构造出许多重要分布的样本。这里介绍两个非常经典且实用的例子。

第一个例子展示如何从均匀分布生成伽马分布。假设我们有 a 个独立的标准均匀分布随机变量 u_1， u_2， \dots， u_a，每个都服从 \mathcal{U}(0，1)。我们定义变换:

Y = \frac{a}{\beta} \sum_{j=1}^{a} \log(u_j)

那么这个随机变量 Y 服从参数为 (a， \beta) 的伽马分布，即 Y \sim Ga(a， \beta)。

首先对于每个均匀分布的随机变量 u_j，-\log(u_j) 服从参数为 1 的指数分布（不证明）。

其次独立指数分布随机变量的和服从伽马分布（不证明），这是伽马分布的一个重要性质。因此 \sum_{j=1}^{a} (-\log(u_j)) 服从 Ga(a， 1)，再乘以缩放因子 \frac{a}{\beta}，就得到了 Ga(a， \beta)。

第二个例子展示如何生成贝塔分布。假设我们有 a+b 个独立的均匀分布随机变量，定义变换:

Y = \frac{\sum_{j=1}^{a} \log(u_j)}{a+b \sum_{j=1}^{b} \log(u_j)}

那么 Y \sim Be(a，b)，即服从参数为 (a，b) 的贝塔分布。

这个构造利用了贝塔分布与伽马分布之间的关系。如果有两个独立的伽马分布随机变量 X_1 \sim Ga(a， \theta) 和 X_2 \sim Ga(b， \theta)，注意它们有相同的尺度参数，那么比值 \frac{X_1}{X_1 + X_2} 服从 Be(a，b)。

这些经典方法需要生成多个均匀随机数并进行对数运算和求和，计算成本相对较高。其次这里的方法有一个重要的限制条件: 参数 a 和 b 必须是正整数，即 a， b \in \mathbb{N}^*。这是因为这些构造方法依赖于对有限个随机变量进行求和，而求和的项数必须是整数。当伽马分布或贝塔分布的形状参数不是整数时，这些简单的变换方法就不再适用，我们需要寻求其他的采样技术。

Box-Müller变换

Box-Müller变换是从均匀分布生成正态分布样本的一个经典方法。这个方法由George Box和Mervin Müller在1958年提出，至今仍被广泛使用。

这个算法是一个精确的算法，不依赖于任何近似（不是通过中心极限定理），意思是我们不需要通过对大量随机变量求和来近似正态分布，而是通过一个精确的数学变换直接得到严格服从正态分布的样本。

Box-Müller方法的核心是一个巧妙的变换。假设我们有两个独立的标准均匀分布随机变量 U_1， U_2 \sim \mathcal{U}(0，1)，通过以下变换可以得到两个独立的标准正态分布随机变量:

X_1 = \sqrt{-2\log(U_1)} \cos(2\pi U_2) \sim \mathcal{N}(0,1)

X_2 = \sqrt{-2\log(U_1)} \sin(2\pi U_2) \sim \mathcal{N}(0,1)

这两个公式本质上是在做极坐标到直角坐标的转换。U_2 通过乘以 2\pi 被转换为角度，而 -2\log(U_1) 被转换为径向距离的平方。对数变换确保了径向分量服从正确的分布，即卡方分布，而角度的均匀分布确保了方向的各向同性。

从两个均匀分布的样本，我们一次性生成了两个独立的正态分布样本。这种一次生成两个样本的特性使得方法在计算效率上有优势。X_1 和 X_2 不仅各自服从标准正态分布，而且它们是相互独立的，这是通过正弦和余弦的正交性保证的。

接受-拒绝采样法

方法的动机与核心思想

到目前为止，我们学习的逆变换法和变量变换法都有各自的局限性。现在我们介绍一种更加通用的采样方法——接受-拒绝采样法，它能够在很多其他方法失效的情况下依然有效。

接受-拒绝采样法的提出是为了解决以下几种实际困境。

第一种情况是我们无法获得目标分布的显式形式。这里的无显式形式具体指的是：我们没有概率密度函数 g 的解析表达式，或者无法计算累积分布函数 F，更不用说求出其逆函数 F^{-1} 了。在贝叶斯推断中，这种情况非常常见，因为后验分布的表达式包含一个难以计算的归一化常数，即证据 p(y)，因此我们只能得到后验分布的未归一化形式。

第二种情况是，即使我们知道分布的数学形式，也可能不存在直接的采样算法。逆变换法要求累积分布函数有解析的逆函数，变量变换法需要找到合适的变换关系，这不总是可行的。

面对这些困难，接受-拒绝采样法提供了一个巧妙的解决方案：我们不直接从困难的目标分布 g 中采样，而是利用另一个更简单、更容易采样、或者我们已经掌握了的辅助分布。这个辅助分布被称为提议分布或候选分布。

等价性定理

接受-拒绝采样法建立在一个等价性定理之上：从一维分布 f 中采样，等价于从一个二维区域中进行均匀采样。该定理的核心陈述是：从概率密度函数 f(x) 中采样一个随机变量 X，等价于从如下定义的二维区域中均匀采样一个点对 (X, u)：

(X, u) \sim \mathcal{U}\left(\{(x, u): 0 < u < f(x)\}\right)

换句话说，我们考虑二维平面上所有满足条件 0 < u < f(x) 的点构成的区域。这个区域正是概率密度函数曲线 f(x) 下方的面积。定理表明，如果我们在这个区域内均匀地抽取点对 (X, u)，那么 X 的边缘分布恰好就是我们想要的目标分布 f。

从几何角度理解：一维的概率分布 f(x) 被提升到了二维空间。x 轴上某个点 x_0 处的密度值 f(x_0) 决定了在该点上方可以取到的 u 值的范围 [0, f(x_0)]。密度越大的地方，u 方向的范围越大，因此在二维均匀采样时，该 x 值被选中的概率也越大。这正是一维概率密度的本质含义。

基于这个等价性定理，接受-拒绝采样法的核心想法就自然浮现了。我们在一个比目标密度曲线 f(x) 更高的某个区域中模拟二维点对 (X, u)，然后只保留那些满足约束条件 u < f(x) 的样本点。具体来说，如果我们能够找到一个更简单的、高度更大的包络函数来覆盖目标密度 f(x)，我们就可以在这个更大的区域内轻松地生成候选样本点。对于每个生成的候选点 (x, u)，我们检查它是否落在目标密度曲线 f(x) 下方，即是否满足 u < f(x)。如果满足，我们接受这个样本的 x 坐标；如果不满足，我们拒绝它并重新生成新的候选点。通过这个生成-检验-筛选的过程，最终保留下来的 x 样本就精确地服从目标分布 f。

贝塔分布的例子

为了更具体地理解接受-拒绝采样法，让我们看一个实际例子：从贝塔分布中采样。图示展示了如何为贝塔分布构造一个简单的包络区域。贝塔分布的支撑集是有界的，记为区间 D，通常是 [0,1]。我们首先找到贝塔分布密度函数的最大值，记为 m = \max(\text{Beta})。这个最大值的位置取决于贝塔分布的参数。

基于这个最大值 m，我们可以构造一个高度为 m、底边为支撑区间 D 的矩形区域。这个矩形完全包含了贝塔密度曲线下方的区域。密度函数 f(x) 下方的目标区域将会被接收，，黄色斜线部分的候选样本会被拒绝。

接受-拒绝采样的实现过程非常直接。我们在矩形区域 [0,1] \times [0,m] 内均匀地生成二维点对 (x, u)：

(x, u) \sim \mathcal{U}_{[0,1] \times [0,m]}(x, u)

这个均匀采样可以通过分别生成两个独立的均匀随机数来实现：x 从 [0,1] 上均匀采样，u 从 [0,m] 上均匀采样。接下来我们保留那些满足条件 0 < u < f(x) 的样本。换句话说，我们检查生成的点 (x,u) 是否落在目标密度函数曲线下方。如果 u < f(x)，这个点在曲线下方，我们就接受它，保留其 x 坐标作为有效样本；如果 u \geq f(x)，这个点在曲线上方，我们就拒绝它，重新生成新的候选点。

接受-拒绝采样法的一个基本要求是我们需要能够计算目标密度函数 f(x) 的值。虽然我们不需要知道 f(x) 的归一化常数，即不需要知道精确的概率密度函数，但我们必须能够对任意给定的 x 值计算出与 f(x) 成比例的某个量，以便判断候选点是否应该被接受。

正确性证明

现在我们需要从数学上严格证明接受-拒绝采样法的正确性，也就是说通过这个生成-检验-筛选的过程得到的样本确实精确地服从目标分布 f。证明的核心是计算 X 的累积分布函数，并验证它恰好等于目标分布的累积分布函数。

假设目标密度函数 f 在区间 [a,b] 上有界，即存在一个有限的常数 m 使得对所有 x \in [a,b] 都有 f(x) \leq m。我们在矩形区域 [a,b] \times [0,m] 内均匀地生成样本对 (Y, u) \sim \mathcal{U}_{[a,b] \times [0,m]}，其中 Y 是候选样本的横坐标，u 是用于判断接受或拒绝的纵坐标。我们保留那些满足约束条件 0 < u < f(Y) 的样本，将满足条件的 Y 值记作 X。

对于任意实数 x，X 的累积分布函数为：

\mathbb{P}(X < x) = \mathbb{P}(Y < x \mid u < f(Y))

X < x 的概率等于在接受条件 u < f(Y) 成立的前提下，Y < x 的条件概率。这正是接受-拒绝采样的核心思想：我们只关心那些被接受的样本。根据条件概率的定义，这可以改写为：

\mathbb{P}(X < x) = \int_a^x f(y \mid u < f(y)) dy

为了计算这个条件密度，我们需要利用贝叶斯法则。根据贝叶斯法则，条件密度函数可以表示为：

f(y|u) = \frac{f(y,u)}{\int_a^b f(y,u) dy}

在我们的设定中，(Y,u) 在矩形区域 [a,b] \times [0,m] 上服从均匀分布，因此它们的联合密度是一个常数：

f(y,u) = \frac{1}{(b-a) \times m}

将这个联合密度代入到累积分布函数的计算中，利用条件概率公式和联合分布的性质，累积分布函数可以表示为：

= \int_a^x \int_0^{f(y)} f(y, u) \, du \, dy \times \frac{1}{\int_a^b \int_0^{f(y)} f(y, u) \, dy \, du}

我们需要计算在区间 [a,x] 上、密度曲线 f(y) 下方的区域内均匀采样的概率，然后除以在整个区间 [a,b] 上、密度曲线 f(y) 下方的区域内均匀采样的总概率作为归一化。由于 (Y,u) 的联合密度是常数 \frac{1}{(b-a) \times m}，这个常数在乘积和积分中都出现，可以提取出来。

现在关键的一步是计算内层积分。对于内层积分 \int_0^{f(y)} du，由于被积函数是常数 1，积分结果就是积分上限 f(y)。因此上式化简为：

= \frac{\int_a^x f(y) dy}{\int_a^b f(y) dy}

现在我们需要处理分母。由于 f(y) 是一个概率密度函数，根据概率分布的归一化条件，密度函数在整个支撑区间上的积分必须等于 1：

\int_a^b f(y) dy = 1

这是概率密度函数的基本性质。将这个结果代入，我们得到：

= \int_a^x f(y) dy = F(y \leq x)

这正是目标分布 f 的累积分布函数。根据微积分的基本定理，我们还知道：

\frac{\partial}{\partial x} \int_a^x f(y) dy = f(x)

这意味着对累积分布函数求导就得到概率密度函数本身。这就完整地证明了：通过接受-拒绝采样法得到的随机变量 X 的累积分布函数恰好是 F(y \leq x)，因此 X 精确地服从目标分布 f。整个证明的精妙之处在于，虽然我们从一个更大的矩形区域中均匀采样，但通过接受-拒绝机制的筛选，最终保留下来的样本在 x 方向上的分布恰好按照目标密度函数 f(x) 的形状分布。

一般化的接受-拒绝采样

前面我们看到的接受-拒绝采样使用矩形作为包络区域，但这种方法在很多情况下效率较低。现在我们来看一个更加通用且高效的接受-拒绝采样定理，它允许我们使用任意的提议分布作为包络。

假设 f(x) 和 g(x) 都是概率密度函数，并且存在一个常数 \Pi > 1 使得对所有 x 都满足：

f(x) < \Pi \times g(x)

目标密度函数 f(x) 在任何点处都被放大后的提议密度 \Pi \times g(x) 所覆盖。常数 \Pi 被称为包络常数，它衡量了提议分布需要被放大多少倍才能完全覆盖目标分布。

白色曲线是目标密度 f(x)，黄色的三角形包络线表示放大后的提议密度 \Pi \times g(x)，它完全覆盖了目标密度。

在这个条件下，为了从目标分布 f 中采样，我们只需要按照以下步骤生成样本。首先从提议分布 g 中生成一个候选样本 Y \sim g。然后在给定这个 Y 的条件下，从区间 [0, \Pi g(Y)] 上生成一个均匀随机数 U|Y \sim \mathcal{U}(0, \Pi g(Y))。接下来我们重复这个生成过程，直到找到一个满足接受条件 0 < u < f(Y) 的样本对。一旦找到这样的样本对，我们就保留对应的 Y 值作为从目标分布 f 中抽取的有效样本。

具体算法流程如下。

步骤1：从提议分布中生成候选样本 Y \sim g，同时生成一个标准均匀随机数 U \sim \mathcal{U}(0,1)。

步骤2：检查接受条件。如果满足 U \leq \frac{f(Y)}{\Pi g(Y)}，则接受这个候选样本，令 X = Y 并输出。

步骤3：如果步骤2的条件不满足，即候选样本被拒绝，则返回步骤1，重新生成新的候选样本，继续尝试。

这个算法的核心在于接受概率 \frac{f(Y)}{\Pi g(Y)}。由于我们已经保证了 f(x) < \Pi \times g(x)，因此这个比值总是小于1，可以被正确地解释为一个概率。当候选点 Y 位于目标密度 f 较高的区域时，这个比值较大，接受的可能性就高；反之，在目标密度较低的区域，接受概率也相应较低。

算法的重要性质

第一个重要的理论保证是：这个算法会在有限时间内终止，且终止的概率为1。这意味着虽然理论上我们可能需要尝试多次才能找到一个被接受的样本，但算法必然会在有限步内成功。更进一步，算法期望需要等待的迭代次数恰好等于 \Pi

这个结果可以从几何分布的性质得出。当我们生成一个候选样本 Y \sim g 和相应的均匀随机数 U 时，该样本被接受的概率恰好是 \frac{1}{\Pi}。这是因为接受条件 U < \frac{f(Y)}{\Pi g(Y)} 在 Y 的分布下平均满足的概率就是 \frac{1}{\Pi}，而这些尝试是独立的，因此直到第一次成功所需的尝试次数服从几何分布。几何分布的期望值是成功概率的倒数，因此期望迭代次数为 \frac{1}{1/\Pi} = \Pi。因此包络常数 \Pi 越接近1，即提议分布 g 越接近目标分布 f，算法效率越高，期望迭代次数越少。

接受-拒绝采样法有一个严重的限制：维度诅咒。当我们需要从高维分布中采样，即 X \in \mathbb{R}^d 时，情况会急剧恶化。在 d 维空间中，样本的接受概率变为：

\frac{1}{\Pi^d}

接受概率随着维度 d 呈指数级下降。即使包络常数 \Pi 只是略大于1，比如 \Pi = 2，在10维空间中接受概率就只有 \frac{1}{2^{10}} \approx 0.001，意味着平均需要约1000次尝试才能得到一个有效样本。这就是为什么接受-拒绝采样法在高维问题中通常不可行，我们需要寻求其他更适合高维空间的采样方法，比如马尔可夫链蒙特卡洛方法。

马尔可夫链

蒙特卡洛估计的前提条件

在开始学习马尔可夫链之前，我们需要先回顾一下如何用蒙特卡洛方法估计期望值。具体来说，假设我们想要计算某个函数 h(x) 关于随机变量 x 的期望值 I：

I = \mathbb{E}[h(x)] = \int h(x) f(x) \, dx

这个积分往往非常难以直接计算，尤其当维度很高或者积分区域复杂的时候。蒙特卡洛方法的核心思想是：如果我们能够从目标分布 f(x) 中抽取 N 个独立同分布的样本 x_1, x_2, \ldots, x_N，那么我们就可以用这些样本的函数值的算术平均来近似真实的期望值：

\hat{I} = \frac{1}{N} \sum_{n=1}^{N} h(x_n), \quad x_n \sim f(x)

这个估计量有三个重要性质保证了蒙特卡洛方法的有效性。

一致性：当样本数量 N 趋向于无穷大时，估计量 \hat{I}_N 会依概率收敛到真实值 I。

无偏性：估计量的期望值恰好等于真实值，即 \mathbb{E}[\hat{I}_N] = I。虽然单次估计可能有误差，但重复多次估计并取平均时，这个平均值会围绕真实值波动。

估计的分布特性：当 N 足够大时，估计误差经过标准化后会服从标准正态分布：

\frac{\hat{I} - I}{\sqrt{V_{\hat{I}}}} \sim \mathcal{N}(0, 1)

然而这里有一个关键的前提条件：我们需要能够从目标分布 f(x) 中抽取独立同分布的样本，但是我们有些时候很难抽取样本，因此我们才有了上面的蒙特卡洛法、重采样法，本质上都是避开对原目标分布的直接采样。但是逆变换法要求我们必须写出累积分布函数的解析形式，并且能够计算逆函数，变量变换法中，不一定能找到变换公式且高纬空间下计算雅可比行列式太难算，接收拒绝法的问题在于维度灾难。

因此我们重新审视一下问题，之前的所有方法都试图直接生成独立同分布的样本，每个样本的生成都是独立的，互不影响。这种独立性在理论分析时很方便，但在实际采样时却成了沉重的负担——它要求我们必须完全理解目标分布，要么知道如何精确地变换到它，要么知道如何有效地包络它。我们是否可以放松这个独立性要求？我们真正需要的并不是独立的样本，而是来自目标分布的样本，如果我构造的样本序列并非独立的，但是随着序列的演进，样本分布逼近目标分布，那么在迭代足够多次后，后续的样本就可以近似看作从目标分布中抽取的。

换句话说，我们不再试图一次性生成一个完美的样本，而是迭代的方式，将样本序列逐步收敛到目标分布，从长期来看，后续样本会按照目标分布在空间中分布。

这个样本序列就是马尔可夫链。

马尔可夫链的基本思想

马尔可夫链方法提供了一个巧妙的解决方案来应对无法直接从目标分布采样的困境。它的基本思想是：既然我们无法直接生成服从分布 f(x) 的独立样本，那么我们可以构造一个随机变量序列 X^{(0)}, X^{(1)}, X^{(2)}, \ldots, X^{(m)}，使得这个序列的分布逐渐趋向于我们的目标分布 f(x)。形象地说，这就像一个渐进的过程：

X^{(0)} \rightarrow X^{(1)} \rightarrow X^{(2)} \rightarrow \ldots \rightarrow X^{(m)} \sim f(x)

这里的关键是，当迭代次数 m 足够大时，X^{(m)} 的分布就会非常接近目标分布 f(x)。虽然序列中前面几个变量的分布可能与目标分布相差甚远，但只要我们运行足够长的时间，最终得到的样本就可以近似看作是从目标分布中抽取的。

马尔可夫链的数学定义

一个马尔可夫链是一个随机变量序列 (X^{(0)}, X^{(1)}, \ldots, X^{(n)})，这个序列的概率分布会随着时间或者说随着迭代步数的推移而演化。马尔可夫链的特征可以用三个核心要素来完整刻画。

第一个要素是转移核。转移核描述的是从序列的一个状态如何转移到下一个状态的概率机制。数学上，转移核定义为给定前一个状态 X^{(n-1)} 的条件下，下一个状态 X^{(n)} 的条件概率密度函数：

K\big(X^{(n-1)}, X^{(n)}\big) = f_{X^{(n)}|X^{(n-1)}}\big(X^{(n)} \mid X^{(n-1)}\big)

这个转移核是马尔可夫链的核心，它完全决定了序列如何从一个状态演化到下一个状态。可以把它理解成一个规则或机制：当我们处于状态 X^{(n-1)} 时，这个规则告诉我们下一个状态 X^{(n)} 应该如何随机生成，以及生成每个可能值的概率密度是多少。

第二个要素是马尔可夫性，也叫作条件独立性质。这是马尔可夫链最重要的特征，它说的是：下一个状态只依赖于当前状态，而与更早的历史状态无关。用数学语言表达就是：

f\big( X^{(n)} \mid X^{(n-1)}, X^{(n-2)}, \ldots, X^{(0)} \big) = f\big( X^{(n)} \mid X^{(n-1)} \big)

等式左边是给定所有历史信息，从初始状态 X^{(0)} 到当前状态 X^{(n-1)}，的条件下下一个状态的条件分布。等式右边则只依赖于当前状态 X^{(n-1)}。这个等式告诉我们，当前状态已经包含了预测未来所需要的所有信息，更早的历史状态不再提供额外的预测价值。这就是所谓的无记忆性。这个性质大大简化了马尔可夫链的分析，因为我们不需要追踪整个历史轨迹，只需要知道当前状态就够了。

第三个要素是初始分布。由于马尔可夫链是一个动态演化的过程，我们需要指定这个过程从哪里开始，也就是序列的第一个随机变量 X^{(0)} 服从什么分布：

f^0 = f\big( X^{(0)} \big)

初始分布决定了链的起点。不同的初始分布会导致链在初期有不同的行为，但在许多情况下，当链满足某些良好性质时，无论从哪里开始，最终链的分布都会收敛到同一个稳态分布，这就是我们想要的目标分布 f(x)。

随机游走：最简单的马尔可夫链

图中展示的是一个典型的随机游走过程，黄色和白色两条轨迹代表了两个不同的实现路径。随机游走是理解马尔可夫链的绝佳起点，因为它的演化规则极其简单：

X^{(n)} = X^{(n-1)} + \varepsilon^{(n)}

这个公式描述了状态的演化机制：第 n 步的位置等于前一步的位置加上一个随机扰动。这里的 \varepsilon^{(n)} 是一个服从正态分布的随机变量：

\varepsilon^{(n)} \sim \mathcal{N}(0, \sigma^2)

这个扰动项 \varepsilon^{(n)} 独立于过去的所有状态，它的均值为零，方差为 \sigma^2。正是这个独立性保证了马尔可夫性质的成立。为了理解这一点，我们来看看下一个状态是如何依赖于历史的。给定所有历史状态 X^{(0)}, X^{(1)}, \ldots, X^{(n-1)} 时，X^{(n)} 的条件分布完全由 X^{(n-1)} 和独立的噪声 \varepsilon^{(n)} 决定。由于 \varepsilon^{(n)} 与过去无关，因此 X^{(n)} 的条件分布只依赖于 X^{(n-1)}，这正是马尔可夫性质。

在图中我们能观察到，虽然两条轨迹从相近的位置出发，但随着时间推移它们的路径完全不同。黄色轨迹整体呈上升趋势，而白色轨迹则在零附近震荡。这种行为完全符合随机游走的特征：每一步的方向都是随机的，长期行为具有很大的不确定性。

离散有限状态的马尔可夫链

现在我们转向一个更加结构化的情形：离散有限状态的马尔可夫链。在这种情况下，随机变量不再在连续空间中取值，而是只能取有限个离散状态。我们将状态空间记为：

E = \{\alpha_1, \ldots, \alpha_K\} = \{1, \ldots, q\}

这里 K 或 q 表示状态的总数。为了简化记号，我们通常直接用整数 1, 2, \ldots, q 来标记这些状态。在离散情形下，转移核变成了一个矩阵，称为转移矩阵 K，其中每个元素定义为：

P_{xy} = P\big(X^{(n)} = y \mid X^{(n-1)} = x\big)

这里需要特别注意索引的含义：x 对应矩阵的行，表示当前所处的状态；y 对应矩阵的列，表示将要转移到的状态。因此 P_{xy} 表示从状态 x 转移到状态 y 的概率。

转移矩阵必须满足两个基本性质。首先，由于 P_{xy} 表示概率，它必须非负：

P_{xy} \geq 0 \quad \forall x, y

其次，从任意状态 x 出发，必须转移到某个状态，因此所有可能的转移概率之和必须为 1：

K(x, E) = \sum_{y \in E} P_{xy} = 1

这个条件意味着转移矩阵的每一行的和都等于 1。满足这两个条件的矩阵称为随机矩阵或马尔可夫矩阵。这是一个 q \times q 的方阵，完全刻画了离散马尔可夫链的动力学特性。

随机矩阵的具体例子

为了更好地理解随机矩阵，我们来看一个具体的例子：

K = \begin{pmatrix} \frac{1}{2} & \frac{1}{4} & \frac{1}{4} \\ \frac{1}{4} & \frac{1}{4} & \frac{1}{2} \\ \frac{1}{4} & \frac{1}{4} & \frac{1}{2} \end{pmatrix}

假设这个矩阵描述的是天气状态的转移，状态空间为：

E = \{\text{晴天}, \text{雨天}, \text{雪天}\}

矩阵的第一行 (\frac{1}{2}, \frac{1}{4}, \frac{1}{4}) 告诉我们：如果今天是晴天，那么明天有 \frac{1}{2} 的概率仍是晴天，有 \frac{1}{4} 的概率是雨天，有 \frac{1}{4} 的概率是雪天。我们可以验证每一行的和确实等于 1，这保证了概率的归一化。

概率分布的时间演化

在马尔可夫链中，一个关键问题是：如果我们知道初始时刻的概率分布，如何计算未来任意时刻的概率分布？答案由以下递推关系给出：

p^{(m)} = p^{(m-1)} K

这里 p^{(m)} 是一个行向量，表示时刻 m 时系统处于各个状态的概率分布。这个公式的推导基于全概率公式。设系统在时刻 m 处于状态 y 的概率为 p_y^{(m)}，那么：

p_y^{(m)} = \sum_{x \in E} p_x^{(m-1)} P_{xy}

要计算时刻 m 处于状态 y 的概率，我们需要考虑所有可能的路径——从时刻 m-1 的每个状态 x 出发，乘以从 x 转移到 y 的概率，然后求和。写成矩阵形式就是 p^{(m)} = p^{(m-1)} K。

通过递推，我们可以得到更一般的结果：

p^{(m)} = p^{(0)} K^m

这个公式揭示了一个深刻的事实：第 m 步的概率分布完全由初始分布 p^{(0)} 和转移矩阵的 m 次幂决定。

数值计算与收敛行为

让我们用前面的天气转移矩阵进行具体计算。假设初始分布是均匀分布：

p^{(0)} = \left(\frac{1}{3}, \frac{1}{3}, \frac{1}{3}\right)

第一步后的分布为：

p^{(1)} = p^{(0)} K = \left(\frac{1}{3}, \frac{1}{4}, \frac{5}{12}\right)

继续迭代计算 p^{(2)}，我们发现：

p^{(2)} = \left(\frac{1}{3}, \frac{1}{4}, \frac{5}{12}\right)

可见，p^{(2)} = p^{(1)}，分布不再变化了，这意味着系统已经达到了平衡状态。实际上，如果我们继续迭代，会发现：

p^{(\infty)} = \left(\frac{1}{3}, \frac{1}{4}, \frac{5}{12}\right)

平稳分布的概念与性质

我们将这个不再变化的分布记为：

p^* = \left(\frac{1}{3}, \frac{1}{4}, \frac{5}{12}\right)

并称之为平稳分布或平衡分布。平稳分布的定义特征是它满足方程：

p = pK

如果系统的当前分布是 p，那么经过一步转移后，分布仍然是 p。换句话说，平稳分布是转移算子的不动点。从线性代数的角度看，平稳分布是转移矩阵转置的特征值为 1 对应的左特征向量（经过归一化）。

对于离散有限状态的马尔可夫链，有一个重要的定理保证了平稳分布的存在性和唯一性。这个定理告诉我们：任何有限状态的马尔可夫链都至少存在一个平稳分布。在满足某些额外条件（如不可约性和非周期性）时，这个平稳分布是唯一的，而且无论从什么初始分布开始，系统的分布都会收敛到这个唯一的平稳分布。

平稳分布的存在性可以通过Brouwer不动点定理来证明。由于概率分布构成的单纯形是紧凸集，而转移操作是连续映射，因此必然存在不动点。唯一性则需要额外的条件来保证，这涉及到马尔可夫链的遍历理论。

连续状态空间的马尔可夫链

从离散状态转向连续状态是一个重要的扩展，因为在实际中我们经常需要处理连续分布的采样问题。在连续状态空间中，转移核不再是一个矩阵，而是一个条件概率密度函数。

在连续情况下，转移核 K(x', x) 表示从状态 x' 转移到状态 x 的条件概率密度。我们也可以用另一个记号 T(x \mid x') 来表示这个条件密度，强调它是给定 x' 条件下 x 的分布。他们的关系是：

K(x', x) = T(x \mid x')

作为一个条件概率密度函数，转移核必须满足归一化条件：对于状态空间 E 中的任意状态 x'，以 x 为变量对转移核积分必须等于1：

\forall x' \in E, \quad \int T(x \mid x') \, dx = 1

给定我们当前处于状态 x'，下一个状态 x 必然会落在整个状态空间的某个位置，所以对所有可能的 x 的概率密度积分应该等于1。这是概率密度函数的基本要求。

如果我们知道时刻 m-1 的概率密度函数 f^{(m-1)}(x')，那么时刻 m 的密度函数 f^{(m)}(x) 是什么？这个问题的答案类似于离散情况，但现在求和变成了积分。密度演化的公式为：

f^{(m)}(x) = \int T(x \mid x') \, f^{(m-1)}(x') \, dx'

或者等价地写成：

f^{(m)}(x) = \int K(x', x) \, f^{(m-1)}(x') \, dx'

时刻 m 系统处于状态 x 附近的概率密度，等于对所有可能的前一时刻状态 x' 进行积分：前一时刻处于 x' 的密度 f^{(m-1)}(x') 乘以从 x' 转移到 x 的条件密度 K(x', x)。这就是全概率公式在连续情况下的体现。这个积分对所有可能的来源状态 x' 进行加权平均，权重是该来源状态的概率密度。

类似于离散情况的平稳分布，在连续情况下我们有不变分布的概念。一个概率密度函数 f^*(x) 被称为不变分布，如果它满足：

f^*(x) = \int K(x', x) \, f^*(x') \, dx'

如果系统的密度函数是 f^*(x)，那么经过一步转移后密度函数仍然是 f^*(x)。这是离散情况下 p = pK 的连续版本。找到不变分布就是要解这个积分方程，这通常比离散情况更困难，因为涉及到积分而不是简单的矩阵运算。

随机游走的详细分析

为了更好地理解连续马尔可夫链的行为，我们用随机游走这个经典例子来分析。考虑如下的随机游走过程：

X^{(n)} = X^{(n-1)} + \varepsilon^{(n)}, \quad \varepsilon \sim \mathcal{N}(0, \nu)

这里 \varepsilon^{(n)} 是均值为0、方差为 \nu 的正态分布随机变量。假设初始分布为标准正态分布：

p^{(0)} = \mathcal{N}(0,1)

首先我们需要确定这个过程的转移核。给定当前状态 X^{(n-1)} = x'，下一个状态是 X^{(n)} = x' + \varepsilon^{(n)}，其中 \varepsilon^{(n)} \sim \mathcal{N}(0, \nu)。这意味着给定 x'，X^{(n)} 的分布是均值为 x'、方差为 \nu 的正态分布。因此转移核为：

K(x', x) = T(x \mid x') = \frac{1}{\sqrt{2\pi \nu}} \exp\left( -\frac{1}{2} \frac{(x - x')^2}{\nu} \right)

这个表达式就是均值为 x'、方差为 \nu 的正态分布的概率密度函数。注意这里 (x - x') 表示新状态相对于旧状态的偏差，这个偏差服从 \mathcal{N}(0, \nu)。

现在我们要计算经过一步转移后的分布 p^{(1)}(x)。根据密度演化公式：

p^{(1)}(x) = \int K(x', x) \, p^{(0)}(x') \, dx'

将转移核和初始分布的具体表达式代入：

p^{(1)}(x) = \int \frac{1}{\sqrt{2 \pi \nu}} \exp\left( -\frac{1}{2} \frac{(x - x')^2}{\nu} \right) \frac{1}{\sqrt{2\pi}} \exp\left( -\frac{1}{2} (x')^2 \right) dx'

我们可以把 X^{(1)} 看作两个独立正态随机变量的和：X^{(1)} = X^{(0)} + \varepsilon^{(1)}，其中 X^{(0)} \sim \mathcal{N}(0,1)，\varepsilon^{(1)} \sim \mathcal{N}(0,\nu)。根据正态分布的可加性，两个独立正态随机变量的和仍然是正态分布

p^{(1)} = \mathcal{N}(0, 1 + \nu)

这意味着经过一步随机游走后，分布仍然是正态分布，但变得更加分散了。

通过类似的推理，我们可以得出一般性的结论。在每一步中，我们都加上一个方差为 \nu 的独立正态随机变量。经过 m 步后，X^{(m)} 就是初始变量加上 m 个独立的正态扰动：

X^{(m)} = X^{(0)} + \varepsilon^{(1)} + \varepsilon^{(2)} + \cdots + \varepsilon^{(m)}

由于 X^{(0)} \sim \mathcal{N}(0,1)，每个 \varepsilon^{(i)} \sim \mathcal{N}(0,\nu)，根据正态分布的可加性：

p^{(m)} = \mathcal{N}\big(0, 1 + m\nu\big)

随机游走的分布方差会随着时间线性增长。这意味着随机游走会越来越分散，不会收敛到某个固定的分布。这个例子说明，并不是所有的马尔可夫链都有有意义的平稳分布。对于随机游走，由于方差持续增长，不存在一个有限方差的不变分布。这也说明了为什么我们需要对马尔可夫链施加某些条件，如不可约性和常返性，才能保证平稳分布的存在和收敛性。

马尔可夫链的收敛性理论

核心问题是：我们如何确保链会收敛到目标分布？这就需要研究马尔可夫链的收敛性。

对于离散有限状态的马尔可夫链，收敛性问题相对容易处理。理论告诉我们，平稳分布具有三个良好的性质：它存在、它是唯一的在满足一定条件下，并且我们可以通过解线性方程组来找到它。具体来说，我们需要求解方程组：

p = pK

这是一个齐次线性方程组，再加上概率归一化条件 \sum_i p_i = 1，我们就可以用标准的线性代数方法，如高斯消元法，来求解平稳分布。这种方法在理论上是完美的，但在实践中，当状态空间很大时，直接求解可能计算量很大。

除了直接求解线性方程组，还有另一种更加直接的策略：反复应用转移矩阵。也就是说，从任意一个初始分布 p^{(0)} 出发，不断计算 p^{(1)} = p^{(0)} K，p^{(2)} = p^{(1)} K，直到 p^{(m)} = p^{(m-1)} K。如果链满足一定的条件，那么当 m 足够大时，p^{(m)} 会收敛到平稳分布 p^*。这种方法的优势在于它模拟了实际运行马尔可夫链的过程，而且不需要求解大型线性方程组。

但是不是所有的马尔可夫链都会收敛到平稳分布。为了保证收敛性，链必须满足一些重要的性质。两个最关键的性质是不可约性和常返性。

不可约性描述的是链的连通性：从任意一个状态出发，都有可能在有限步内到达任意另一个状态。

P_x(\tau_y < \infty) > 0, \quad \forall x, y \in E

这里 \tau_y 表示首次到达状态 y 的时刻，即 \tau_y = \min\{n \geq 1 : X^{(n)} = y\}。条件 P_x(\tau_y < \infty) > 0 的含义是：从状态 x 出发，在有限时间内到达状态 y 的概率大于0。这个条件必须对所有状态对 (x,y) 都成立。

不可约性的直观含义是：状态空间中没有相互隔离的孤岛。如果链是不可约的，那么从任何地方出发，经过足够多步的随机游走，都有可能访问到状态空间的任何其他地方。这个性质对于收敛性至关重要，因为如果链不是不可约的，那么不同的状态子集可能有不同的平稳分布，唯一性就不成立了。

常返性描述的是链的另一个重要性质：链不会永远逃离某个状态。用直观的话说，如果一个状态是常返的，那么链一旦访问过这个状态，将来一定会以概率1再次回到这个状态，并且会无限多次回访。为了刻画常返性，我们定义一个状态 \omega 的访问次数。考虑随机变量 \eta_\omega，它计算链在整个运行过程中访问状态 \omega 的总次数：

\eta_\omega = \sum_{n=1}^{\infty} \mathbb{1}_\omega(X^{(n)})

这里 \mathbb{1}_\omega(X^{(n)}) 是指示函数：当 X^{(n)} = \omega 时取值1，否则取值0。所以 \eta_\omega 就是把所有时刻访问状态 \omega 的次数加起来。

现在关键问题是，这个访问次数的期望值是多少？根据期望值的大小，我们可以对状态进行分类。如果期望访问次数是无穷大，即 \mathbb{E}[\eta_\omega] = \infty，那么状态 \omega 被称为常返的。这意味着链会无限多次地回到这个状态，这是一个好的性质。相反，如果期望访问次数是有限的，即 \mathbb{E}[\eta_\omega] < \infty，那么这个状态是暂态的。暂态意味着链可能只访问这个状态有限次，然后就再也不回来了。

对于一个不可约的马尔可夫链，要么所有状态都是常返的，要么所有状态都是暂态的。如果链是不可约且所有状态都是常返的，再加上非周期性，那么链一定会收敛到唯一的平稳分布，而且这个收敛与初始状态无关。这就是马尔可夫链蒙特卡洛方法的理论基础。

Harris常返性

在前面我们讨论了不可约性和常返性这两个重要性质。现在我们要引入一个更强的概念：Harris常返性。Harris常返性是这两个性质的结合和强化。定义是：如果一个马尔可夫链是不可约的，并且其所有状态都是常返的，那么这个链就被称为Harris常返的。

因为在连续状态空间中，简单的常返性概念变得复杂且不够用。Harris常返性提供了一个更强、更适合连续情况的条件。它保证了链不仅会无限次访问每个状态的邻域，而且这种访问具有某种均匀性。在离散情况下，不可约加常返基本等价于Harris常返。但在连续情况下，Harris常返性是一个更强的条件，它确保了链在状态空间中的彻底混合。

非周期性

除了不可约性和常返性，还有第三个重要的性质需要满足，那就是非周期性。对于一个状态 i，我们考虑从这个状态出发经过多少步可以回到自己，定义集合 \{n \geq 1 \mid P^n(i,i) > 0\}，这个集合包含所有可能的回归时间：即经过 n 步从状态 i 回到 i 的概率大于0的那些时间步 n。状态 i 的周期定义为这个集合的最大公约数：

d(i) = \operatorname{pgcd}\{ n \ge 1 \mid P^n(i, i) > 0 \}

通过一个例子来理解周期的含义。假设从状态 i 出发，可能在第2步、第4步、第6步回到 i，但不可能在奇数步回到 i。那么可能的回归时间集合是 \{2, 4, 6, 8, \ldots\}，这个集合的最大公约数是2，所以状态 i 的周期是 d(i) = 2。这意味着链的行为有一种周期性模式：从 i 出发只能在偶数步回到 i。这种周期性会阻碍收敛，因为链的分布会在不同的时刻呈现不同的模式，不会稳定下来。

如果一个状态的周期 d(i) = 1，那么这个状态被称为非周期的。周期为1意味着最大公约数是1，也就是说存在某个足够大的时刻 n_0 之后，对于所有 n \geq n_0 都有 P^n(i,i) > 0。换句话说，经过足够多步后，每一步都有可能回到状态 i，不存在周期性的限制。

如果链中的所有状态都是非周期的，那么整个马尔可夫链被称为非周期的。对于不可约链，实际上只需要检查一个状态的周期性，如果一个状态是非周期的，那么所有状态都是非周期的。非周期性是保证收敛到平稳分布的重要条件，因为它排除了分布在不同时刻之间周期性振荡的可能性。

遍历性

现在我们可以定义遍历性了，这是马尔可夫链理论中最重要的概念之一。如果一个马尔可夫链同时满足Harris常返性和非周期性，那么被称为遍历的。这两个条件的结合产生了强大的收敛保证。

当一个链是遍历的时，我们有如下的收敛定理。设 \pi 是链的平稳分布，那么无论从什么初始分布 p^{(0)} 出发，经过 m 步转移后得到的分布 K^m(x, \cdot) 都会收敛到 \pi。用数学语言表达，这个收敛在全变差距离意义下成立：

\lim_{m \to \infty} \left\| \int K^m(x, \cdot)\, p^{(0)}(x)\, dx - \pi(\cdot) \right\|_{TV} = 0

左边的积分 \int K^m(x, \cdot)\, p^{(0)}(x)\, dx 表示从初始分布 p^{(0)} 出发，经过 m 步转移后得到的分布。这个分布与目标平稳分布 \pi 的全变差距离 \|\cdot\|_{TV}，当 m 趋向无穷时收敛到0。全变差距离是衡量两个概率分布差异的一种方式，当它趋向0时，意味着两个分布在各个方面都越来越接近。

遍历性给我们带来三个重要的结论。

第一，存在唯一的平稳分布 \pi。唯一性非常重要，因为它保证了不同的方法、不同的初始条件都会导向同一个目标分布。

第二，链会收敛到这个平稳分布 \pi。这意味着只要我们运行链足够长的时间，我们得到的样本分布会任意接近目标分布。

第三，这个收敛与初始分布 p^{(0)} 的选择无关。无论我们从哪里开始，最终都会到达同一个地方。这个性质使得马尔可夫链蒙特卡洛方法非常稳健：我们不需要精心选择初始状态，只要运行足够长的时间，链自然会忘记它的起点并收敛到目标分布。

期望的收敛

遍历性保证的分布收敛进一步导致期望的收敛，这对于我们最初的目标——估计期望值——至关重要。如果我们想要估计某个函数 h(X) 关于平稳分布 \pi 的期望 \mathbb{E}_{\pi}[h(X)]，那么遍历性告诉我们，使用链在时刻 n 的状态 X^{(n)} 计算的期望 \mathbb{E}_{p^{(n)}}[h(X^{(n)})] 会收敛到真实期望：

\lim_{n \to +\infty} \left| \mathbb{E}_{p^{(n)}}[h(X^{(n)})] - \mathbb{E}_{\pi}[h(X)] \right| = 0

当链运行足够长时间后，n 很大时，X^{(n)} 的分布已经非常接近 \pi，所以用 X^{(n)} 计算的期望也会非常接近真实的期望值。这就把分布的收敛转化为我们真正关心的期望值的收敛。

这个结果为马尔可夫链蒙特卡洛方法提供了理论保证：只要我们等待足够长的时间，让链达到平稳状态，我们就可以用链上的样本来估计期望值，而这个估计会是准确的。

遍历定理：时间平均与空间平均的等价性

遍历定理建立了时间平均和空间平均之间的联系。定理表述如下：如果马尔可夫链 (X^{(n)}) 是遍历的，那么沿着链的轨迹对函数 h 求时间平均，其极限等于 h 关于平稳分布的空间平均，即期望：

\lim_{n \to +\infty} \frac{1}{n} \sum_{i=1}^{n} h(X_i) = \int h(x)\, \pi(x)\, dx

等式左边是时间平均：我们运行马尔可夫链，得到一系列状态 X_1, X_2, \ldots, X_n，对每个状态计算函数值 h(X_i)，然后求这些函数值的算术平均。这是一个沿着时间轴的平均。

等式右边是空间平均：对整个状态空间按照平稳分布 \pi 加权平均函数 h，这就是期望值 \mathbb{E}_{\pi}[h(X)]。

遍历定理告诉我们，这两种看似不同的平均方式在极限情况下给出相同的结果。这个结果意味着我们可以通过运行一条单一的马尔可夫链轨迹，用时间平均来估计期望值，而不需要生成多个独立的样本。

更深层次地说，遍历定理体现了遍历系统的一个基本特征：系统在长时间的演化过程中会访遍状态空间的各个区域，并且访问的频率与平稳分布成正比。单个轨迹的长时间行为包含了整个系统的统计性质。这也解释了为什么我们不需要从平稳分布生成大量独立样本，而是可以用一条长链来达到同样的估计效果。

需要注意的是，遍历定理成立的前提是链确实是遍历的，即Harris常返且非周期。如果这些条件不满足，时间平均可能不收敛，或者收敛到错误的值。因此验证或确保马尔可夫链的遍历性是至关重要的。

平衡条件：构造平稳分布的工具

最后我们介绍一个重要的充分条件，它可以用来验证某个分布是否为平稳分布。如果一个分布 \pi 和转移核 K 满足如下等式：

\pi(x) K(x, x') = \pi(x') K(x', x)

那么 \pi 一定是这个马尔可夫链的平稳分布。

等式左边 \pi(x) K(x, x')：表示处于状态 x 的概率乘以从 x 转移到 x' 的条件概率，这给出了从 x 流向 x' 的联合概率流。
等式右边 \pi(x') K(x', x)：表示处于状态 x' 的概率乘以从 x' 转移到 x 的条件概率，这给出了从 x' 流向 x 的联合概率流。

平衡条件说的是：从 x 到 x' 的概率流与从 x' 到 x 的概率流相等。这种双向的平衡保证了系统的稳态：如果系统按照分布 \pi 分布，那么任意两个状态之间的概率流是平衡的，没有净的概率流动，因此整体分布不会改变。

为什么成立呢？对上式对 x 积分（或离散情形求和）：

\int \pi(x)K(x, x')\,dx = \pi(x') \int K(x', x)\,dx = \pi(x')

其中 \int K(x', x)\,dx = 1

因此左右积分求和的结果都是 \pi(x')

需要特别注意：平衡条件只是平稳分布的充分条件，不是必要条件。满足平衡条件的分布一定是平稳分布，但平稳分布不一定满足平衡条件。更重要的是，平衡条件本身并不能告诉我们关于唯一性或收敛性的任何信息。即使我们验证了某个分布满足平衡条件，我们还需要检验其他条件，如不可约性、非周期性，来确保链会收敛到这个分布。

尽管如此，平衡条件在实践中非常有用，因为它提供了一个相对容易验证的充分条件。在设计马尔可夫链蒙特卡洛算法时，比如Metropolis-Hastings算法，我们经常通过构造满足平衡条件的转移核来确保目标分布是平稳分布。

平衡条件给了我们一个清晰的设计原则：只要我们能设计出一个转移核，使得任意两个状态之间的概率流是平衡的，我们就知道这个马尔可夫链会以目标分布作为平稳分布。这个思想在接下来学习具体的MCMC算法时会反复出现。

马尔可夫链蒙特卡洛方法：Gibbs采样器

Gibbs采样器的基本思想

在前面的学习中，我们已经掌握了马尔可夫链的理论基础，特别是遍历性如何确保链最终会收敛到唯一的平稳分布。但现在摆在我们面前的实际问题是：当我们想要从某个复杂的目标分布中采样时，应该如何具体构造这样一个马尔可夫链呢？更进一步说，我们要如何设计转移核，使得这个链的平稳分布恰好就是我们想要采样的那个目标分布？Gibbs采样器正是回答这个问题的一个巧妙方案。

让我们先建立问题的数学框架。假设我们面对的是一个 P 维的随机向量

X = [X_1, \ldots, X_P]

我们的目标分布记为 \pi(X)，这是一个多维的联合分布。在许多实际应用中，这个联合分布可能非常复杂，直接从中采样几乎是不可能的。比如在贝叶斯推断中，\pi(X) 可能是一个高维的后验分布，它涉及多个相互关联的参数，没有解析的表达式，也没有现成的采样算法。

Gibbs采样器的思路是这样的：虽然联合分布很复杂，但如果我们固定其他所有变量，只看某一个变量的条件分布，这个条件分布往往会简单很多，甚至可能是我们熟悉的标准分布，那么从中采样并不困难。为了表达这个想法，我们需要引入一个记号。对于向量 X 的第 k 个分量 X_k，我们用 X_{-k} 表示除了 X_k 之外的所有其他分量

X_{-k} = [X_1, \ldots, X_{k-1}, X_{k+1}, \ldots, X_P]

那么，X_k 在给定 X_{-k} 条件下的条件分布就记为

\pi_k(X_k \mid X_{-k})

Gibbs采样器的策略就是：不直接从联合分布采样，而是轮流从每个变量的条件分布中采样。通过这种逐个更新的方式，我们构造出一条马尔可夫链，而这条链的平稳分布恰好就是我们想要的联合分布 \pi(X)。

算法的具体流程

现在让我们把Gibbs采样器的算法步骤具体化。整个算法分为两个阶段。第一个阶段是初始化。我们需要给出一个起始点 X^{(0)} = [X_1^{(0)}, X_2^{(0)}, \ldots, X_P^{(0)}]。这个初始值可以随意选择，比如从均匀分布中随机抽取，也可以根据我们对问题的了解给出一个合理的猜测。根据马尔可夫链的遍历性理论，只要链满足遍历条件，无论从哪里出发，最终都会收敛到同一个平稳分布。当然，初始值的选择会影响收敛的快慢：如果初始值恰好在目标分布的高概率区域附近，收敛会快一些，否则，可能需要更多的迭代才能进入目标分布的典型集。

第二个阶段是迭代更新。对于第 m 次迭代（m 从1到 N），我们要依次更新每一个分量。当轮到更新第 k 个分量时，我们从下面这个条件分布中采样

X_k^{(m)} \sim \pi_k\big(X_k \mid X_1^{(m)}, \ldots, X_{k-1}^{(m)}, X_{k+1}^{(m-1)}, \ldots, X_P^{(m-1)}\big)

当我们更新 X_k 时，排在 X_k 前面的那些变量（X_1 到 X_{k-1}），我们使用的是当前迭代中已经更新过的新值，它们的上标都是 (m)；而排在 X_k 后面的那些变量（X_{k+1} 到 X_P），我们使用的还是上一轮迭代的旧值，上标是 (m-1)。这样设计原因很简单：当我们按顺序从 X_1 更新到 X_P 时，前面的变量已经被更新了，我们自然应该使用最新的信息；而后面的变量还没轮到更新，所以只能暂时使用它们的旧值。这种更新方式叫做"系统扫描"或"顺序更新"。

二维情况下的几何直观

为了更直观地理解Gibbs采样的行为，我们来看二维的情况。图中的等高线代表目标分布 \pi(X_1, X_2) 的概率密度，等高线越密集的地方概率密度越大。从某个初始点出发，Gibbs采样的轨迹一会儿垂直移动，一会儿水平移动。

当我们更新 X_1 而保持 X_2 固定时，采样点只能沿着平行于 X_1 轴的方向移动，也就是垂直方向。反过来，当我们更新 X_2 而保持 X_1 固定时，采样点只能沿着平行于 X_2 轴的方向移动，也就是水平方向。每次移动的距离和方向不是随机游走，而是根据条件分布来决定的：我们更可能移向在该条件下概率密度更高的位置。

通过这种交替的垂直和水平移动，采样点会逐渐探索整个概率分布的支撑区域。更重要的是，由于每步都倾向于移向高概率区域，采样点在概率密度高的区域会停留更长时间，在概率密度低的区域则快速通过。经过足够长的时间后，采样点访问各个区域的频率将正比于该区域的概率密度，这正是我们想要的采样效果。

为什么目标分布是平稳分布

Gibbs采样器能够工作的理论保证是：目标分布 \pi 确实是这个马尔可夫链的平稳分布。现在我们来严格证明这个结论。我们考虑二维情况，即 X = (X_1, X_2)。高维情况的证明逻辑完全相同，只是符号更繁琐一些。

首先我们要明确Gibbs采样器的转移核是什么。在二维情况下，一次完整的迭代包含两个步骤：先从条件分布 \pi_{X_1|X_2}(X_1 \mid X_2') 中采样得到新的 X_1，然后从条件分布 \pi_{X_2|X_1}(X_2 \mid X_1) 中采样得到新的 X_2。注意第二步中，我们使用的是刚刚采样得到的新 X_1 作为条件。因此，从状态 X' = (X_1', X_2') 转移到状态 X = (X_1, X_2) 的转移核可以写成

T(X \mid X') = \pi_{X_1 \mid X_2}(X_1 \mid X_2') \, \pi_{X_2 \mid X_1}(X_2 \mid X_1)

这是两个条件分布的乘积，体现了先后两步的复合作用。

要证明 \pi 是平稳分布，我们需要验证不变性条件，即如果链当前的分布是 \pi，经过一步转移后分布仍然是 \pi :

\int T(X \mid X') \, \pi(X') \, dX' = \pi(X)

现在让我们开始推导。把转移核和目标分布的表达式代入左边

\int T(X \mid X') \, \pi(X') \, dX' = \int \pi_{X_1 \mid X_2}(X_1 \mid X_2') \, \pi_{X_2 \mid X_1}(X_2 \mid X_1) \, \pi(X_1', X_2') \, dX_1' dX_2'

这是一个对 (X_1', X_2') 的二重积分，我们的策略是利用边缘化和条件概来逐步化简这个积分。

首先观察到，第二个条件分布 \pi_{X_2 \mid X_1}(X_2 \mid X_1) 中不包含积分变量 X_1' 和 X_2'，所以可以提到积分号外面

= \pi_{X_2 \mid X_1}(X_2 \mid X_1) \int_{X_2'} \pi_{X_1 \mid X_2}(X_1 \mid X_2') \left( \int_{X_1'} \pi(X_1', X_2')\, dX_1' \right) dX_2'

现在看内层的积分 \int_{X_1'} \pi(X_1', X_2') dX_1'。这是对联合分布 \pi(X_1', X_2') 关于第一个变量的积分，根据边缘化的定义，这给出的是 X_2' 的边缘分布

\int_{X_1'} \pi(X_1', X_2')\, dX_1' = \pi_{X_2}(X_2')

将这个结果代入，我们得到

= \pi_{X_2 \mid X_1}(X_2 \mid X_1) \int_{X_2'} \pi_{X_1 \mid X_2}(X_1 \mid X_2')\, \pi_{X_2}(X_2')\, dX_2'

接下来处理这个积分 \int_{X_2'} \pi_{X_1 \mid X_2}(X_1 \mid X_2') \pi_{X_2}(X_2') dX_2'。根据条件概率的基本关系，条件分布乘以条件变量的边缘分布等于联合分布，即

\pi_{X_1 \mid X_2}(X_1 \mid X_2') \times \pi_{X_2}(X_2') = \pi(X_1, X_2')

因此这个积分变成了对联合分布关于 X_2' 的边缘化

\int_{X_2'} \pi(X_1, X_2')\, dX_2' = \pi_{X_1}(X_1)

这给出了 X_1 的边缘分布。将这个结果代回，我们有

= \pi_{X_2 \mid X_1}(X_2 \mid X_1) \times \pi_{X_1}(X_1)

最后一步，再次使用条件概率关系，条件分布乘以条件变量的边缘分布等于联合分布

\pi_{X_2 \mid X_1}(X_2 \mid X_1) \times \pi_{X_1}(X_1) = \pi(X_1, X_2) = \pi(X)

这正是我们要证明的结果。因此我们确实验证了

\int T(X \mid X') \, \pi(X') \, dX' = \pi(X)

证明完成。Gibbs采样器之所以保持目标分布不变，正是因为它的转移核很好的保持了目标分布的条件概率结构。当我们从各个条件分布中采样时，这些条件分布本身就是从目标联合分布推导出来的，因此通过适当的组合，我们又回到了原来的联合分布。我们最终证明了如果链从分布 \pi 开始，经过一步Gibbs转移后仍然是分布 \pi。

Gibbs采样器：贝叶斯参数估计

一个具体的统计推断问题

为了真正理解Gibbs采样器如何在实践中发挥作用，我们来看一个经典的贝叶斯统计问题：从观测数据中推断正态分布的参数。假设我们有 n 个观测数据 Y_1, \ldots, Y_n，这些数据是从同一个正态分布中独立抽取的

Y_1, \ldots, Y_n \overset{iid}{\sim} \mathcal{N}(y \mid \mu, \sigma^2)

这里的 \mu 和 \sigma^2 分别是这个正态分布的均值和方差，它们是未知的。在频率学派的框架下，我们会把这两个参数看作固定但未知的常数，通过最大似然估计等方法来求出它们的点估计。但在贝叶斯框架下，我们采取完全不同的视角：我们把这些参数也视为随机变量，它们有自己的概率分布。我们的目标是在看到数据之后，计算出参数的后验分布，从而得到关于参数的完整概率信息。

从数据到参数的后验分布，需要经过三个关键步骤。第一步是写出似然函数，也就是给定参数值时观测到这些数据的概率。由于数据是独立同分布的，联合似然就是各个数据点似然的乘积

p(y \mid \mu, \sigma^2) \propto \prod_{i=1}^n \frac{1}{\sqrt{\sigma^2}} \exp\left( -\frac{1}{2\sigma^2}(y_i - \mu)^2 \right)

这里我们用正比符号 \propto 是因为在贝叶斯推断中，归一化常数会在后面通过后验的归一化自动处理，所以我们可以暂时忽略那些不依赖于参数的常数因子。将这个乘积展开整理，所有的 \frac{1}{\sqrt{\sigma^2}} 合并起来给出 (\sigma^2)^{-n/2}，指数项中的平方和可以提取到一个总的指数中

p(y \mid \mu, \sigma^2) \propto (\sigma^2)^{-n/2} \exp\left( -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu)^2 \right)

也可以把 (\sigma^2)^{-n/2} 写成 \sigma^{-n} 的形式，这在后面会更方便

p(y \mid \mu, \sigma^2) \propto \sigma^{-n} \exp\left( -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu)^2 \right)

第二步是指定先验分布。在这个问题中，我们采用无信息先验，也就是说我们对参数几乎不施加任何先验假设，让数据尽可能地主导后验分布。对于均值 \mu，我们使用平坦先验

p(\mu) = 1

这意味着在看到数据之前，我们认为 \mu 可以是实数轴上的任何值，每个值都同等可能。严格来说这不是一个真正的概率密度（因为它在整个实数轴上积分是无穷大），但在贝叶斯计算中这种不当先验是可以使用的，只要后验分布是正常的。

对于方差 \sigma^2，我们使用Jeffreys先验，这是一种对尺度变换具有不变性的无信息先验

p(\sigma^2) = \frac{1}{\sigma^2}

第三步是应用贝叶斯公式得到后验分布。根据贝叶斯定理，后验分布正比于似然函数乘以先验

p(\mu, \sigma^2 \mid y) \propto p(y \mid \mu, \sigma^2)\, p(\mu)\, p(\sigma^2)

把前面得到的似然和先验表达式代入

p(\mu, \sigma^2 \mid y) \propto \sigma^{-n} \exp\left( -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu)^2 \right) \times 1 \times \frac{1}{\sigma^2}

合并关于 \sigma 的幂次，\sigma^{-n} \times \sigma^{-2} = \sigma^{-(n+2)}

p(\mu, \sigma^2 \mid y) \propto \sigma^{-(n+2)} \exp\left( -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu)^2 \right)

这就是我们的后验分布。虽然这个联合后验分布看起来比较复杂，直接从中采样很困难，但我们可以推导出条件分布，然后使用Gibbs采样器。条件分布 p(\mu \mid \sigma^2, y) 和 p(\sigma^2 \mid \mu, y) 通常具有标准的形式（在这个例子中，前者是正态分布，后者是逆伽马分布），从中采样相对容易。这样，通过Gibbs采样器，我们就能从复杂的联合后验分布中生成样本，用于参数的贝叶斯推断。

推导条件分布

要使用Gibbs采样器，我们需要知道如何从两个条件分布中采样。

条件分布1：给定方差时均值的条件分布

首先我们考虑在已知方差 \sigma^2 和数据 y 的条件下，均值 \mu 的条件分布。从联合后验分布出发，我们固定 \sigma^2，只关注与 \mu 相关的部分：

p(\mu \mid y, \sigma^2) \propto \exp\!\left( -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu)^2 \right)

所有不包含 \mu 的因子（包括 \sigma^{-(n+2)}）都可以看作归一化常数，暂时不管。现在关键是识别这个指数型表达式对应什么分布。指数中的 \sum_{i=1}^n (y_i - \mu)^2 是关于 \mu 的二次型，这个平方和可以分解为

\sum_{i=1}^n (y_i - \mu)^2 = n(\mu - \bar{y})^2 + \sum_{i=1}^n (y_i - \bar{y})^2

其中 \bar{y} = \frac{1}{n}\sum_{i=1}^n y_i 是样本均值。这个分解的第二项不含 \mu，可以归入归一化常数。因此条件分布正比于

\exp\!\left( -\frac{n}{2\sigma^2}(\mu - \bar{y})^2 \right)

这正是正态分布 \mathcal{N}(\bar{y}, \frac{\sigma^2}{n}) 的未归一化形式。因此

p(\mu \mid y, \sigma^2) \sim \mathcal{N}\!\left( \frac{1}{n}\sum_{i=1}^n y_i,\ \frac{\sigma^2}{n} \right)

均值是样本均值 \bar{y} = \frac{1}{m}\sum_{i=1}^m y_i，方差是 \frac{\sigma^2}{m}。这个结果非常直观：给定方差，均值的后验分布以样本均值为中心，方差随着样本量增加而减小（反映了我们对均值估计的不确定性随样本量增加而降低）。从这个正态分布中采样是非常容易的，这就是Gibbs采样器第一步的操作。

条件分布2：给定均值时方差的条件分布

接下来考虑在已知均值 \mu 和数据 y 的条件下，方差 \sigma^2 的条件分布。从联合后验分布出发，固定 \mu，只保留与 \sigma^2 相关的部分：

p(\sigma^2 \mid y, \mu) \propto \sigma^{-n-2}\, \exp\!\left( -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu)^2 \right)

我们可以把 \sigma^{-n-2} 改写成 (\sigma^2)^{-(n/2+1)} 的形式

p(\sigma^2 \mid y, \mu) \propto (\sigma^2)^{-(n/2+1)} \exp\!\left( -\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu)^2 \right)

这个表达式的结构对应逆伽马分布。逆伽马分布 \mathcal{IG}(\alpha, \beta) 的概率密度正比于 x^{-(\alpha+1)} \exp(-\beta/x)。通过对比，我们可以确定参数：\alpha = \frac{n}{2} -1，而 \beta = \frac{1}{2}\sum_{i=1}^n (y_i - \mu)^2。

因此，\sigma^2 的条件分布是逆伽马分布：

p(\sigma^2 \mid y, \mu) \sim \mathcal{IG}\!\left( \frac{m}{2}-1,\ \frac{1}{2}\sum_{i=1}^m (y_i-\mu)^2 \right)

从逆伽马分布采样也有现成的算法。至此，我们得到了Gibbs采样器所需的全部条件分布。实际操作时，我们从某个初始值 (\mu^{(0)}, \sigma^{2(0)}) 出发，然后在每次迭代中我们交替地从这两个分布采样：先固定当前的 \sigma^2 采样新的 \mu，然后固定刚采样的 \mu 采样新的 \sigma^2，如此反复。经过足够多次迭代后，样本就会收敛到联合后验分布 p(\mu, \sigma^2 \mid y)。

Metropolis-Hastings算法：更通用的MCMC框架

Gibbs采样器虽然简洁优雅，但它的应用范围受到一个关键限制：我们必须能够推导出所有的条件分布，并且这些条件分布必须是可以直接采样的。在许多复杂的实际问题中，这个要求很难满足。条件分布可能没有解析形式，或者即使有解析形式也不是标准分布，从中采样同样困难。Metropolis-Hastings算法提供了一个更加通用的解决方案：它只需要我们能够计算目标分布的值（甚至只需要未归一化的值），而不需要知道任何条件分布的形式。

Metropolis-Hastings算法的核心机制是接受-拒绝采样（在马尔可夫链框架下）。算法的基本思路是这样的：我们从某个容易采样的建议分布中生成候选点，然后通过一个精心设计的接受概率来决定是否移动到这个候选点。这个接受概率的设计确保了马尔可夫链的平稳分布恰好是我们想要的目标分布。

算法从一个初始状态 X^{(0)} 开始。然后对于每一次迭代 n，我们执行两个步骤。

第一步是生成候选点，我们从建议分布 q(Y \mid X^{(n)}) 中采样一个新的候选状态

Y^{(n)} \sim q(Y \mid X^{(n)})

建议分布 q(Y \mid X) 描述的是：给定当前状态 X，我们如何提议一个候选状态 Y。这个分布完全由我们自己设计，唯一的要求是它应该容易采样。常见的选择包括在当前点附近进行随机游走，或者使用一些启发式规则来生成候选点。建议分布的设计会显著影响算法的效率，但不会影响最终收敛到的平稳分布，这是算法的一个重要性质。

第二步是接受-拒绝决策。我们不是直接接受候选点 Y^{(n)}，而是计算一个接受概率 \alpha^{(n)}，然后以这个概率接受候选点。具体来说，下一个状态的确定规则是

X^{(n+1)} = \begin{cases} Y^{(n)} & \text{以概率 } \alpha^{(n)} \\ X^{(n)} & \text{以概率 } 1 - \alpha^{(n)} \end{cases}

如果候选点被接受，链移动到新位置 Y^{(n)}；如果被拒绝，链停留在原位置 X^{(n)}。关键问题是：接受概率 \alpha^{(n)} 应该如何设计？

接受概率

接受概率 \alpha^{(n)} 的公式是Metropolis-Hastings算法的核心，它被精心设计以确保目标分布 \pi 是链的平稳分布：

\alpha^{(n)} = \min \left( 1,\ \frac{\pi(Y^{(n)})}{\pi(X^{(n)})} \frac{q(X^{(n)} \mid Y^{(n)})}{q(Y^{(n)} \mid X^{(n)})} \right)

首先，\frac{\pi(Y^{(n)})}{\pi(X^{(n)})} 是目标分布在候选点和当前点的比值。这个比值衡量的是：相对于当前位置，候选位置在目标分布下的"受欢迎程度"。如果候选点的概率密度更高（比值大于1），说明候选点是一个"更好"的位置，我们倾向于接受它。如果候选点的概率密度更低（比值小于1），说明候选点不如当前位置"好"，但我们也不是完全拒绝它，而是以一定概率接受它，因为如果我们只接受"上坡"的移动而拒绝所有"下坡"的移动，链会陷入局部的高概率区域，无法探索整个分布。通过以一定概率接受下坡移动，链保持了足够的随机性，能够全面探索整个概率空间。

第二个比值 \frac{q(X^{(n)} \mid Y^{(n)})}{q(Y^{(n)} \mid X^{(n)})} 是一个校正项，它考虑了建议分布的不对称性。

分子 q(X^{(n)} \mid Y^{(n)}) 表示如果我们在候选点 Y^{(n)}，提议回到当前点 X^{(n)} 的概率密度，分母 q(Y^{(n)} \mid X^{(n)}) 表示从当前点 X^{(n)} 提议到候选点 Y^{(n)} 的概率密度。如果建议分布使得从 X 到 Y 很容易（分母很大），但从 Y 回到 X 很难（分子很小），并且我们不考虑这种不对称性，链会过度地停留在 Y 而不是 X，这样它们在目标分布下的权重就不相同了。校正项的作用就是平衡这种不对称：如果从 X 到 Y 容易而反向困难，我们就降低接受从 X 到 Y 的概率，反之则提高接受概率。这个校正项确保了马尔可夫链满足平衡条件，从而保证目标分布是平稳分布。

最后，外层的 \min(1, \cdot) 操作起到截断的作用。如果括号内的值大于1（意味着无论从概率密度还是提议的对称性来看，移动到候选点都是有利的），我们就以概率1接受，即确定性地接受。如果小于1，我们就以该值本身作为接受概率。这确保了接受概率始终在0到1之间。

这个公式还有一个巨大优势：它只需要计算目标分布的比值，而不需要知道归一化常数。在贝叶斯推断中，后验分布可以写成

\pi(x \mid y) = \frac{p(x, y)}{p(y)} = \frac{1}{K} p(y \mid x) p(x)

这里 K = p(y) 是边缘似然，它是一个对所有参数积分的高维积分，通常难以计算。但在计算比值时

\frac{\pi(Y^{(n)})}{\pi(X^{(n)})} = \frac{p(y \mid Y^{(n)}) p(Y^{(n)}) / K}{p(y \mid X^{(n)}) p(X^{(n)}) / K} = \frac{p(y \mid Y^{(n)}) p(Y^{(n)})}{p(y \mid X^{(n)}) p(X^{(n)})}

归一化常数 K 在分子分母中抵消了。我们只需要计算似然和先验的乘积，这通常是容易做到的。这使得Metropolis-Hastings算法在贝叶斯推断中具有极大的实用价值。

随机游走建议分布

在实际应用中，最简单也最常用的建议分布是基于随机游走的。我们将候选点设为当前点加上一个随机扰动

Y^{(n)} = X^{(n)} + \varepsilon^{(n)}, \quad \varepsilon^{(n)} \sim \mathcal{N}(0, \sigma^2)

扰动项 \varepsilon^{(n)} 服从均值为0、方差为 \sigma^2 的正态分布。这意味着候选点在当前点周围进行局部的随机游走，游走的步长由 \sigma^2 控制。对应的建议分布是

q(Y \mid X^{(n)}) \sim \mathcal{N}(X^{(n)}, \sigma^2)

这是一个以当前点为中心的正态分布。这种建议分布的一个重要性质是对称性：从 X 到 Y 的提议概率等于从 Y 到 X 的提议概率

q(\bar{X} \mid X) = q(X \mid \bar{X})

这是因为正态分布关于均值对称，\mathcal{N}(Y; X, \sigma^2) = \mathcal{N}(X; Y, \sigma^2)。当建议分布是对称的时候，接受概率公式大大简化。校正项 \frac{q(X^{(n)} \mid Y^{(n)})}{q(Y^{(n)} \mid X^{(n)})} 的分子分母相等，比值为1，因此

\alpha = \min \left( 1,\ \frac{\pi(Y^{(n)})}{\pi(X^{(n)})} \right)

校正项消失了，接受概率只取决于目标分布（或后验分布）的比值。这个简化版本就是经典的Metropolis算法，它是Metropolis-Hastings算法在对称建议分布下的特殊情况。

随机游走建议分布虽然简单，但参数 \sigma^2 的选择会显著影响算法的性能。如果 \sigma^2 太小，每步的移动距离很短，链需要很长时间才能探索整个状态空间，收敛会很慢。如果 \sigma^2 太大，候选点经常会跳到概率密度很低的区域（离目标分布的主要质量区域很远），导致很高的拒绝率，链大部分时间停留在原地，同样导致探索缓慢。

Metropolis-Hastings算法的这种 提议 - 接受 / 拒绝 机制说明了：通过容易采样的建议分布来探索状态空间，通过精心设计的接受概率来保证收敛到正确的目标分布。这种机制使得我们能够从几乎任何复杂的高维分布中采样，只要我们能够计算这个分布的未归一化概率密度。这正是MCMC方法在现代统计推断和机器学习中如此重要的原因。

分享

数值统计（完整版）：从极大似然到生成模型、逆问题的贝叶斯方法

Part I : 从极大似然到生成模型

1.课程概览

2. MLE/MAP方法回顾

基本框架设定

判别模型与生成模型的区分

2.1 极大似然估计的定义

MLE的核心思想

可识别性的前提条件

从MLE到MAP

2.2 MLE与KL散度的联系

KL散度的定义与性质

MLE如何最小化KL散度

例2.1：多元高斯观测的MLE

问题设定

构建对数似然函数

求解均值的MLE

求解协方差矩阵的MLE

MLE闭式解的局限性

例2.2：学习固定模型的最优混合

混合模型的问题设定

混合模型MLE的困难

高斯混合模型的现实挑战

隐变量的启发

EM算法的动机

2.3 通过随机梯度下降求解MLE

梯度下降法的基本思路

随机梯度下降的优势

小批量梯度估计的构造

SGD算法流程

3 EM算法

3.1 EM算法的定义与性质

3.1.1 EM算法的基本定义

完全数据的概念

EM算法的两步迭代结构

E步骤的深层含义

M步骤的优化目标

迭代收敛的过程

EM算法的收敛性质

3.1.2 完全数据的性质要求

完全数据的可优化性

马尔可夫关系的必要性

确定性函数的特例

3.1.3 独立同分布样本情况下的EM形式

i.i.d.假设下的问题简化

引理3.1：i.i.d.情况下Q函数的分解

证明过程：建立样本的相互独立性

证明过程：推导后验分布的分解

完成Q函数分解的证明

3.2 EM算法的收敛性分析

EM算法的单调性保证

定理3.2：EM算法的收敛性

定理的前提条件

定理的结论

证明的整体思路

引入辅助分布进行重写

应用Jensen不等式

利用贝叶斯定理和马尔可夫性质

分解对数并展开期望

推导单调性结论

注释3.3：单调性的深层含义

单调性与收敛性的区别

有界性对收敛性的作用

注释3.4：在EM算法中引入先验

为何需要引入先验

MAP-EM算法的形式

MAP-EM的单调性

3.3 通过EM算法学习隐变量模型

隐变量模型的基本框架

隐变量模型中的完全数据定义

深度隐变量模型的表达能力

3.4 学习固定模型的最优混合

问题的具体设定

完全数据的选择与直觉

E步：计算后验概率

M步：推导更新公式

使用拉格朗日乘数法处理约束

EM算法的完整流程

4 超越标准EM算法