前言

$EM$算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。

$EM$算法的每次迭代由两步组成：

• $E$步：求期望
• $M$步：求极大

$EM$算法引入

首先介绍一个使用$EM$算法的例子：三硬币模型

假设有3枚硬币，分别记作$A,B,C$，这些硬币正面出现的概率分别是$\pi,p,q$，进行如下试验：
先掷硬币$A$根据其结果选出硬币$B$或者硬币$C$，正面选硬币$B$，反面选硬币$C$；然后掷出选出的硬币，出现正面记为1，反面记为0；独立地重复$n$次试验（这里$n$取10），观测结果如下：
1101001101
假设只能观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三个硬币正面出现的概率，即三硬币模型的参数。

三硬币模型可写为：

$$p(x|\theta) = \sum_z p(x,z|\theta) = \sum_z p(z|\theta)p(x|z, \theta) \\ = \pi p^x(1-p)^{1-x}+ (1-\pi)q^x(1-q)^{1-x}$$
这里，随机变量 $x$是观测变量，表示一次试验观测的结果1或0；随机变量 $z$是隐变量，表示未观测到的掷硬币 $A$的结果； $\theta=(\pi,p,q)$是模型参数。

下面将观测数据表示为$X = (x_1,x_2,...,x_n)^T$，未观测数据表示为$Z = (z_1,z_2,...,z_n)^T$，则观测数据的似然函数为：

$$p(X|\theta) = \sum_z p(Z|\theta)p(X|Z,\theta)$$
即：
$$p(X|\theta) = \prod_{j=1}^{n}[\pi p^{x_j}(1-p)^{1-x_j}+ (1-\pi)q^{x_j}(1-q)^{1-x_j}]$$
考虑求模型参数 $\theta=(\pi,p,q)$的极大似然估计，取对数似然：

$$\theta_{MLE} = arg \max_\theta log \ p(X|\theta)$$

这个问题没有解析，只有通过迭代的方法来求解。$EM$算法就是可以用于求解这个问题的一种迭代算法。

$EM$算法通过迭代求$\mathcal{L}(\theta) = log \ p(X|\theta)$的极大似然估计。每次迭代包含两步：
E步：求期望；M步：求极大化。

$EM$算法

输入： 观测变量数据$X$，隐变量数据$Z$，联合分布$p(X,Z|\theta)$，条件分布$p(Z|X,\theta)$;

输出：模型参数$\theta$

（1）选择参数的初值$\theta^{(0)}$，开始迭代；
（2）E步：记$\theta^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的E步，计算：

$$Q(\theta,\theta^{(i)} )= E_z[log \ p(X,Z|\theta)|X,\theta^{(i)}] \\ =\sum_z p(X,Z|\theta)p(Z|X,\theta^{(i)})$$

这里，$p(Z|X,\theta^{(i)})$表示在给定观测数据$X$和当前的参数估计$\theta^{(i)}$下隐变量数据$Z$的条件概率分布；
（3）M步：求使$Q(\theta,\theta^{(i)} )$极大化的$\theta$，确定第$i+1$次迭代的参数估计值$\theta^{(i+1)}$

$$\theta^{(i+1)} = arg \max_\theta Q(\theta,\theta^{(i)} )$$
（4）重复第（2）步和第（3）步，直到收敛。

注：
函数$Q(\theta,\theta^{(i)} )$是$EM$算法的核心，称为Q函数。
参数的初值可以任意选择，但$EM$算法对初值是敏感的
$Q(\theta,\theta^{(i)} )$中$\theta$表示要极大化的参数，$\theta^{(i)}$表示当前估计值，每次迭代实际在求Q函数及其极大化
M步求Q函数的极大化，得到$\theta^{(i+1)}$，完成一次迭代
第（4）步给出停止迭代的条件，一般是对较小的正数$\mathcal{E}_1,\mathcal{E}_2$，若满足：

$$||\theta^{(i+1)}-\theta^{(i)}|| \leqslant \mathcal{E}_1,或 ||Q(\theta^{(i+1)},\theta^{(i)} )-Q(\theta^{(i)} ,\theta^{(i)} )|| \leqslant \mathcal{E}_2$$
则停止迭代。

导出$EM$算法

上面我们给出了$EM$算法，那么为什么$EM$算法能近似实现对观测数据的极大似然估计呢？下面通过近似求解观测数据的对数似然函数的极大化问题来导出$EM$算法。

由上面内容可以知道，极大化观测数据$X$关于参数$\theta$的对数似然函数，即极大化：

$$\mathcal{L}(\theta) = log\ p(X|\theta) = log \ \sum_z p(X,Z|\theta)= log \big(\sum_z p(X|Z,\theta)p(Z|\theta)\big)$$

注意到：上式中有未观测数据并包含和的对数

事实上，$EM$算法通过迭代逐步近似极大化$\mathcal{L}(\theta)$的。

假设在第$i$次迭代后$\theta$的估计值是$\theta^{(i)}$，我们希望新估计值$\theta$能使$\mathcal{L}(\theta)$增加，即$\mathcal{L}(\theta) > \mathcal{L}(\theta^{(i)})$，并逐步达到极大值，为此考虑两者的差：

$$\mathcal{L}(\theta) - \mathcal{L}(\theta^{(i)}) = log \big(\sum_z p(X|Z,\theta)p(Z|\theta)\big) - log\ p(X|\theta^{(i)})$$

利用$Jensen$不等式，得到其下界：

$$log \sum_j \lambda_jy_i \geqslant \sum_j \lambda_j log\ y_j$$
其中 $\lambda_j \geqslant 0, \sum_j \lambda_j = 1$

$$\mathcal{L}(\theta) - \mathcal{L}(\theta^{(i)}) = log \big(\sum_z p(X|Z,\theta^{(i)})\frac{p(X|Z,\theta)p(Z|\theta)}{p(X|Z,\theta^{(i)})}\big) - log\ p(X|\theta^{(i)}) \\ \geqslant \sum_z p(Z|X,\theta^{(i)})log\frac{p(X|Z,\theta)p(Z|\theta)}{p(Z|X,\theta^{(i)})} - log\ p(X|\theta^{(i)}) \\ = \sum_z p(Z|X,\theta^{(i)})log\frac{p(X|Z,\theta)p(Z|\theta)}{p(Z|X,\theta^{(i)})p(X|\theta^{(i)})}$$

令：

$$B(\theta,\theta^{(i)}) = \mathcal{L}(\theta^{(i)})+\sum_z p(Z|X,\theta^{(i)})log\frac{p(X|Z,\theta)p(Z|\theta)}{p(Z|X,\theta^{(i)})p(X|\theta^{(i)})}$$
则：
$$\mathcal{L}(\theta) \geqslant B(\theta,\theta^{(i)})$$
即函数 $B(\theta,\theta^{(i)})$是 $\mathcal{L}(\theta)$的一个下届，因此任何使 $B(\theta,\theta^{(i)})$增大的 $\theta$都会使 $\mathcal{L}(\theta)$的下届增大，也即使 $\mathcal{L}(\theta)$增大，为了使 $\mathcal{L}(\theta)$尽可能的大，自然选择使 $B(\theta,\theta^{(i)})$达到最大的 $\theta^{(i+1)}$，即：
$$\theta^{(i+1)} = arg \max_\theta B(\theta,\theta^{(i)})$$

现在求$\theta^{(i+1)}$的表达式，省去对$\theta$的极大化而言的常数的项，有：

$$\theta^{(i+1)} = arg \max_\theta \big(\mathcal{L}(\theta^{(i)})+\sum_z p(Z|X,\theta^{(i)})log\frac{p(X|Z,\theta)p(Z|\theta)}{p(Z|X,\theta^{(i)})p(X|\theta^{(i)})} \big) \\ = arg \max_\theta \big(\sum_z p(Z|X,\theta^{(i)})log \ (p(X|Z,\theta)p(Z|\theta))\big) \\ = arg \max_\theta \big(\sum_z p(Z|X,\theta^{(i)})log \ (p(X,Z|\theta)\big) \\ = arg \max_\theta Q(\theta,\theta^{(i)})$$

$EM$算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

上图给出$EM$算法的直观解释。图中上方曲线为$\mathcal{L}(\theta)$，下方曲线为$B(\theta,\theta^{(i)})$，图中两个函数在点$\theta = \theta^{(i)}$处相等，但此时$\mathcal{L}(\theta)$并没用达到极大值，这时使$B(\theta,\theta^{(i)})$极大化，由于$\mathcal{L}(\theta) \geqslant B(\theta,\theta^{(i)})$，所以$\mathcal{L}(\theta)$在每次迭代中也是增加的。
在这个过程中，$\mathcal{L}(\theta)$不断增大，从图中可以看出$EM$算法不能保证找到全局最优解。

使用迭代必须保证$p(X|\theta^{(i)}）$单调递增：

$$log \ p(X|\theta^{(i+1)}) \geqslant log \ p(X|\theta^{(i)})$$
证明：
由：
$$p(X) = \frac{p(X,Z)}{p(Z|X)}$$
得到：
$$E[log\ p(X|\theta)] = E[log \ p(X,Z|\theta)-log \ p(Z|X,\theta)]$$

求期望，连续变量即求积分，离散变量即求和，上面我们用的都是求和，这个使用求积分，是为了让大家理解上面求和符号怎么来的。

上式左端：

$$= \int_z log \ p(X|\theta)p(Z|X,\theta^{(i)})dz = log \ p(X|\theta)$$
上式右端：
$$= \int_z log \ p(X,Z|\theta)p(Z|X,\theta^{(i)})dz-\int_z log \ p(Z|X,\theta)p(Z|X,\theta^{(i)}dz$$
上式减号左边的式子与上面 $\theta^{(i+1)}$中的式子一样，我们将其设为 $Q(\theta, \theta^{(i)})$，减号后面一项设为 $H(\theta, \theta^{(i)})$

那么对数似然函数（等式）可写为：

$$log \ p(X|\theta) = Q(\theta, \theta^{(i)}) - H(\theta, \theta^{(i)})$$

因为求得是极大似然函数估计，那么$Q(\theta^{(i+1)}, \theta^{(i)}) \geqslant Q(\theta, \theta^{(i)})$；

假如对于所有的$\theta$都有$H(\theta^{(i)}, \theta^{(i)}) \geqslant H(\theta, \theta^{(i)})$
=> 那么$H(\theta^{(i)}, \theta^{(i)}) \geqslant H(\theta^{(i+1)}, \theta^{(i)})$

这样一来，对数似然函数就会逐渐增大。

那么我们又该如何证明对于所有的$\theta$都有$H(\theta^{(g)}, \theta^{(g)}) \geqslant H(\theta, \theta^{(g)})$

证明：

$$H(\theta^{(i)}, \theta^{(i)}) - H(\theta, \theta^{(i)}) \geqslant 0 \\ = \int_z log \ p(Z|X,\theta^{(i)})p(Z|X,\theta^{(i)})dz - \int_z log \ p(Z|X,\theta)p(Z|X,\theta^{(i)})dz \\ = \int_z log(\frac{p(Z|X,\theta^{(i)})}{p(Z|X,\theta)})p(Z|X,\theta^{(i)})dz \\ = - \int_z log(\frac{p(Z|X,\theta)}{p(Z|X,\theta^{(i)})})p(Z|X,\theta^{(i)})dz$$
利用 $Jensen$不等式（函数的期望大于等于期望的函数）：
$$\geqslant -log \int_z\frac{p(Z|X,\theta)}{p(Z|X,\theta^{(i)})}p(Z|X,\theta^{(i)})dz = 0$$

总结

1. $EM$算法对初值比较敏感，初值的选择非常重要，通常选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

2. 一般条件下$EM$算法是收敛的，但不能保证收敛到全局最优

参考资料：
李航《统计学习方法》