em算法python代码_EM算法原理和python简单实现

EM算法推导

假设有样本集

,其中是m个独立的样本,每个样本对应的类别z(i)是不知道的(不知道是好鸡蛋还是坏鸡蛋)，要估计概率模型p(x,)的参数(估计坏鸡蛋的重量),但是由于z我们也不知道，那么就不能使用最大似然估计。那如果知道了z，就好计算了。

目标是求等式(1)左边的最大值,

等式(1): 对每一个样本i的所有可能类别z求等式右边的联合概率密度函数和，也就得到等式左边为随机变量x的边缘概率密度.

如果对等式(1)进行求导,即log(f1(x)+log(f2(x))+.....)，很麻烦，于是有等式(2)的变形

等式(2):分子分母同乘一个不为0的数

等式(3):变成和的对数

插播

Jensen不等式:

如上图是一个凸函数,凸函数定义是什么:

设f是定义域为实数的函数，如果对于所有的实数x。如果对于所有的实数x，f(x)的二次导数大于等于0，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的，那么f是凸函数。如果只大于0，不等于0，那么称f是严格凸函数

如果f是凸函数，X是随机变量，那么：E[f(X)]>=f(E[X])

如果f是凹函数，X是随机变量，那么：E[f(X)]<=f(E[X])

辅助记忆:y=x2是凸函数，y=-x2是凹函数

如等式(3)中

是概率p(x),

是x，

是log(x)，这是个凹函数。

因此等式(2)中的

就是期望(E(X)=∑x*p(x)，f(X)是X的函数，则E(f(X))=∑f(x)*p(x)),而且概率和

，因此对照凹函数的E[f(X)]<=f(E[X])

以及(2)式，可以得到(3)式,模式如下。

但是.......

由于我们要得到(2)式的最大值，(2)和(3)式之间却是个不等式。咋办？我不断的增加(3)式的最大值，是不是相应的让(2)取得的值会更大？水涨船高的原理，你想一想是不是这么个理儿！

是的！

似然函数L(θ)>=J(z,Q)，那么我们可以通过不断的最大化这个下界J，来使得L(θ)不断提高，最终达到它的最大值，如下图所是:

图1

我先让θ固定，下界J(z,Q)它只能跑到上图蓝色的地方，因为他的最大值是红色的L(θ)，然后我再固定Q(z),这个时候调整θ，使得J(z,Q)最大，如下如到了θt+1;然后重复上述动作，直至收敛。

在上述(1)式到(2)式中，我们凭空添加了个Q！

那么这个Q根据什么来选择？

对于上图而言，当θ固定，下界J(z,Q)它只能跑到上图蓝色的地方，(2)式和(3)式是相等的。按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：

分母乘过去,

至此，我们计算出了的概率(在知道θ和样本时，该样本来自于z(i)的概率)，即在固定θ时，为了让(2)式和(3)式相等，即下界拉升(跑到上图1蓝色位置)。这是第一步即E，接下来就是L(θ)的最大化下界J(z,Q)，即调整θ，使得J(z,Q)最大，即步骤M:

θ = maxJ(z,Q)

将J(z,Q)展开就是

上面E步骤和M-步骤，就是EM算法。

对于图1而言，

我们为啥能知道一定能收敛？

当

Jensen不等式相等，这是在θt时刻，那么在θt+1时刻呢？我们要最大化的极大似然是

L(θt+1)有:

上面不等式(4)的右边，就是上式(1)-上式(3)么？

我们在θt时刻，让下界J(z,Q)最大化的时候，让θ变化，从θt跳到θt+1时刻么，这步不就是M的步骤么？于是有:

上面不等式的右边就是L(θt)。

合起来就是:

So:到这儿，EM算法的推导算是结束了。