全面理解EM算法

EM算法的引入

EM算法用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计。和其他估计方法不同之处就在于，最大化观测数据的似然函数变成最大化完全数据的似然函数。

如果概率模型的变量都是观测变量，那给定数据就能直接用极大似然估计或贝叶斯估计来估计模型参数。但是，当模型含有隐变量时，就不能简单使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法或极大后验概率估计法。

EM算法

《统计学习方法》中用三硬币模型引出了EM算法，该模型描述如下：

三硬币模型用表达式可以写成如下形式
$$\begin{array}{rl}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{array}$$
y是观测变量（0或1），z是隐变量，表示未观测到的掷硬币A的结果，$\theta =(\pi ,p,q)$是模型参数。这个模型就是以上数据的生成模型。最后要求的就是联合的分布，也就是把z也考虑在内。
将观测数据表示为$Y={\left(Y_1,Y_2,\cdots ,Y_n\right)}^{\mathrm{T}}$，未观测数据表示为$Z={\left(Z_1,Z_2,\cdots ,Z_n\right)}^{\mathrm{T}}$。则观测数据似然函数为：
$P(Y|\theta )={\sum }_{Z}P(Z|\theta )P(Y|Z,\theta )$
即
$$P(Y|\theta )=\prod _{j=1}^n\left[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}\right]$$
于是就考虑求模型参数$\theta =(\pi ,p,q)$的极大似然估计，即
$$\hat{\theta }=\arg \underset{\theta }{\max }\log P(Y|\theta )$$
这个问题没有解析解，只有通过迭代的方法求解。由此就引出了EM算法。在给出针对以上问题的EM算法之前，先对整个问题进行一个求解思路的概括。

最终目标是估计出模型参数求得具体概率预测模型。要是没隐变量那直接用极大似然估计问题就解决了，但是有隐变量。所以我们的目标是最大化完全数据（观测变量+隐变量）的似然函数，完全数据的似然函数就是联合分布，也就是求联合分布$P(Y,Z|\theta )$的极大似然估计就可以估计出模型的参数了。观测变量y的值是知道的，但是Z的情况我们是不知道的，所以EM算法思想就是用Z的期望来替代所有的$Z_i$。那么怎么求期望？需要先知道Z的分布才能求期望，但是这个Z的分布又取决于参数$\theta =(\pi ,p,q)$以及随机变量Y的值（Z的分布就是联合分布的边缘分布）。随机变量Y是观测到的属于已知，但是参数$\theta =(\pi ,p,q)$未知，所以又要求参数，那就回到了一开始的目标。这样就陷入了死循环了。

如何打破这个死循环让思路可以继续下去？我们可以先给参数$\theta =(\pi ,p,q)$一个任意的初值，在这个任意初值的情况下先求得一个Z的期望值（E步），再利用Y和Z组成的已知的似然函数来极大似然估计参数$\theta =(\pi ,p,q)$，回过头来更正这个参数（M步）。但是这样的一次更正肯定达不到目标，那就不断迭代反复更正这个参数$\theta$，直到最后$\theta$收敛至某个值（有EM算法的收敛性证明），这个值很有可能就会是真正的模型参数了（也有可能是局部最优解）。

上面就是EM算法的核心思想，接下来就用这个思想来解决一开始提到的三硬币模型吧。

首先选取参数的初值，记作${\theta }^{(0)}=\left({\pi }^{(0)},p^{(0)},q^{(0)}\right)$，然后通过下面的步骤迭代计算参数的估计值，直至收敛。第i次迭代参数的估计值为${\theta }^{(i)}=\left({\pi }^{(i)},p^{(i)},q^{(i)}\right)$。EM算法的第i+1次迭代如下。
E步：计算在模型参数${\theta }^{(i)}=\left({\pi }^{(i)},p^{(i)},q^{(i)}\right)$下第j个观测数据$y_j$来自硬币B的概率
$${\mu }_j^{(i+1)}=\frac{{\pi }^{(i)}{\left(p^{(i)}\right)}^{y_j}{\left(1-p^{(i)}\right)}^{1-y_j}}{{\pi }^{(i)}{\left(p^{(i)}\right)}^{y_j}{\left(1-p^{(i)}\right)}^{1-y_j}+\left(1-{\pi }^{(i)}\right){\left(q^{(i)}\right)}^{y_j}{\left(1-q^{(i)}\right)}^{1-y_j}}$$
分子是A扔出正面后接着扔B的概率分布，分母是A扔出正面后接着扔B的概率分布 + A扔出反面后接着扔C的概率分布。
M步：计算模型参数的新估计值
$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}$$
掷硬币A掷出正面的概率，也就是所有样本结果里来自B的结果占的比率，也就是单个样本来自B的概率的总和的平均值。
$$p^{(i+1)}=\frac{{\sum }_{j=1}^n{\mu }_j^{(i+1)}{y}_j}{{\sum }_{j=1}^n{\mu }_j^{(i+1)}}$$
B扔正面的概率就是属于B的结果里，观测变量Y=1的个数
$$q^{(i+1)}=\frac{{\sum }_{j=1}^n\left(1-{\mu }_j^{(i+1)}\right)y_j}{{\sum }_{j=1}^n\left(1-{\mu }_j^{(i+1)}\right)}$$
C扔正面的概率就是不属于B的结果（只能属于C）里，观测变量Y=1的个数。

值得注意的是，EM算法与初值的选择有关，不同初值可能得到不同的参数估计值。《统计学习方法》中对EM算法的概括如下：

可以看到，求期望的Q函数是EM算法的核心，书中对其定义以及对所给出的算法概括的说明如下：

EM算法的导出

上面是通过实例对EM算法思路的叙述。那么为什么EM算法能近似实现对观测数据的极大似然估计？下面内容就通过近似求解观测数据的对数似然函数的极大化问题来导出EM算法，看清其作用。
我们面对一个含有隐变量的概率模型，目标是极大化观测数据Y关于参数的对数似然函数，即极大化
$$\begin{array}{rl}L(\theta )& =\log P(Y|\theta )=\log \\ & \end{array}$$