EM算法原理详解

参考http://blog.****.net/zouxy09/article/details/8537620 
参考 http://www.cnblogs.com/jerrylead

转载于 http://blog.****.net/yzheately/article/details/51164441

之前介绍了EM算法在混合高斯模型中的应用，现在让我们来看看问什么EM算法可以用于这类问题。 
首先介绍一下Jensen 不等式

Jensen 不等式

我们知道，如果设  是定义域为实数的函数，如果对于所有的实数，，那么  是凸函数。 
显然我们的样本,是有很多属性的，也就是说函数的输入是一个向量。这时是凸函数就等价为为的 hessian 矩阵 H 是半正定的（ H ≥ 0）。

begin-补充-hessian矩阵

对于一个实值多元函数 ，如果函数 的二阶偏导数都存在，则定义  的hessian矩阵为:

其中 D_i表示对第 个变量的微分算子，,上式展开成矩阵形式如下： 

可见如果hessian矩阵存在那么它必然是对称的因为求偏导数时的求导顺序并不影响最终结果： 

利用hessian进行多元函数极值的判定:

如果实值多元函数二阶连续可导，我们可以利用某个临界点M处的hessian矩阵判断该临界点是否为极值： 
如果H(M)是正定矩阵，则临界点M处是一个局部的极小值。 
如果H(M)是负定矩阵，则临界点M处是一个局部的极大值。 
如果H(M)是不定矩阵，则临界点M处不是极值。

end-补充-hessian矩阵

如果,那么它的hessian矩阵必然是半正定的，因为hessian矩阵中的每个元素都是有的二阶导数。 
下面给出jensen不等式定理： 
如果 f 是凸函数， X 是随机变量，那么 

特别地，如果  是严格凸函数，那么(以后都将f(E[X])表示成f(EX))当且仅当，即当且仅当 是常量时。

为了便于理解咱们先看下面： 
凸函数的概念： 
【定义】如果函数满足对定义域上任意两个数都有 ，那么为凸函数。 
注意哦开口向下的是凸,开口向上的是凹。 
如果不等式中等号只有 时才成立，我们分别称它们为严格的凹凸函数. 
推广下就是： 
对于任意的凹函数以及其定义域上个数那么都有 
对于任意的凸函数以及其定义域上个数,那么都有 
如果上面凹凸是严格的，那么不等式的等号只有才成立.

其实上面的结论就是我们的jensen不等式，相信大家都见过。 
可将jensen用图形表示如下： 
 
其中就是,就是,显然，在凸函数中。

另外，  是（严格）凹函数当且仅当 是（严格）凸函数。 
Jensen 不等式应用于凹函数时，不等号方向反向，也就是.

EM算法

假如我们有训练样本集{},我们之前求模型的参数的方式是利用似然值： 
 
但是在上一篇中我们提到由于的值是不确定的随机变量，因此不能通过求最大似然值的方式获得参数。但是一般确定了  后，求解就容易了。

EM 是一种解决存在隐含变量优化问题的有效方法。其思想是：不断地建立的下界（ 步），然后优化下界（  步）。这就话没看懂吧！看不懂正常，下面来详细介绍EM算法： 
对于每一个样例 ，让 表示该样例隐含变量  的某种分布，  满足的条件(如果  是连续性的，那么是概率密度函数，需要将求和符号换做积分符号)。 
这样我们可以得到： 
 
注： 
1、（1）到（2）比较直接，就是分子分母同乘以一个相等的函数。 
2、（2）到（3）利用了 Jensen不等式；首先log函数是凹函数。其次根据lazy Statistician规则，可知其实就是的数学期望，可以看作jensen不等式中的,同理可以看作jensen不等式中的；此时根据jensen不等式我们可得： 

begin-补充-Lazy Statistician规则

设  是随机变量 的函数， （  是连续函数），那么 
（1）  是离散型随机变量，它的分布律为。若绝对收敛，则有 

（2） 是连续型随机变量，它的概率密度为，若绝对收敛，则有 

end-补充

因此，对于任何一种分布，上面的式子(3)都给的值确定了一个下限。但是对于的选择，有多种可能，那种更好的？ 
（我们知道，在EM算法中的E步中，我们的是已知的，即在当前条件下可获得的）假设已经给定，那么的值就决定于了。我们可以通过调整这两个概率使下界尽可能大（即使得(3)的值尽可能大）， 以逼近的真实值,显然当(3)等于(2)时的下限最大。根据jensen不等式我们知道当且仅当是常量时等号成立。这里X，我们假设 
其中是一个不依赖于z^{(i)}的常数。我们知道，因此可得. 
进而可得： 
 
再利用条件概率公式可得:

上面的推导有点乱，现在把他们压缩下就是： 

现在我们知道该如何选择啦，的计算公式就是后验概率. 
这一步就是 步，建立的下界。接下来的步，就是在给定后，调整，去极大化的下界（在固定后，下界还可以调整的更大）。 一般的算法的步骤如下： 

下面让我们来证明EM算法的收敛性： 
假定θ(t)和θ(t+1)是 EM 第 t 次和 t+1 次迭代后的结果。 如果我们证明了ℓ(θ(t)) ≤ ℓ(θ(t+1))，也就是说极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。 下面来证明，选定后，我们得到 E 步，为了能取等号，选取如下： 
 
进而等号满足： 
 
然后我们通过最大化上面等式的式右面获得了新的参数. 
此时必然有： 
 
上面第一行是由式（3）得到，即基于 
第（ 5）步利用了 M 步的定义。第（ 5）步利用了 M 步的定义， M 步就是将θ(t)调整到θ(t+1)，即的选择是基于： 
因此（ 5）成立，（ 6）是之前的等式结果。 
这样就证明了ℓ(θ)会单调增加。因此EM算法是收敛的。 
如果我们定义： 
 
从前面的推导中我们知道ℓ(θ) ≥ J(Q, θ)， EM 可以看作是 J 的坐标上升法， E 步固定θ，优化Q， M 步固定Q优化θ。 
下面从别人那偷了张图，很好的说明了EM算法的优化过程： 

添加一个别人的EM算法代码

#include "math.h"
	#include "stdio.h"
	enum {
	NUMBER,
	LABLE,
	P,
	ALL
	};
	double x[][ALL] = {{1.5,0,0},{1.8,0,0},{1.9,1,0},{2.0,1,0},{2.1,1,0},{2.2,1,0},{2.2,1,0},{2.4,0,0},{2.5,1,0},{2.7,0,0},{2.9,1,0},{3.0,1,0},{3.3,0,1},{3.9,0,0},{3.9,1,0},{3.9,1,0},{4.0,0,0},{4.1,1,0},{4.1,0,0},{4.2,0,0},{5.0,1,0}};
	double mean[2] = {0,0};
	double var[2]={0,0};
	int count[2] = {0,0};
	int sample_count = sizeof(x) / (sizeof(double)*ALL);
	void intialize()
	{
	for(int i=0;i<sample_count;++i)
	{
	count[ int(x[i][LABLE]) ]++;;
	}
	for(int i=0;i<sample_count;++i)
	{
	x[i][P] = (count[int(x[i][LABLE])]*1.0)/sample_count;
	}
	}
	void update_mean_var()
	{
	mean[0] = 0; mean[1] = 0;
	for(int i=0;i<sample_count;++i)
	{
	mean[int(x[i][LABLE])] += x[i][NUMBER];
	}
	mean[0] /= count[0];
	mean[1] /= count[1];
	var[0] = 0; var[1] = 0;
	for(int i=0;i<sample_count;++i)
	{
	var[int(x[i][LABLE])] += ( x[i][NUMBER] - mean[int(x[i][LABLE])] ) * ( x[i][NUMBER] - mean[int(x[i][LABLE])] );
	}
	var[0] /= count[0];
	var[1] /= count[1];
	}
	double pi = 3.1415926;
	double e = 2.718281828;
	double Gos(double x, double mean, double var)
	{
	return (1/sqrt(var*pi*2)) * pow(e, -1*(x-mean)*(x-mean)/(2*var));
	}
	void update_p()
	{
	for(int i=0;i<sample_count;++i)
	{
	x[i][P] = 0.0;
	double p1 = Gos(x[i][NUMBER],mean[int(x[i][LABLE])],var[int(x[i][LABLE])]);
	double p2 = Gos(x[i][NUMBER],mean[(int(x[i][LABLE])+1)%2],var[(int(x[i][LABLE])+1)%2]);
	x[i][P] = p1 / (p1 + p2) ;
	}
	}
	int main(void)
	{
	intialize();
	int iter_max = 10;
	printf("before:");
	for(int k=0;k<sample_count;++k)
	{
	printf("\t#%d#,%f",int(x[k][LABLE]),x[k][NUMBER]);
	}
	printf("\n");
	for(int i = 0,j=1;i<iter_max&&j>0;++i)
	{
	j=0;
	update_mean_var();;
	update_p();
	printf("iter:%d",i);
	for(int k=0;k<sample_count;++k)
	{
	x[k][LABLE] = x[k][P]<0.5 ? (int(x[k][LABLE])+1)%2:x[k][LABLE] ;
	j = x[k][P]<0.5 ? j+1:j;
	printf("\t#%d#,%f",int(x[k][LABLE]),x[k][NUMBER]);
	}
	printf("\n");
	}
	return 0;
	}