EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计。
使用EM算法的原因首先举李航老师《统计学习方法》中的例子来说明为什么要用EM算法估计含有隐变量的概率模型参数。
假设有三枚硬币，分别记作A， B， C。这些硬币正面出现的概率分别是$\pi,p,q$。进行如下掷硬币试验：先掷硬币A，根据其结果选出硬币B或C，正面选硬币B，反面边硬币C；然后掷选出的硬币，掷硬币的结果出现正面记作1，反面记作0；独立地重复$n$次试验，观测结果为$\{y_1,y_2,...,y_n\}$。问三硬币出现正面的概率。
三硬币模型（也就是第二枚硬币正反面的概率）可以写作
$ \begin{aligned} &P(y|\pi,p,q) \\ =&\sum\limits_z P(y,z|\pi,p,q)\\ =&\sum\limits_z P(y|z,\pi,p,q)P(z|\pi,p,q)\\ =&\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y} \end{aligned} $
其中$z$表示硬币A的结果，也就是前面说的隐变量。通常我们直接使用极大似然估计，即最大化似然函数
$ \begin{aligned} &\max\limits_{\pi,p,q}\prod\limits_{i=1}^n P(y_i|\pi,p,q) \\ =&\max\limits_{\pi,p,q}\prod\limits_{i=1}^n[\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}]\\ =&\max\limits_{\pi,p,q}\sum\limits_{i=1}^n\log[\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}]\\ =&\max\limits_{\pi,p,q}L(\pi,p,q) \end{aligned} $
分别对$\pi,p,q$求偏导并等于0，求解线性方程组来估计这三个参数。但是，由于它是带有隐变量的，在获取最终的随机变量之前有一个分支选择的过程，导致这个$\log$的内部是加和的形式，计算导数十分困难，而待求解的方程组不是线性方程组。当复杂度一高，解这种方程组几乎成为不可能的事。以下推导EM算法，它以迭代的方式来求解这些参数，应该也算一种贪心吧。
算法导出与理解对于参数为$\theta$且含有隐变量$Z$的概率模型，进行$n$次抽样。假设随机变量$Y$的观察值为$\mathcal{Y} = \{y_1,y_2,...,y_n\}$，隐变量$Z$的$m$个可能的取值为$\mathcal{Z}=\{z_1,z_2,...,z_m\}$。
写出似然函数：
$ \begin{aligned} L(\theta) &= \sum\limits_{Y\in\mathcal{Y}}\log P(Y|\theta)\\ &=\sum\limits_{Y\in\mathcal{Y}}\log \sum\limits_{Z\in \mathcal{Z}} P(Y,Z|\theta)\\ \end{aligned} $
EM算法首先初始化参数$\theta = \theta^0$，然后每一步迭代都会使似然函数增大，即$L(\theta^{k+1})\ge L(\theta^k)$。如何做到不断变大呢？考虑迭代前的似然函数（为了方便不用$\theta^{k+1}$）：
$ \begin{gather} \begin{aligned} L(\theta)=&\sum\limits_{Y\in \mathcal{Y}} \log\sum\limits_{Z\in \mathcal{Z}} P(Y,Z|\theta)\\ =&\sum\limits_{Y\in \mathcal{Y}} \log\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^k)}\\ \end{aligned} \label{} \end{gather} $
至于上式的第二个等式为什么取出$P(Z|Y,\theta^k)$而不是别的，正向的原因我想不出来，马后炮原因在后面记录。
考虑其中的求和
$ \sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)=1$
且由于$\log$函数是凹函数，因此由Jenson不等式得
$ \begin{gather} \begin{aligned} L(\theta) \ge&\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^k)}\\ =&B(\theta,\theta^k) \end{aligned}\label{} \end{gather} $
当$\theta = \theta^k$时，有
$ \begin{gather} \begin{aligned} L(\theta^k) \ge& B(\theta^k,\theta^k)\\ =&\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log\frac{P(Y,Z|\theta^k)}{P(Z|Y,\theta^k)}\\ =&\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log P(Y|\theta^k)\\ =&\sum\limits_{Y\in \mathcal{Y}}\log P(Y|\theta^k)\\ =&L(\theta^k)\\ \end{aligned} \label{} \end{gather} $
也就是在这时，$(2)$式取等，即$L(\theta^k) = B(\theta^k,\theta^k)$。取
$ \begin{gather} \theta^*=\text{arg}\max\limits_{\theta}B(\theta,\theta^k)\label{} \end{gather} $
可得不等式
$L(\theta^*)\ge B(\theta^*,\theta^k)\ge B(\theta^k,\theta^k) = L(\theta^k)$
所以，我们只要优化$(4)$式，让$\theta^{k+1} = \theta^*$，即可保证每次迭代的非递减势头，有$L(\theta^{k+1})\ge L(\theta^k)$。而由于似然函数是概率乘积的对数，一定有$L(\theta) < 0$，所以迭代有上界并且会收敛。以下是《统计学习方法》中EM算法一次迭代的示意图：
进一步简化$(4)$式，去掉优化无关项：
$ \begin{aligned} \theta^*=&\text{arg}\max\limits_{\theta}B(\theta,\theta^k) \\ =&\text{arg}\max\limits_{\theta}\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^k)} \\ =&\text{arg}\max\limits_{\theta}\sum\limits_{Y\in \mathcal{Y}}\sum\limits_{Z\in \mathcal{Z}} P(Z|Y,\theta^k)\log P(Y,Z|\theta) \\ =&\text{arg}\max\limits_{\theta}Q(\theta,\theta^k) \\ \end{aligned} $
$Q$函数使用导数求极值的方程与没有隐变量的方程类似，容易求解。
综上，EM算法的流程为：
1. 设置$\theta^0$的初值。EM算法对初值是敏感的，不同初值迭代出来的结果可能不同。