EM算法整理及其python实现

此文转载以供以后查阅学习方便，参考原文请移步链接 http://www.hankcs.com/ml/em-algorithm-and-its-generalization.html

另外，参考李航老师的《统计学习方法》和周志华老师的《机器学习》

如果概率模型的变量都是观测变量（数据中可见的变量），则可以直接用极大似然估计，或者用贝叶斯估计模型参数。但是，当模型含有隐变量（数据中看不到的变量）时，就不能简单地使用这些估计方法，而应该使用含有隐变量的概率模型参数的极大似然估计法，也即EM算法。

EM算法的引入

引入EM算法有两个常见的例子，一个是三硬币模型，一个是双硬币模型。《统计学习方法》介绍的是三硬币模型，本文引述该模型后，再补充双硬币模型。

三硬币模型

有ABC三枚硬币，单次投掷出现正面的概率分别为π、p、q。利用这三枚硬币进行如下实验：

1、第一次先投掷A，若出现正面则投掷B，否则投掷C

2、记录第二次投掷的硬币出现的结果，正面记作1，反面记作0

独立重复1和2十次，产生如下观测结果：

1 1 0 1 0 0 1 0 1 1

假设只能观测到掷硬币的最终结果，无法观测第一次投掷的是哪一枚硬币，求π、p、q，即三硬币模型的参数。

记模型参数为θ=(π,p,q)，无法观测的第一次投掷的硬币为随机变量z，可以观测的第二次投掷的硬币为随机变量y，则观测数据的似然函数为：

这是个一目了然的式子，两个事件，第一个事件选出那枚看不到的硬币，第二个事件利用这枚硬币进行一次投掷。利用硬币结果只可能是0或1这个特性，可以将这个式子展开为：

y的观测序列给定了，怎么找出一个模型参数，使得这个序列的概率（似然函数的值）最大呢，也就是求模型参数的极大似然估计：

这个问题我认为是个NP问题，一方面，给定模型参数，可以在多项式时间求出似然函数的值，然而模型参数的组合是无穷的，谁也不知道它是否是最优的。

EM算法简单理解

EM算法是求解这个问题的一种迭代算法（我认为并非精确算法，而是近似算法），它有3步：

初始化：选取模型参数的初值：，循环如下两步迭代

E步：计算在当前迭代的模型参数下，观测数据y来自硬币B的概率：

这个式子也是一目了然的，分子代表选定B并进行一次投掷试验，分母代表选定B或C并进行一次投掷试验，两个一除就得到试验结果来自B的概率。

M步：估算下一个迭代的新的模型估算值：

这个也好说，把这n个{试验结果来自B的概率}求和得到期望，平均后，得到B出正面的似然估计，同理有p和q。

重复迭代，直到收敛为止。

这个模型中，观测数据Y和隐数据Z组合在一起称为完全数据，单独的观测数据Y称为不完全数据。在隐数据未知的情况，无法直接估计Y的概率分布。但当模型概率给定时，就可以估计Y的条件概率分布了。

Y的条件概率分布估计出来后有什么用呢？利用Y的条件概率分布，又可以更新模型参数……那问题来了，为什么要这么做，这么做能否找到最优解，原理是什么？

带着这些问题啃书本稍微有趣一些，在探索这个问题之前，有必要规范地描述EM算法，并引入一些正规的符号和定义：

EM算法的标准定义

输入：观测变量数据Y，隐变量数据Z,联合分布，条件分布；

输出：模型参数θ。

(1)  选择参数的初值,开始迭代；

(2)  E步：记为第i次迭代参数θ的估计值，在第i+1次迭代的E步，计算

这里，是在给