首先明白EM算法是什么，它是一个优化算法，而且大多数情况针对的是最大似然问题。其实从名字也可以看出EM算法的用途，最大期望算法，“期望”二字表明算法主要针对概率问题，“最大”说明这是在求解优化问题。
我们从哪儿说起呢，既然提到概率问题，那我们就从概率问题中参数估计最常见的方法最大似然估计说起吧。最大似然估计基于一个假设，观察到的数据出现概率最大，这个假设其实很符合我们的直观认识，因为你都观察到这组数据了那么它们的出现概率自然也不小，使其最大也挺合理的。其实，有人说过统计学就是猜上帝的游戏，那我们现在开始和上帝来玩一个游戏：上帝手上一枚硬币，但这么硬币不是均匀的。现在上帝抛掷这么硬币，然后告诉你其中有10次正面向上，有90次反面向上，问你上帝再抛掷一次是正面向上还是反面向上。我们学过最大似然估计，所以我们打算用这种方法来进行估计，首先建立概率模型：$f(k;n,p)=Pr(k=k)=\dbinom{n}{k}p^k(1-p)^{n-k}$，然后进行参数估计：$L(\theta)=L(x_1,…,X_n;\theta)=\prod_{i=1}^np(x_i;\theta),{\theta}{\in}{\Theta}$，$\hat{\theta}=argmax_{\theta}L(\theta)$。针对我们和上帝玩的游戏：$\hat{\theta}=argmax_{\theta}\dbinom{100}{10}{\theta}^{10}(1-\theta)^{90}$。
好了，现在上帝发现我们还挺聪明，决定提升游戏难度，Game V2.0：上帝有一台神奇的机器，这台机器类似于我们玩的老虎机，每摇一下手柄，机器就会产生一个1到6的随机数字（概率不相同）。这台机器神奇的地方在于它不仅会产生一个数字还会产生噪声（Noise），而且噪声的分贝是服从高斯分布的，产生的数字不同高斯分布也不同。现在上帝用开始玩游戏，你有仪器记录下每次噪声的分贝，你获得了100次的数据。然后上帝再玩一次，你依旧得到噪声的分贝，这时上帝问你机器产生的数字是多少？我们一想，这其实也不是什么很难的问题啊，继续使用最大似然估计。构建概率模型：$p(x^{(i)};\phi,\mu,\Sigma)=\sum_{z^{(i)=1}}^k{\quad}p(x^{(i)}{\mid}z^{(i)};\mu,\Sigma)p(z^{(i)};\phi)$，由于$p(x^{(i)};\phi,\mu,\Sigma)$的概率有隐含变量$z^{(i)}$和一个依赖$z^{(i)}$的高斯分布决定，所以这个也叫混合高斯模型。概率模型构建好了，我们进行参数估计：$\ell(\phi,\mu,\Sigma)=\sum_{i=1}^mlogp(x^{(i)};\phi,\mu,\Sigma)=\sum_{i=1}^mlog\sum_{z^{(i)=1}}^k{\quad}p(x^{(i)}{\mid}z^{(i)};\mu,\Sigma)p(z^{(i)};\phi)$。我们一看这个对数似然函数没有办法直接用求导然后令倒数等于零来求解参数啊。那怎么办，没事儿，我们还有一招——EM算法。
下面我们就来看一看如何使用EM算法对混合高斯模型进行参数估计。为了便于书写，我们把所有参数都表示成$\theta$。那么我们的优化目标是：$$\hat{\theta}=argmax\ell(\theta)=argmax\sum_{i=1}^mlogp(x^{(i)};\theta)\\
=argmax\sum_{i=1}^mlog\sum_{z}p(x^{(i)},z^{(i)};\theta)$$
对于这样一个函数，我们可不可以考虑找到它的一个下界函数，然后求解它的下界函数的最大值来不断逼近原函数的最大值，这其实就是EM算法的思想，要找函数的下界需要用到另外一个工具—— Jensen不等式(对于凹函数，有$E[f(x)]{\le}f(E[x])$)。利用Jensen不等式对原函数进行放缩：$$\sum_{i=1}^mlogp(x^{(i)};\theta)=\sum_{i=1}^mlog\sum_{z^{(i)}}{\quad}p(x^{(i)},z^{(i)};\theta)\\
=\sum_{i=1}^mlog\sum_{z^{(i)}}{\quad}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\\
\ge\sum_{i=1}^m\sum_{z^{(i)}}{\quad}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$$
这里要求$Q_i(z^{(i)})$是一个概率分布，这样$\sum_{z^{(i)}}{\quad}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$就是函数$\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$的期望。那么现在问题来了，$Q_i(z^{(i)})$的选择有很多，到底应该选择怎么样概率分布啦。假设𝜃已经给定，那么$\ell(\theta)$的值就决定于$𝑄_𝑖(𝑧^{(𝑖)})$和$p(𝑥^{(𝑖)},𝑧^{(𝑖)})$了。我们可以通过调整这两个概率使下界不断上升，以逼近$\ell(\theta)$的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率能够等价于$\ell(\theta)$了，如下图所示
按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：$\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c$,
推导一下，可得，
$$Q_i(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta)}{\sum_zp(x^{(i)},z^{(i)};\theta)}\\
=\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)}\\
=p(z^{(i)}{\mid}x^{(i)};\theta)$$
综上所述，E步骤：$Q_i(z^{(i)})=p(z^{(i)}{\mid}x^{(i)};{\theta}^t)$,
M步骤：${\theta}^{t+1}=argmax_{\theta}\sum_i\sum_zQ_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$
回到我们之前的和上帝的游戏中，具体优化算法为：
为了更加直观的演示EM算法的优化过程，我们可以写一个程序来模型混合高斯模型产生数据的过程，再用EM算法对产生的数据做最大似然估计，github传送门：EM Algorithm
下面看一下随着迭代次数的增加，对数似然值得变化情况吧，
我们可以放大1800到2300次迭代之间的对数似然值变化情况，发现2000次知乎对数似然值几乎就不怎么变化了，说明已达到最大值，