最大期望算法（EM）

1. 数学基础

1.1 最大似然估计

当我们从一个总体中拿到了 $n$ 个样本的时候，我们可以认为，我们之所以取到了它们，是因为它们出现的概率比较大． 并且，我们已知它的概率分布模型，但是不知道其中的一个参数 $\theta$．这时候我们可以把 $\theta$ 设为变量，因为样本显然是独立同分布的（这跟废话一样，从同一个总体中取出的样本，肯定是同分布的），所以取到这个样本序列的概率理所当然地等于分别取到每一个值的概率之积，即

$$P = \prod_{i=1}^{n}f(x_i;\theta)$$

其中$f(x_i;\theta)$ 表示当参数取值为 $\theta$ 的时候，模型的概率密度函数在 $x_i$ 处的取值．

当我们把样本值$(x_1, x_2,\cdots, x_n)$代入进去，那么 $P$ 就是关于 $\theta$ 的一元函数． 但既然是函数，我们就要用函数的方法表示它．我们将这个函数记为 $L(\theta)$，叫做样本的 似然函数．

刚才我们提到，我们认为取到这个样本序列的原因是它们出现的概率比较大.那么我们就干脆让上面这个样本序列的概率，也就是似然函数，达到最大值. 这个最大值的意义是，未知参数 $\theta$ 的取值要能够让取到这个样本序列的概率尽可能地大．

通过求一个函数导数的零点来求它的极值，这是坠吼的！所以我们对$L(\theta)$求导，令$\frac{dL(\theta)}{d\theta}=0$即可．但是前面也说过了，$L(\theta)$是个连乘函数，一弄不好就搞个 $n$ 次方，实在是太坑爹了．如果能把乘法变成加法，那就好得多了．那么如何把乘法变加法呢？没错，取对数．然而问题来了，取了对数，又怎么保证取了对数的函数和原函数有相同的极值点呢？这个是可以推出来的：由复合函数的求导公式可得$\frac{d\ln L(\theta)}{d\theta} = \frac{d\ln \theta}{d\theta} \frac{dL(\theta)}{d\theta} = \frac{1}{\theta} \frac{dL(\theta)}{d\theta}$， 因此 $\ln L(\theta)$ 和 $L(\theta)$ 具有相同的增减性. 因此，问题转换为求方程 $\frac{d\ln L(\theta)}{d\theta} = 0$ 的解．解得的 $\theta$ 就是要求的参数. 求得的估计量 $\hat{\theta}$ 叫做最大似然估计量．

当参数有多个时，比如正态分布有两个参数 $\mu$ 和 $\sigma$，我们就把导数变成偏导数，方程变为方程组．当有 $n$ 个参数$\theta_1, \theta_2,\cdots,\theta_n$ 的时候，方程组变为

$$\begin{cases} \cfrac{\partial \ln L(\boldsymbol{\theta})}{\partial \theta_1} &= 0 \\ \cfrac{\partial \ln L(\boldsymbol{\theta})}{\partial \theta_2} &= 0 \\ &\vdots\\ \cfrac{\partial \ln L(\boldsymbol{\theta})}{\partial \theta_n} &= 0 \\ \end{cases}$$

解这个方程组即可．

1.2 琴生不等式

琴生不等式属于凸优化的领域．我们介绍一下它的内容：

如果一个函数$f(x)$是凸函数，那么

$$\forall \lambda_1, \lambda_2, ... \lambda_n > 0,\\ \sum_{i=1}^{n} \lambda_i = 1 \implies f(\sum_{i=1}^{n} \lambda_i x_i) \le \sum_{i=1}^{n} f(\lambda_i x_i).$$

如果是在概率论领域，那么可以把$X$和$f(x)$看做离散型随机变量,$x_i$和$f(x_i)$是可以取到的值，而$\lambda_i$是取到$x_i$或$f(x_i)$的概率，那么上面的不等式就写成

$$f[E(X)] \le E[f(x)].$$

从几何角度来理解，不等式的含义是，过凸函数任意两点的线段，都包含在该凸函数的上境图中．（上境图是指函数曲线正上方的点集）

下面我们分别从数学分析和概率论的角度来证明这个不等式．

数学分析方法

如果f(x)二阶可导，那么先令$x_0 = \sum_{i=1}^{n}\lambda_i x_i$，在$x_0$点进行泰勒展开可得

$$\begin{aligned}f(x) &= \frac{1}{0!}f(x_0) + \frac{1}{1!}f'(x_0)+\frac{1}{2!}f''(\xi)(x-x_0)^2\\ &= f(x_0) + f'(x_0)(x-x_0) + \frac{1}{2!}f''(\xi)(x-x_0)^2 \end{aligned}$$

其中$\frac{1}{2!}f''(\xi)(x-x_0)^2$是拉格朗日余项，$\xi$ 在 $x$ 和 $x_0$ 之间．

凸函数的二阶导数大于零，所以 $\frac{1}{2!}f''(\xi)(x-x_0)^2 \ge 0$，因此有

$$\begin{aligned} f(x) &\ge f(x_0) + f'(x_0)(x-x_0)\\ f(x_i) &\ge f(x_0) + f'(x_0)(x_i-x_0)\\ \lambda_i f(x_i) &\ge \lambda_i f(x_0) + \lambda_i f'(x_0)(x_i-x_0)\\ \sum_{i=1}^{n} \lambda_i f(x_i) &\ge \sum_{i=1}^{n}\lambda_i f(x_0) + \sum_{i=1}^{n}\lambda_i f'(x_0)(x_i-x_0)\\ \begin{matrix}(\text{注意}\sum_{i=1}^{n}\lambda_i = 1)& & \end{matrix} & =f(x_0) + f'(x_0)(\sum_{i=1}^{n}\lambda_i x_i - \sum_{i=1}^{n} \lambda_i x_0)\\ &=f(x_0) + f'(x_0)(\sum_{i=1}^{n}\lambda_i x_i - x_0)\\ \begin{matrix}\text{先前定义了}(x_0 = \sum_{i=1}^{n}\lambda_i x_i) & & \end{matrix} & = f(x_0) = f(\sum_{i=1}^{n}\lambda_i x_i) \end{aligned}$$

证毕．

概率论方法

凸函数有一条性质：过凸函数一点的切线不在凸函数图象的上方．

已知$X$的期望$E(X) = \sum_{i=1}^{n}\lambda_i x_i$，过$(E(X),f[E(X)])$做函数切线$\ell(x)=ax+b$，易得

$$\begin{aligned} f(x)\ge \ell(x) \iff E[f(X)] &\ge E[\ell(X)] = aE(X) + b = \ell[E(x)] = f[E(X)]. \end{aligned}$$

证毕．

2. 最大期望算法

2.1 引入

以三硬币模型为例，介绍最大期望算法的使用场景.

我们指定一个游戏规则：$A,B,C$是三枚硬币，正面朝上的概率各不相同($\pi, p, q$).如果抛$A$的结果是正面，则抛$B$；否则抛$C$. 正面记作 1，反面记作0. 用数学语言描述，是这样的：

$$A\sim B(n,\pi), B\sim B(n, p), C\sim B(n,q).$$

根据以上的游戏规则，我们定义最终抛出来的正反面值随机变量为 $Y$，抛出硬币 $A$ 的正反面值随机变量是 $Z$，这个游戏的参数$\theta=(\pi,p,q)$.

现在我们添加一个条件：$Z$的值和参数$\theta$都是不可见的，这好比另一个人在另一个屋子里抛硬币，每次只告诉你这次抛出 $0$ 还是 $1$，却不告诉你他抛的是 $B$ 还是 $C$.

这是我们生活中经常遇到的问题类型. 我们知道未知参数$\pi$的分布模型（二项分布），但是我们不知道它的参数. 如果我们能拿到数据的话，就可以用最大似然估计来估计参数，但是我们连数据也没有，只有一些不由这个变量所唯一决定的数据（在这个模型中，就是我们最终抛出的结果）. 另外，有些时候我们还需要知道 $p$ 和 $q$，但这些问题都可以归为一类.

现在我们想根据那个人给出的试验数据，来估计变量$Z$满足的分布模型（二项分布）中未知参数 $\pi,p,q$ ，也就是 $\theta$. 这是一个很复杂的过程，我们先从头开始分析.

首先，令 $Y$ 的一个可以取到的值为 $y$，则有

$$\begin{aligned} P(y;\theta) &= \sum_z P(y,z; \theta) = \sum_z P(z;\theta)P(y|z;\theta)\\ &= \pi p^y (1-p)^{(1-y)} + (1-\pi) q^y (1-q)^{(1-y)} \end{aligned}$$

解释一下这个公式. 因为我们不知道这一次试验的 $z$ 是多少，所以对于每一种可能取到的 $z$ 值，都要计算“取到 $y$ ”和“取到这个 $z$ ”两个事件同时发生的概率. 考虑到所有可能的 $z$ 的取值之后，再将这些计算出的概率求和.

从一个数据点推广到数据序列，也是一样的.

令观测数据$Y=(Y_1,Y_2,\cdots,Y_n)^T$，隐藏的数据$Z=(Z_1,Z_2,\cdots,Z_n)^T$我们得到公式

$$P(Y;\theta) = \sum_Z P(Z;\theta)P(Y|Z;\theta).$$

显然似然函数

$$L(\theta)=P(Y;\theta)=\prod_{j=1}^{n}\bigl[ \pi p^{y_j} (1-p)^{1-y_j}+(1-\pi) q^{y_j} (1-q)^{1-y_j} \bigr].$$

这个函数是没有解析解的．那么我们该怎么解决这个问题呢？此时就应该请出 EM 算法了．

世界上没有十全十美的事情．即使是EM算法，也算不出它的解析解．但是它可以通过迭代的方式使得解逐渐趋于局部最优．

2.2 算法需求及推导

EM，EM，顾名思义，一个 E(Expectation，期望) 一个 M(Maximization，最大化)．算法的迭代部分也主要由这两步组成．

上面的方程里，如果$\theta$是已知的，那么求$Z$的期望就很容易；如果$Z$是已知的，那么用最大似然估计求$\theta$也很简单．所以EM算法的思想是：先选定一个初始的$\theta_0$,用它来求$Z$的期望，把它当作$Z$；然后用这个$Z$放到公式里求$\theta$的最大似然估计，得到新$\theta$．周而复始，直到结果收敛为止．

还要说明一点，很多时候，我们不求$Z$的期望，而是去求$Z$基于$\theta$的概率分布$P(Z|Y;\theta)$的期望．目的都是为了用它求$\theta$的最大似然估计．

下面我们主要说说为什么要用 $P(Z|Y;\theta)$ 的期望来估计 $\theta$．

推导

现在我们把问题抽象一下：我们要求似然函数

$$\begin{aligned} L(\theta) &= \log(Y|\theta) = \log\sum_Z P(Y,Z;\theta) \end{aligned}$$

的极大值，来进行极大似然估计．

要通过迭代的方法，使得新的$\theta$能让$L(\theta)>L(\theta^{(i)})$，因此要求 $\arg \max_\theta \left(L(\theta) - L(\theta^{(i)}\right)$．

我们遇到的主要困难是上式中的隐变量$Z$，以及给求导造成困难的和（或积分）的对数．

先写出

$$\begin{aligned} \left(L(\theta) - L(\theta^{(i)}\right) &= \log \sum_Z {\rm P}(Y,Z;\theta) - \log {\rm P}(Y;\theta)\\ \begin{matrix} \text{乘}{\rm P}(Z|Y;\theta^{(i)})\text{，变形．}\\ \text{注意到}\sum_ZP(Z|Y;\theta^{(i)})=1\text{．}& & \end{matrix} &=\log\left[\sum_ZP(Z|Y;\theta^{(i)})\frac{{\rm P}(Y,Z;\theta)}{{\rm P}(Z|Y;\theta^{(i)})}\right] - \log{\rm P}(Y;\theta)\\ \begin{matrix}\text{应用琴生不等式}& & \end{matrix} &\ge \sum_Z P(Y|Z;\theta^{(i)})\log\frac{{\rm P}(Y,Z;\theta)}{{\rm P}(Z|Y;\theta^{(i)})} - \log{\rm P}(Y;\theta)\\ &= \sum_Z P(Y|Z;\theta^{(i)})\log\frac{{\rm P}(Y,Z;\theta)}{{\rm P}(Z|Y;\theta^{(i)}){\rm P}(Y;\theta^{(i)})} \end{aligned}$$

所以我们找到了 $L(\theta)$ 的下界．令下界

$$B(\theta,\theta^{(i)}) \hat{=} L(\theta^{(i)}) + \sum_Z P(Y|Z;\theta^{(i)})\log\frac{{\rm P}(Y,Z;\theta)}{{\rm P}(Z|Y;\theta^{(i)}){\rm P}(Y;\theta^{(i)})}.$$

另外把$\theta=\theta^{(i)}$带入$B(\theta,\theta^{(i)})$得到$L(\theta)=B(\theta,\theta^{(i)})$．

所以对于任意的$\theta$满足$B(\theta,\theta^{(i)})$，都有$L(\theta) >L(\theta^{(i)})$．也就是说，若$\theta$可以使$B(\theta,\theta^{(i)})$增大，则$\theta$也可以使 $L(\theta)$增大．注意$\theta$可使$B(\theta,\theta^{(i)})$减小时，$L(\theta)$不一定减小．

为了使$L(\theta)$尽可能大，我们可以使$B(\theta,\theta^{(i)})$达到极大．求$\frac{\partial B}{\partial \theta}$时，$\theta^{(i)}$可看作常数而被省去．故有

$$\begin{aligned} \theta^{(i+1)} &= \arg\max_\theta B(\theta,\theta^{(i)})\\ &= \arg\max_\theta \left( L(\theta^{(i)}) + \sum_Z P(Y|Z;\theta^{(i)})\log\frac{{\rm P}(Y,Z;\theta)}{{\rm P}(Z|Y;\theta^{(i)}){\rm P}(Y;\theta^{(i)})} \right)\\ &= \arg\max_\theta \left( \sum_Z P(Z|Y;\theta^{(i)})\log P \left(Y,Z;\theta \right) \right)\\ &= \arg\max_\theta E_Z\left[\log P(Y,Z;\theta)|Y;\theta\right] \end{aligned}$$

令

$$Q(\theta,\theta^{(i)}) = E_Z\left[\log P(Y,Z;\theta)|Y;\theta\right]$$

则我们就求 $\arg \max_\theta Q(\theta,\theta^{(i)})$，即可得到新的 $\theta$．

在该公式里面的$Q(\theta,\theta^{(i)})$函数是这样的形式：

$$Q(\theta^{(i)},\theta^{(i+1)}) = E_Z\bigl[\log P(Y,Z;\theta)|Y,\theta^{(i)}\bigr]$$

它的意思是：在给定观测数据$Y$和当前参数$\theta^{(i)}$的情况下，$\log P(Y,Z;\theta)$的期望．

至于公式的展开，首先让我们复习一下离散随机变量的期望公式：$E(X)=\sum_{x\in X} x\cdot p(x)$．所以

$$\sum_Z \log \left(P(Y,Z;\theta)P(Z|Y;\theta^{(i)})\right) = E_Z\left[\log P(Y,Z;\theta)|Y;\theta\right]$$

也就是理所当然了．

总之，求$Q(\theta, \theta^{(i)})$就是 E 步骤，求$\arg \max_\theta Q(\theta,\theta^{(i)})$就是 M 步骤．

2.3 算法描述

当结果不收敛时，循环以下两步

{

（当前步骤计数为$i$,从$0$开始）

E 步骤： 以 $\theta^{(i)}$ 和 $Y$ 计算

$$\begin{aligned} Q(\theta,\theta^{(i)}) &= E_Z\bigl[\log P(Y,Z;\theta)|Y,\theta^{(i)}\bigr]\\ &= \sum_Z \log \left(P(Y,Z;\theta)P(Z|Y;\theta^{(i)})\right). \end{aligned}$$

其中$P(Z|Y;\theta^{(i)})$是给定观测数据 $Y$ 和 参数估计 $\theta^{(i)}$ 下，隐变量 $Z$ 的条件概率分布；$Q$函数的第一个变量是要求的参数，第二个变量是我们上一步估计的参数．

M 步骤： 求$Q(\theta,\theta^{(i)})$关于参数$\theta$的最大似然估计，即

$$\theta^{(i+1)}=\arg \max_\theta Q(\theta,\theta^{(i)}).$$

}

停止迭代的条件是，对较小的整数$\varepsilon_1, \varepsilon_2$，满足 $\lVert \theta^{(i+1)}-\theta^{(i)}\rVert<\varepsilon_1$ 或 $\lVert Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta{^(i)})\rVert<\varepsilon_2$，则停止迭代．

2.4 算法分析

根据推导过程，我们只能保证当$B(\theta,\theta^{(i)})$增加时， $L(\theta)$ 一定增加．但是当 $B(\theta,\theta^{(i)})$ 达到极值，开始减小时，$L(\theta)$不一定随之减小．换言之，二者同增不同减．也就是说，我们有可能忽略掉$L(\theta)$的全局最优解，只选择了它的局部最优解．但其实这也是没办法的事情．如果可以在不穷举的情况下找到它的全局最优解，那也就证明了命题 ${\rm P=NP}$．然而人类到现在也没能证明或证伪这个命题．在数学家们摘取硕果之前，我们只能退而求其次了．