EM算法1-原理
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|\theta)=\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的模型产生，模型参数为\(\theta\)，我们要最大化下面这个似然:
\(l(\theta)=\displaystyle\sum_{i}^{m}logp(x_i;\theta)=\displaystyle\sum_{i}^{m}log\sum_{z_i}p(x_i,z_i;\theta)\)。
直接求解这个优化问题非常困难。EM算法是通过迭代的方式求解，分为Expectation步和Maximization步。它的主要思想先找到目标函数的下边界，然后逐步提高这个下边界，进而得到一个最优解，但是这个最优解不一定是全局最优的。
下面看一下这个下界是怎么推导出来的---
\(\displaystyle\sum_{i}^{m}logp(x_i;\theta)\)
\(=\displaystyle\sum_{i}^{m}log\sum_{z_i}p(x_i,z_i;\theta)\)
对于i,假设\(Q_i\)是Z上的某个概率分布
\(=\displaystyle\sum_{i}^{m}log\sum_{z_i}Q_i(z_i)\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}\)
\(>=\displaystyle\sum_{i}^{m}\sum_{z_i}Q_i(z_i)log\frac{p(x_i,z_i;\theta)}{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;\theta)}{Q_i(z_i)}\)看做是期望，然后可以套用Jensen不等式[2]，即可得到上述结果。
现在有了下限，但是里面的Qi还不知道。怎么确定Qi呢？如果我们已经有了\(\theta\)的一个猜测值，那么这里自然让下限在\(\theta\)处值和似然函数在\(\theta\)处的值越接近越好，让不等式eq1在\(\theta\)处取得等号。因为log函数是严格凹函数，所以只有在E(X)==X(恒等于)的时候等号才会成立，比如当X是个常数的时候。基于上面的性质，令
\(\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}=c\)
基于此，可以推出
\(\frac{\sum_zp(x_i,z;\theta)}{\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;\theta)}{\sum_zp(x_i,z;\theta)}\)
\(=\frac{p(x_i,z_i;\theta)}{p(x_i;\theta)}\)
\(=p(z_i|x_i;\theta)\)
因此，Qi是给定xi和\(\theta\)下zi的后验概率。
这就是E步骤，总结一下，假设已知\(\theta\)，先求出似然函数的下限，然后求出隐含变量的分布Qi。
在接下来的M步骤，因为E步骤已经得到了一个Qi，这个步骤求最大化eq1的\(\theta\)值，也就是求下限的最大值点。
然后把M步骤求得的\(\theta\)输入E步骤，循环往复，直至收敛。
Repeat until convergence{
E-step:for each i,set
\(Q_i(z_i):=p(z_i|x_i;\theta)\)
M-step:set
\(\theta:=argmax_{\theta}\displaystyle\sum_{i}^{m}\sum_{z_i}Q_i(z_i)log\frac{p(x_i,z_i|\theta)}{Q_i(z_i)}\)
}
下面图片更直观的描述了EM过程，图片来自[4],E步挪动下界到\(\theta\)值与目标函数相同，M步求下界函数的最大值点做为新的\(\theta\)
如果定义\(J(Q,\theta)=\displaystyle\sum_{i}^{m}\sum_{z_i}Q_i(z_i)log\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}\)
那么，EM算法可以看做函数J的坐标轴下降过程，E步最大化Q，M步最大化\(\theta\)。
EM算法是会收敛的，具体的证明参见参考[3]，但是EM算法有可能陷入局部最优的，它对初始值敏感。
下面尝试用EM算法解决文章开始的三硬币问题。
假设已经经过了j步的迭代，现在已经有了\(\theta^j=(r^j,c^j,q^j)\) (为了避免写起来混淆，把参数里面B的正面概率p变成了c)
E步:
这里要求\(p(z_i|x_i;\theta^j)\)，因为是二分问题，为了描述简便，可以直接求正面的概率，根据贝叶斯概率公式：
(为了书写简单，下面把表示迭代次数的上角标j去掉了，心里记得，r,c,q是已知的)
\(p(z_i=1|x_i;\theta)=\frac{p(x_i|z_i=1;\theta)p(z_i=1;\theta)}{p(x_i|z_i=1;\theta)p(z_i=1;\theta)+p(z_i=0;\theta)p(x_i|z_i=0;\theta)}\)
\(=\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;\theta)\)记做\(\mu^{(j+1)}\)表示是第j+1次迭代得到的值，为了书写清楚一下(cnblog对公式支持的有些差啊)，还是把上角标去掉了
M步骤:
现在\(p(z_i|x_i;\theta)\)已经知道，就开始解决下面这个优化问题了
\(J(\theta)=\sum\mu_ilog\frac{p(x_i,z_i=1;\theta)}{\mu_i}+(1-\mu_i)log\frac{p(x_i,z_i=0;\theta)}{1-\mu_i}\)
\(=\sum\mu_ilog\frac{rc^{x_i}(1-c)^{(1-x_i)}}{\mu_i}+(1-\mu_i)log\frac{(1-r)q^{x_i}(1-q)^{(1-x_i)}}{1-\mu_i}\)
令\(\frac{\partial J(\theta)}{r}=0\)
很容易可以得到\(r=\frac{1}{m}\sum\mu_i\)
令\(\frac{\partial J(\theta)}{c}=0\)
同样容易得到\(c=\frac{\sum\mu_ix_i}{\sum\mu_i}\)
令\(\frac{\partial J(\theta)}{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
