8.5 最大期望（EM）算法

最大期望（expectation–maximization，EM） 算法是简化复杂的最大似然问题的常见工具。本节先以一个简单的混合模型为例介绍这个算法。

8.5.1 二成分混合模型

本节介绍一个简单混合模型的密度估计，并用最大期望算法来实现最大似然估计。这与贝叶斯推断中的吉布斯采样（Gibbs sampling）方法有自然的联系。在本书中很多次用混合模型进行演示，例如第 6.8 节、第 12.7 节、第 13.2.3 节。

表 8.1：图 8.5 中所用的二成分混合模型的 20 个模拟数据。

图 8.5 的左图展示了表 8.1 中列出的 20 个模拟数据的直方图。

我们想要为数据的密度函数建模，由于直方图中明显呈现双峰，可能用单独的高斯分布是不恰当的。看起来似乎有潜在生成机制，所以用两个正态分布的混合模型来描述 $Y$：

其中 $\Delta\in\{0,1\}$，并且 $\operatorname{Pr}(\Delta=1)=\pi$。这是一个显式的生成表达式：先根据概率 $\pi$ 生成 $\Delta\in\{0,1\}$，然后再根据其结果最终给出 $Y_1$ 和 $Y_2$ 的其中一个。令 $\phi_\theta(x)$ 表示参数为 $\theta=(\mu,\sigma^2)$ 的正态分布密度函数。则 $Y$ 的密度函数为：

现若要在图 8.5 的数据上用最大似然来拟合这个模型。参数为：

基于 $N$ 个训练样本的对数似然函数为：

由于对数中存在求和，直接对 $\ell(\theta;\mathbf{Z})$ 求最大化在数值计算上非常困难。不过存在一个简化的方法。式 8.36 中不可观测的隐变量 $Delta_i$ 取值为 0 或 1：若 $\Delta_i=1$ 则由模型 2 产生 $Y_i$，否则由模型 1 产生 $Y_i$。假设已知 $\Delta_i$ 的值，则对数似然函数为：

那么 $\mu_1$ 和 $\sigma_1^2$ 的最大似然估计就是 $\Delta_i= 0$ 的那部分数据的样本均值和方差，与之类似，$\mu_1$ 和 $\sigma_1^2$ 的估计就是 $\Delta_i=1$ 的那部分数据的样本均值和方差。而 $\pi$ 的估计便是样本中 $\Delta_i=1$ 的比例。

由于实际上并不知道 $\Delta_i$ 的值，可采取迭代的方式，用其期望值代替式 8.40 中的每个 $\Delta_i$：

这也被称为模型 2 对观测值 $i$ 的 责任值（responsibility）。算法 8.1 列出了在这个高斯混合模型特例中的最大期望（EM）算法的过程。在期望（E）步骤将每个观测“软分配”（soft assignment）给每个模型，即基于当前的参数估计值按照训练样本点在每个模型下的相对密度计算/分配责任值。在最大化（M）步骤根据这些责任值计算加权最大似然拟合，从而更新参数的估计。

一个简单而有效的构建 $\mu_1$ 和 $\mu_2$ 的初始猜测值的方法是随机从样本中抽取两个 $y_i$。初始的 $\sigma_1^2$ 和 $\sigma_2^2$ 可设置为整体样本的方差 $\sum_{i=1}^N(y_i-\bar{y})^2/N$。混合比例 $\hat{\pi}$ 可初始化为 0.5。

算法 8.1：二成分高斯混合模型的最大期望（EM）算法

1. 确定参数的初始猜测值 $\hat{\mu}_1$、$\hat{\sigma}_1^2$、$\hat{\mu}_2$、$\hat{\sigma}_2^2$、$\hat{\pi}$（见正文描述）。
2. 期望（E）步骤：计算责任值 $$\hat{\gamma}_i = \frac{\hat{\pi}\phi_{\hat{\theta}_2}(y_i)} {(1-\hat{\pi})\phi_{\hat{\theta}_1}(y_i) + \hat{\pi}\phi_{\hat{\theta}_2}(y_i)}, i=1,2,\dots,N \tag{8.42}$$
3. 最大化（M）步骤：计算加权均值和方差 $$\begin{matrix} \hat{\mu}_1 = \frac {\sum_1^N (1-\hat{\gamma}_i)y_i} {\sum_{i=1}^N (1-\hat{\gamma}_i)} & \hat{\sigma}_1^2 = \frac {\sum_{i=1}^N (1-\hat{\gamma}_i)(y_i-\hat{\mu}_1)^2} {\sum_{i=1}^N (1-\hat{\gamma}_i)} \\ \hat{\mu}_2 = \frac {\sum_{i=1}^N \hat{\gamma}_i y_i} {\sum_{i=1}^N \hat{\gamma}_i} & \hat{\sigma}_2^2 = \frac {\sum_{i=1}^N \hat{\gamma}_i (y_i-\hat{\mu}_2)^2} {\sum_{i=1}^N \hat{\gamma}_i} \end{matrix}$$ 以及混合概率：$\hat{\pi}=\sum_{i=1}^N\hat{\gamma}_i/N$。
4. 重复步骤 2 和 3，直到结果收敛。

注意当将所有权重赋予任一数据点 $i$，相应地 $\hat{\mu_1}=y_i$ 和 $\hat{\sigma}_1^2=0$，这实际上会得到似然函数的最大值解。这时的似然函数是无穷大，然而这个解并不可用。因此实际上需要找到似然函数的局部最大值解，满足 $\hat{\sigma}_1^2,\hat{\sigma}_2^2>0$。在更复杂的场景中，可能存在多个局部最大值解，并且都满足 $\hat{\sigma}_1^2,\hat{\sigma}_2^2>0$。在这个例子中，我们在最大期望算法中使用了很多不同的参数初始猜测值，使其全都满足 $\hat{\sigma}_k^2>0.5$，然后从中选择得出了最大的对数函数最大值的一组。图 8.6 展示了最大期望算法的最大化对数似然函数过程。表 8.2 展示了在选定的某几次算法迭代过程中的 $\hat{\pi}=\sum_i\hat{\gamma_i}/N$，即观测样本中类别 2 的比例的最大似然估计值。

表 8.2：混合模型示例的最大期望（EM）算法过程中选定的某几次迭代。

最终得到最大似然估计为：

图 8.5 的右图展示了这个过程对高斯混合密度函数的估计（红色实线），以及相应的责任值（绿色虚线）。另外，混合模型也可被用在监督学习中；第 6.7 节中展示了如何从高斯混合模型得到一组径向基函数。

8.5.2 一般场景下的最大期望算法 😱

上述过程是一个应用最大期望（EM）（或鲍姆-韦尔奇，Baum-Welch）算法解决某个特定类型问题中的最大似然估计的例子。在这些问题中，似然函数的最大化往往比较困难，但可以通过引入（无法观测的）隐变量来扩大样本维度，将问题简化。这种方法称为 数据增强（data augmentation）。这里的隐变量为模型的成分指示变量 $\Delta_i$。在其他场景中，隐变量为本应该被观测到但却缺失的真实数据。

算法 8.2：最大期望（EM）算法

1. 确定参数的初始猜测值 $\hat{\theta}^{(0)}$。
2. 期望（E）步骤：在第 j 步，计算虚（dummy）参数 $\theta’$ 的函数： $$Q(\theta', \hat{\theta}^{(j)}) = \operatorname{E}(\ell_0(\theta'; \mathbf{T})|\mathbf{Z}, \hat{\theta}^{(j)}) \tag{8.43}$$
3. 最大化（M）步骤：将最大化 $Q(\theta’,\hat{\theta}^{(j)})$ 的 $\theta’$ 解作为新的估计值 $\hat{\theta}^{(j+1)}$。
4. 重复步骤 2 和 3，直到结果收敛。

算法 8.2 给出了一般情况下的最大期望（EM）算法过程。观测数据是 $\mathbf{Z}$，对数似然函数 $\ell(\theta;\mathbf{Z})$ 同时依赖于参数 $\theta$。隐变量或缺失数据为 $\mathbf{Z}^m$，因此全部数据为 $\mathbf{T}=(\mathbf{Z},\mathbf{Z}^m)$，对数似然函数为 $\ell_0(\theta;\mathbf{T})$，函数 $\ell_0$ 基于整体的密度。在混合模型问题中，$(\mathbf{Z},\mathbf{Z}^m)=(\mathbf{y},\mathbf{\Delta})$，而 $\ell_0(\theta;\mathbf{T})$ 即为等式 8.40。

在混合模型例子中，$\operatorname{E}(\ell_0(\theta’;\mathbf{T})|\mathbf{Z},\hat{\theta}^{(j)})$ 即为用责任值 $\hat{\gamma}_i(\hat{\theta})$ 代替了 $\Delta_i$ 的等式 8.40，并且第三步中的最大化解即为加权均值和方差。

下面将说明在一般的场景中最大（EM）算法有效的原因。

由于：

所以可得：

从对数似然函数的角度，$\ell(\theta’;\mathbf{Z})=\ell_0(\theta’;\mathbf{T})-\ell_1(\theta’;\mathbf{Z}^m|\mathbf{Z})$，其中 $\ell_1$ 基于条件密度函数 $\operatorname{Pr}(\mathbf{Z}^m|\mathbf{Z},\theta’)$。对其取关于参数为 $\theta$ 的 $\mathbf{T}|\mathbf{Z}$ 条件分布的期望，得到：

在第 M 步，最大期望（EM）算法选择使 $Q(\theta’,\theta)$ 最大的 $\theta’$，而不是使实际的目标函数 $\ell(\theta’;\mathbf{Z})$ 最大。那为何它会成功地最大化 $\ell(\theta’;\mathbf{Z})$？注意到 $R(\theta^*,\theta)$ 为由参数 $\theta^*$ 索引的密度函数的对数似然函数的期望，而取期望的密度函数是当参数为 $\theta$ 时的该密度函数；因此根据延森（Jensen）不等式，作为 $\theta^*$ 的函数，$R$ 在 $\theta^*=\theta$ 处取到最大值（见练习 8.1）。所以若 $\theta’$ 最大化 $Q(\theta’,\theta)$，则有：

因此最大期望（EM）的迭代永远不会降低对数似然函数。

从这个论证中也可见，在第 M 步中并不需要完全的最大化：只需要找到可以使 $Q(\theta’,\hat{\theta}^{(j)})$ 是其第一个参数的递增函数的值 $\hat{\theta}^{(j+1)}$，即满足 $Q(\hat{\theta}^{(j+1)},\hat{\theta}^{(j)})>Q(\hat{\theta}^{(j)},\hat{\theta}^{(j)})$。这样的过程也称为 广义最大期望（generalized EM，GEM） 算法。最大期望（EM）算法也可以被视为一个少数化（minorization）过程，见练习 8.7。

8.5.3 EM 是一个最大化-最大化过程 😱

从另一个角度可将最大期望（EM）过程视为一个联合最大化算法。考虑函数：

其中 $\tilde{P}(\mathbf{Z}^m)$ 为在混合模型例子中，$\tilde{P}(\mathbf{Z}^m)$ 即为概率 $\gamma_i=\operatorname{Pr}(\Delta_i=1|\theta,\mathbf{Z})$ 的集合。注意到在 $\tilde{P}(\mathbf{Z}^m)=\operatorname{Pr}(\mathbf{Z}^m|\mathbf{Z},\theta’)$ 处时取值的 $F$ 就是观测数据的对数似然函数，如等式 8.46¹。函数 $F$ 拓展了对数似然函数的范畴，并提供了对其最大化的简便方法。

最大期望（EM）算法可视为函数 $F$ 对 $\theta’$ 和 $\tilde{P}(\mathbf{Z}^m)$ 的联合最大化方法，固定住其中一个对另一个进行最大化。可证明（练习 8.2）对固定的 $\theta’$，$\tilde{P}(\mathbf{Z}^m)$ 的最大化解为：

这是在期望（E）步骤中计算的分布，例如在混合模型例子中的等式 8.42。在最大化（M）步骤中，固定 $\tilde{P}$ 后求使函数 $F(\theta’,\tilde{P})$ 最大化的 $\theta’$：由于第二项中没有出现 $\theta’$，这等同于对第一项 $\operatorname{E}_\tilde{P}[\ell_0(\theta’;\mathbf{T}|\mathbf{Z},\theta)]$ 最大化。

最后，由于当 $\tilde{P}(\mathbf{Z}^m)=\operatorname{Pr}(\mathbf{Z}^m|\mathbf{Z},\theta’)$ 时，$F(\theta’,\tilde{P})$ 与观测数据的对数似然函数相等，因此对前者的最大化也同时会对后者最大化。图 8.7 展示了这个过程的示意图。从中可看出最大期望（EM）算法引导了另一条最大化路径。例如，并不需要同时对所有隐变量的参数最大化，而可以在 M 步骤分别交替地最大化。

本节练习

练习 8.1

Let $r(y)$ and $q(y)$ be probability density functions. Jensen’s inequality states that for a random variable X and a convex function $\phi(x)$, $\operatorname{E}[\phi(X)]geq\phi(\operatorname{E}[X])$$. Use Jensen’s inequality to show that

is maximized as a function of $r(y)$ when $r(y)=q(y)$. Hence show that $R(\theta,\theta)\geq R(\theta’,\theta)$ as stated below equation (8.46).

练习 8.2

Consider the maximization of the log-likelihood (8.48), over distributions $\tilde{\operatorname{P}}(\mathbf{Z}^m)$ such that $\tilde{\operatorname{P}}(\mathbf{Z}^m)\geq0$ and $\sum_{\mathbf{Z}^m}\tilde{\operatorname{P}}(\mathbf{Z}^m)=1$. Use Lagrange multipliers to show that the solution is the conditional distribution $\tilde{\operatorname{P}}(\mathbf{Z}^m)=\operatorname{Pr}(\mathbf{Z}^m|\mathbf{Z},\theta’)$, as in (8.49).

练习 8.7

EM as a minorization algorithm (Hunter and Lange, 2004; Wu and Lange, 2007).

A function $g(x,y)$ to said to minorize a function $f(x)$ if

for all x, y in the domain. This is useful for maximizing $f(x)$ since it is easy to show that $f(x)$ is non-decreasing under the update

There are analogous definitions for majorization, for minimizing a function $f(x)$. The resulting algorithms are known as MM algorithms, for “Minorize- Maximize” or “Majorize-Minimize.”

Show that the EM algorithm (Section 8.5.2) is an example of an MM algorithm, using $Q(\theta’,\theta)+\log\operatorname{Pr}(\mathbf{Z}|\theta)-Q(\theta,\theta)$ to minorize the observed data log-likelihood $\ell(\theta’;\mathbf{Z})$. (Note that only the first term involves the relevant parameter $\theta’$).