EM算法用于含有隐含变量的概率模型参数的极大似然估计。什么是隐含变量的概率模型呢？举个例子，假设有3枚硬币，分别记为A,B,C，它们正面出现的概率分别为r,p,q。每次实验先掷硬币A，如果出现的是正面就投B，如果出现的反面就投C，出现正面记为1，出现反面记为0。独立10次实验，观测结果如下：1101001011。如果只有这个结果，而不知道过程，问如何估计r,q,p?也就是说，我们能看到每次的观测结果，可是这个结果是B产生的还是C产生的我们不知道，也就是A的结果我们不知道，这个就是所谓的隐含变量。如果把观测变量用Y来表示，隐含变量（A的结果）用Z表示，那么观测数据的似然函数为：
(P(Y| heta)=prod_i{rp^{y_i}(1-p)^{1-y_i}+(1-r)q^{y_i}(1-q)^{1-y_i}})
把上面的模型泛化一下可以概括为，有观测数据{(x_1,x_2,...x_m)}，由一个具有观测变量X和隐含变量Z的模型产生，模型参数为( heta)，我们要最大化下面这个似然:
(l( heta)=displaystylesum_{i}^{m}logp(x_i; heta)=displaystylesum_{i}^{m}logsum_{z_i}p(x_i,z_i; heta))。
直接求解这个优化问题非常困难。EM算法是通过迭代的方式求解，分为Expectation步和Maximization步。它的主要思想先找到目标函数的下边界，然后逐步提高这个下边界，进而得到一个最优解，但是这个最优解不一定是全局最优的。
下面看一下这个下界是怎么推导出来的---
(displaystylesum_{i}^{m}logp(x_i; heta))
(=displaystylesum_{i}^{m}logsum_{z_i}p(x_i,z_i; heta))
对于i,假设(Q_i)是Z上的某个概率分布
(=displaystylesum_{i}^{m}logsum_{z_i}Q_i(z_i)frac{p(x_i,z_i; heta)}{Q_i(z_i)})
(>=displaystylesum_{i}^{m}sum_{z_i}Q_i(z_i)logfrac{p(x_i,z_i; heta)}{Q_i(z_i)}) ---- (eq1)
这一步用到了Jensen不等式，因为log函数是concave的(二阶导数小于0)，所以有log(E(X))>=E(log(x))。而(Q_i)是概率分布，所以可以把(sum_{z_i}Q_i(z_i)frac{p(x_i,z_i; heta)}{Q_i(z_i)})看做是期望，然后可以套用Jensen不等式[2]，即可得到上述结果。
现在有了下限，但是里面的Qi还不知道。怎么确定Qi呢？如果我们已经有了( heta)的一个猜测值，那么这里自然让下限在( heta)处值和似然函数在( heta)处的值越接近越好，让不等式eq1在( heta)处取得等号。因为log函数是严格凹函数，所以只有在E(X)==X(恒等于)的时候等号才会成立，比如当X是个常数的时候。基于上面的性质，令
(frac{p(x_i,z_i; heta)}{Q_i(z_i)}=c)
基于此，可以推出
(frac{sum_zp(x_i,z; heta)}{sum_zQ_i(z)}=c) (这个很容易推出来,a1/b1=c,a2/b2=c,a3/b3=c => (a2+a2+a3)/(b1+b2+b3)=c)
有，
(Q_i(z_i)=frac{p(x_i,z_i; heta)}{sum_zp(x_i,z; heta)})
(=frac{p(x_i,z_i; heta)}{p(x_i; heta)})
(=p(z_i|x_i; heta))
因此，Qi是给定xi和( heta)下zi的后验概率。
这就是E步骤，总结一下，假设已知( heta)，先求出似然函数的下限，然后求出隐含变量的分布Qi。
在接下来的M步骤，因为E步骤已经得到了一个Qi，这个步骤求最大化eq1的( heta)值，也就是求下限的最大值点。
然后把M步骤求得的( heta)输入E步骤，循环往复，直至收敛。
Repeat until convergence{
E-step:for each i,set
(Q_i(z_i):=p(z_i|x_i; heta))
M-step:set
( heta:=argmax_{ heta}displaystylesum_{i}^{m}sum_{z_i}Q_i(z_i)logfrac{p(x_i,z_i| heta)}{Q_i(z_i)})
}
下面图片更直观的描述了EM过程，图片来自[4],E步挪动下界到( heta)值与目标函数相同，M步求下界函数的最大值点做为新的( heta)
如果定义(J(Q, heta)=displaystylesum_{i}^{m}sum_{z_i}Q_i(z_i)logfrac{p(x_i,z_i; heta)}{Q_i(z_i)})
那么，EM算法可以看做函数J的坐标轴下降过程，E步最大化Q，M步最大化( heta)。
EM算法是会收敛的，具体的证明参见参考[3]，但是EM算法有可能陷入局部最优的，它对初始值敏感。
下面尝试用EM算法解决文章开始的三硬币问题。
假设已经经过了j步的迭代，现在已经有了( heta^j=(r^j,c^j,q^j)) (为了避免写起来混淆，把参数里面B的正面概率p变成了c)
E步:
这里要求(p(z_i|x_i; heta^j))，因为是二分问题，为了描述简便，可以直接求正面的概率，根据贝叶斯概率公式：
(为了书写简单，下面把表示迭代次数的上角标j去掉了，心里记得，r,c,q是已知的)
(p(z_i=1|x_i; heta)=frac{p(x_i|z_i=1; heta)p(z_i=1; heta)}{p(x_i|z_i=1; heta)p(z_i=1; heta)+p(z_i=0; heta)p(x_i|z_i=0; heta)})
(=frac{rc^{x_i}(1-c)^{(1-x_i)}}{rc^{x_i}(1-c)^{1-x_i}+(1-r)q^{x_i}(1-q)^{(1-x_i)}})
把(p(z_i=1|x_i; heta))记做(mu^{(j+1)})表示是第j+1次迭代得到的值，为了书写清楚一下(cnblog对公式支持的有些差啊)，还是把上角标去掉了
M步骤:
现在(p(z_i|x_i; heta))已经知道，就开始解决下面这个优化问题了
(J( heta)=summu_ilogfrac{p(x_i,z_i=1; heta)}{mu_i}+(1-mu_i)logfrac{p(x_i,z_i=0; heta)}{1-mu_i})
(=summu_ilogfrac{rc^{x_i}(1-c)^{(1-x_i)}}{mu_i}+(1-mu_i)logfrac{(1-r)q^{x_i}(1-q)^{(1-x_i)}}{1-mu_i})
令(frac{partial J( heta)}{r}=0)
很容易可以得到(r=frac{1}{m}summu_i)
令(frac{partial J( heta)}{c}=0)
同样容易得到(c=frac{summu_ix_i}{summu_i})
令(frac{partial J( heta)}{q}=0)
同样容易得到(c=frac{sum(1-mu_i)x_i}{sum(1-mu_i)}) 参考[1][5]
参考：
[1]李航《统计学习方法》
[2]Jensen不等式：http://www.cnblogs.com/naniJser/p/5642288.html
[3]Andrew Ng机器学习课程的讲义:http://cs229.stanford.edu/notes/cs229-notes8.pdf
[4]介绍EM算法的blog:http://blog.csdn.net/zouxy09/article/details/8537620
[5]http://chenrudan.github.io/blog/2015/12/02/emexample.html