EM 算法

因为时间原因实在是来不及一一看那些书籍，经学长指点，直接学习统计翻译的em算法，然后是基于短语的概率翻译表提取，再者就是去学习bp神经网络。 em在ibm model1中的应用，看了很久才看明白，当我以为自己看懂了的时候，学长让我用一句话概括一下em，我竟一时语塞，才惊觉其实并不是很懂。于是老老实实去翻看李航老师的《统计学习方法》，第九章详细介绍了em算法，也加深了自己的印象。

什么是EM

简单来说，em算法是用来解决含有隐变量的概率模型的参数最大似然估计的迭代算法。 用数学语言来说，就是观测数据已知，求参数$\theta$,那么我们建立数学模型$P(y|\theta)$,其似然函数是L(y,$\theta$)=log$P(y|\theta)$;求L的最大值，当无法得到L的解析解时，要用em来迭代求解。

关键词： 隐变量 概率模型 参数估计 迭代算法

Why EM

其实根据上面对em的理解，大致就可以看出来，em用来进行最大似然参数估计，我们之前会使用最小二乘法进行似然估计，但是当无法得到参数解析解时，我们只能一步步减小参数的值，使得概率模型最优。

三枚硬币的例子

李航老师的书中提到的一个例子，就是抛三枚硬币的例子。有三枚硬币，先抛出硬币A，如果A是正面，就选择B，然后抛B，记录结果，反之则选择C，记录结果；重复上述过程10次，得到如下结果：1101100110. 我们需要做的，就是求出三枚硬币抛出正面的概率p,q, π。

解：对于这样一个问题，我们首先来看看隐变量是什么。对于这个问题，直接观测的变量是y，1101100110. 隐变量则是观测变量对应的A的结果。假设参数用$\theta$来表示，那么整个模型就可以表示为：

$P(y|\theta) =\sum_{z}P(y,z|\theta)$

$=\sum_zP(z|\theta)*P(y|z,\theta)$

$=P(y,z=0|\theta)+P(y,z=1|\theta)$

$=π*p^y*(1-p)^{1-y}+(1-π)*q^y*(1-q)^{1-y}$

这里面的y表示的是每一个样本值，假设用Y表示每一次实验的观测变量值，用Z表示隐变量的值，那么 Y={y1,y2,y3.....y10}，Z={在z1,z2,z3.....z10}. 模型的似然估计函数为：

$P(Y|\theta)=\sum_ZP(Y,Z|\theta)$

也就是：

$P(Y|\theta) =\prod_{yi=1}^{10}P(yi|\theta)$

$=\prod_{yi=1}π*p^{yi}*(1-p)^{1-yi}+(1-π)*q^{yi}*(1-q)^{1-yi}$

$L(y,\theta)=logP(Y|\theta)$

利用最大似然估计法，求解最大的 $\theta$, 即 $\hat{\theta}=arg\,max_{\theta}L(Y,\theta)$ 可以看到这个表达式是没有解析解的，只能通过迭代的方式来更新参数，求得最大解。而em算法则正是用来对此进行求解的一种算法。 首先，需要对参数进行初始化，这里假设初始化为0.5,比较符合常识。 $\theta^0=(π^0=0.5, p^0=0.5,q^0=0.5);$ 然后，采用em算法的e步，求出隐变量，也就是A抛出的是正面的概率。

$P(z^{i+1}|\theta^i,y) = \frac{π^i*(p^i)^{y}*(1-p^i)^{1-y}}{π^i*(p^i)^{y}*(1-p^i)^{1-y}+(1-π^i)*(q^i)^{y}*(1-q^i)^{1-y}}$

其次，利用M步来迭代更新参数。

$π^{i+1} = \frac{1}{n}\sum_{yi=0}^{10}P(z^{i+1}|\theta^i,yi)$

$p^{i+1}= \frac{\sum_{yi=0}^{10}P(z^{i+1}|\theta^i,yi)*yi}{\sum_{yi=0}^{10}P(z^{i+1}|\theta^i,yi)}$

$q^{i+1} =\frac{\sum_{yi=0}^{10}(1-P(z^{i+1}|\theta^i,yi))*yi}{\sum_{yi=0}^{10}(1-P(z^{i+1}|\theta^i,yi))}$

最后重复EM步骤，直到收敛为止。

EM算法的推导

看完上面的例子，似乎还是无法明白这个EM算法到底是怎么回事。 来重新理一下思路，根据上面对问题的描述，其实我们的目的就是要通过迭代的思想，来使得$L(Y,\theta)$达到一个最大值。

$L(Y|\theta)-L(Y|\theta^{i+1})=logP(Y|\theta)-logP(Y|\theta^{i+1})$

$=log\sum_zP(Y,z|\theta)-logP(Y|\theta^{i+1})$

$=log\sum_zP(z|y,\theta^{i+1})*\frac{P(Y|z,\theta)*P(z|\theta)}{P(z|y,\theta^{i+1})}-logP(Y|\theta^{i+1})$

$>=\sum_zP(z|y,\theta^{i+1}) log\frac{P(Y|z,\theta)*P(z|\theta)}{P(z|y,\theta^{i+1})}-logP(Y|\theta^{i+1})$

$=\sum_zP(z|y,\theta^{i+1}) log\frac{P(Y|z,\theta)*P(z|\theta)}{P(z|y,\theta^{i+1})*P(Y|\theta^{i+1})}$

那么， $L(Y|\theta)=\sum_zP(z|y,\theta^{i+1}) log\frac{P(Y|z,\theta)*P(z|\theta)}{P(z|y,\theta^{i+1})*P(Y|\theta^{i+1})}+L(Y|\theta^{i+1})$

可以看到，我们使用Jensen不等式，成功地把和的对数转换为了对数的和，可以求导来得到最大值。下面对上式化简：

$\hat\theta_{max} =arg \,max\,\sum_zP(z|y,\theta^{i+1}) log\frac{P(Y|z,\theta)*P(z|\theta)}{P(z|y,\theta^{i+1})*P(Y|\theta^{i+1})}+L(Y|\theta^{i+1})$

$=arg \,max\,\{\sum_zP(z|y,\theta^{i+1})logP(Y,Z|\theta) \}$

$=arg \,max\, Q(\theta,\theta^{i+1})$

实际上我们要最大化的其实就是这个Q函数，但这个是L的下界，因此em算法就是通过一步步的迭代求解下界，来逼近L的最大值，可以参考梯度下降的思想。下面的一幅图解释了具体的过程。

EM算法步骤

(1) 首先需要对参数$\theta^0$进行初始化。一般会采用均匀分布的方法来初始化参数。当然，参数的初始值直接影响到em参数估计的结果。

(2) E步： 用$\theta^i$表示在i次迭代的$\theta$值。 那么， 这个所谓的期望值，其实指的是，完全数据(观测值和隐变量值)在参数的条件概率下的对数似然函数，对上一次迭代求出的参数和当前观测值y的期望，假设用Q来表示，则在第i+1次迭代中：

$Q(\theta,\theta^i)=E_z(logP(Y,Z|\theta)|Y,\theta^i)$

$=\sum_zlogP(Y,Z|\theta)*P(Z|Y,\theta^i)$

可以看到这个Q函数的来源正是上面推导，所要最大化的表达式。

关于这一部分，其实不是很能理解，按照这个公式来看的确是正确的，难点就在于这个Q函数的确定。看了大多分的例子，其实，一般在E步，会先求解出$P(Z|Y,\theta^i)$，也就是根据当前的参数和观测变量来求出可能的隐变量的分布，然后在M步的时候，再根据隐变量和观测变量来写出参数的估计计算方法。

(3)M步： 这一步主要是求最大值的，更新第i+1次迭代的参数： $\theta^{i+1} =arg\,max_\theta Q(\theta,\theta^i)$ (4)重复(2)(3)步骤，直到$\theta^i,\theta^{i+1}$的差值在一个极小的范围之内。

================ 以上就是整个em的介绍了。拖拖拉拉好几天才把这个算法写完，太拖拉了，好嘛原谅我毕竟是在家写，没动力啊╮(╯▽╰)╭。

不过em的入门还是因为IBM1 model，所以在下一篇会以翻译模型中的em作进一步的案例分析。

本文同步在个人博客，欢迎评论交流 开启掘金成长之旅！这是我参与 「掘金日新计划 · 2 月更文挑战」的第 14 天，点击查看活动详情