从似然函数到EM算法(附代码实现)

1. 什么是EM算法

最大指望算法（Expectation-maximization algorithm，又译为指望最大化算法），是在几率模型中寻找参数最大似然估计或者最大后验估计的算法，其中几率模型依赖于没法观测的隐性变量。

最大指望算法通过两个步骤交替进行计算，

第一步是计算指望（E），利用对隐藏变量的现有估计值，计算其最大似然估计值；
第二步是最大化（M），最大化在E步上求得的最大似然值来计算参数的值。M步上找到的参数估计值被用于下一个E步计算中，这个过程不断交替进行。

极大似然估计用一句话归纳就是：知道结果，反推条件θ。

1.1 似然函数

在数理统计学中，似然函数是一种关于统计模型中的参数的函数，表示模型参数中的似然性。“似然性”与“或然性”或“几率”意思相近，都是指某种事件发生的可能性。而极大似然就至关于最大可能的意思。

好比你一位同窗和一位猎人一块儿外出打猎，一只野兔从前方窜过。只听一声枪响，野兔应声到下，若是要你推测，这一发命中的子弹是谁打的？你就会想，只发一枪便打中，因为猎人命中的几率通常大于你那位同窗命中的几率，从而推断出这一枪应该是猎人射中的。

这个例子所做的推断就体现了最大似然法的基本思想。

多数状况下咱们是根据已知条件来推算结果，而最大似然估计是已经知道告终果，而后寻求使该结果出现的可能性最大的条件，以此做为估计值。

1.3 极大似然函数的求解步骤

假定咱们要从10万我的当中抽取100我的来作身高统计，那么抽到这100我的的几率就是(几率连乘)：

\[L(\theta)=L(x_1,...,x_n|\theta)=\prod_{i=1}^{n}p(x_i|\theta),\theta\in\ominus\]

如今要求的就是这个 \(\theta\) 值，也就是使得 \(L(\theta)\) 的几率最大化，那么这时的参数 \(\theta\) 就是所求。

为了便于分析，咱们能够定义对数似然函数，将其变成连加的形式：

\[H(\theta)=lnL(\theta)=ln\prod_{i=1}^{n}p(x_i|\theta)=\sum_{i=1}^{n}lnp(x_i|\theta)\]

对于求一个函数的极值，经过咱们在本科所学的微积分知识，最直接的设想是求导，而后让导数为0，那么解这个方程获得的θ就是了（固然，前提是函数L(θ)连续可微）。但，若是θ是包含多个参数的向量那怎么处理呢？固然是求L(θ)对全部参数的偏导数，也就是梯度了，从而n个未知的参数，就有n个方程，方程组的解就是似然函数的极值点了，最终获得这n个参数的值。

求极大似然函数估计值的通常步骤：

1. 写出似然函数；
2. 对似然函数取对数，并整理；
3. 求导数，令导数为0，获得似然方程；
4. 解似然方程，获得的参数即为所求；

1.4 EM算法

两枚硬币A和B，假定随机抛掷后正面朝上几率分别为PA，PB。为了估计这两个硬币朝上的几率，我们轮流抛硬币A和B，每一轮都连续抛5次，总共5轮：

硬币	结果	统计
A	正正反正反	3正-2反
B	反反正正反	2正-3反
A	正反反反反	1正-4反
B	正反反正正	3正-2反
A	反正正反反	2正-3反

硬币A被抛了15次，在第一轮、第三轮、第五轮分别出现了3次正、1次正、2次正，因此很容易估计出PA，相似的，PB也很容易计算出来(真实值)，以下：

PA = （3+1+2）/ 15 = 0.4
PB= （2+3）/10 = 0.5

问题来了，若是咱们不知道抛的硬币是A仍是B呢（即硬币种类是隐变量），而后再轮流抛五轮，获得以下结果：

硬币	结果	统计
Unknown	正正反正反	3正-2反
Unknown	反反正正反	2正-3反
Unknown	正反反反反	1正-4反
Unknown	正反反正正	3正-2反
Unknown	反正正反反	2正-3反

OK，问题变得有意思了。如今咱们的目标没变，仍是估计PA和PB，须要怎么作呢？

显然，此时咱们多了一个硬币种类的隐变量，设为z，能够把它认为是一个5维的向量（z1,z2,z3,z4,z5)，表明每次投掷时所使用的硬币，好比z1，就表明第一轮投掷时使用的硬币是A仍是B。

• 可是，这个变量z不知道，就没法去估计PA和PB，因此，咱们必须先估计出z，而后才能进一步估计PA和PB。
• 可要估计z，咱们又得知道PA和PB，这样咱们才能用极大似然几率法则去估计z，这不是鸡生蛋和蛋生鸡的问题吗，如何破？

答案就是先随机初始化一个PA和PB，用它来估计z，而后基于z，仍是按照最大似然几率法则去估计新的PA和PB，而后依次循环，若是新估计出来的PA和PB和咱们真实值差异很大，直到PA和PB收敛到真实值为止。

咱们不妨这样，先随便给PA和PB赋一个值，好比：
硬币A正面朝上的几率PA = 0.2
硬币B正面朝上的几率PB = 0.7

而后，咱们看看第一轮抛掷最多是哪一个硬币。
若是是硬币A，得出3正2反的几率为 0.20.20.20.80.8 = 0.00512
若是是硬币B，得出3正2反的几率为0.70.70.70.30.3=0.03087
而后依次求出其余4轮中的相应几率。作成表格以下：

轮数	如果硬币A	如果硬币B
1	0.00512，即0.2 0.2 0.2 0.8 0.8，3正-2反	0.03087，3正-2反
2	0.02048，即0.2 0.2 0.8 0.8 0.8，2正-3反	0.01323，2正-3反
3	0.08192，即0.2 0.8 0.8 0.8 0.8，1正-4反	0.00567，1正-4反
4	0.00512，即0.2 0.2 0.2 0.8 0.8，3正-2反	0.03087，3正-2反
5	0.02048，即0.2 0.2 0.8 0.8 0.8，2正-3反	0.01323，2正-3反

按照最大似然法则：
第1轮中最有可能的是硬币B
第2轮中最有可能的是硬币A
第3轮中最有可能的是硬币A
第4轮中最有可能的是硬币B
第5轮中最有可能的是硬币A

咱们就把几率更大，即更多是A的，即第2轮、第3轮、第5轮出现正的次数二、一、2相加，除以A被抛的总次数15（A抛了三轮，每轮5次），做为z的估计值，B的计算方法相似。而后咱们即可以按照最大似然几率法则来估计新的PA和PB。

PA = （2+1+2）/15 = 0.33
PB =（3+3）/10 = 0.6

就这样，不断迭代 不断接近真实值，这就是EM算法的奇妙之处。

能够期待，咱们继续按照上面的思路，用估计出的PA和PB再来估计z，再用z来估计新的PA和PB，反复迭代下去，就能够最终获得PA = 0.4，PB=0.5，此时不管怎样迭代，PA和PB的值都会保持0.4和0.5不变，因而乎，咱们就找到了PA和PB的最大似然估计。

总结一下计算步骤：

1. 随机初始化分布参数θ

2. E步，求Q函数，对于每个i，计算根据上一次迭代的模型参数来计算出隐性变量的后验几率（其实就是隐性变量的指望），来做为隐藏变量的现估计值：

   \[Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)\]

3. M步，求使Q函数得到极大时的参数取值）将似然函数最大化以得到新的参数值

   \[\theta=argmax\sum_{i}\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\]

4. 而后循环重复二、3步直到收敛。

详细的推导过程请参考文末的参考文献。

2. 采用 EM 算法求解的模型有哪些？

用EM算法求解的模型通常有GMM或者协同过滤，k-means其实也属于EM。EM算法必定会收敛，可是可能收敛到局部最优。因为求和的项数将随着隐变量的数目指数上升，会给梯度计算带来麻烦。

3.代码实现

高斯混合模型 EM 算法

4. 参考文献

如何通俗理解EM算法

做者：@mantchs

GitHub：https://github.com/NLP-LOVE/ML-NLP

欢迎你们加入讨论！共同完善此项目！群号：【541954936】