统计学习笔记九----EM算法

前言

　　EM算法是一种迭代算法，1977年由Dempster等人总结提出，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望（expection）；M步，求极大值（maximization），所以这一算法称为期望极大算法（exception maximization algorithm），简称EM算法。

极大似然估计

　　极大似然估计，只是一种概率论在统计学中的应用，它是参数估计的方法之一。说的是已知某个随机样本满足某种概率分布，但是其中具体的参数不清楚，参数估计就是通过若干次实验，观察其结果，利用结果推出参数的大概值。最大似然估计是建立在这样的思想上：已知某个参数能使这个样本出现的概率值最大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值。

最大似然估计你可以把它看作是一个反推。多数情况下我们是根据已知条件来推算结果，而最大似然估计是已经知道了结果，然后寻求使该结果出现的可能性最大的条件，以此作为估计值。比如，如果其他条件一定的话，抽烟者发生肺癌的危险时不抽烟者的5倍，那么如果现在我已经知道有个人是肺癌，我想问你这个人抽烟还是不抽烟。你怎么判断？你可能对这个人一无所知，你所知道的只有一件事，那就是抽烟更容易发生肺癌，那么你会猜测这个人不抽烟吗？我相信你更有可能会说，这个人抽烟。为什么？这就是“最大可能”，我只能说他“最有可能”是抽烟的，“他是抽烟的”这一估计值才是“最有可能”得到“肺癌”这样的结果。这就是最大似然估计。

求最大似然函数估计值的一般步骤：

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

最大（极大）似然估计也是统计学习中经验风险最小化（RRM）的例子。如果模型为条件概率分布，损失函数定义为对数损失函数，经验风险最小化就等价于最大似然估计。

综之，最大似然估计就是给定模型（含有未知参数）和样本集的情况下，用来估计模型参数的方法。其基本思想是找到最佳的模型参数，使得模型实现对样本的最大程度的拟合，也就是使样本集出现的可能性最大。这一方法是基于这样的思想：我们所估计的模型参数，要使得产生这个给定样本的可能性最大。在最大释然估计中，我们试图在给定模型的情况下，找到最佳的参数，使得这组样本出现的可能性最大。举个极端的反面例子，如果我们得到一个中国人口的样本，男女比例为3:2，现在让你估计全国人口的真实比例，你肯定不会估计为男：女=1:0。因为如果是1:0，不可能得到3:2的样本。我们大多很容易也估计为3:2，为什么？样本估计总体，其实它背后的思想就是最大似然。

如果你还是不是很理解最大似然，可以参考以下实例链接：
参考一
参考二
参考三

EM算法的引入

我们先来看看一下书上的引入实例：

我们把中间的推导公式写一下：

由以上三个等式可以得到：

注意：这里的随机变量y是观测变量，表示一次实验观测的结果是1或者0；随机变量z是隐变量，表示的实未观测到的掷硬币A的结果；是模型参数。这一模型是以上数据的生成模型。注意，随机变量y的数据可以观测，随机变量z的数据不可观测到。

如果我们将这n次独立实验，将观测数据表示为向量，将未观测到的数据表示为向量，则观测数据的似然函数为:

由于观测到的n个向量是独立同分布的，所以概率是连乘的形式：

考虑求模型参数的极大似然估计，也就是说求得一个最优参数使得以上条件概率最大（或者对数条件概率最大）。即：

这个问题没有解析解，只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代算法。求解过程如下所示，其实，理解一下求解过程，你会发现EM算法的思想与坐标上升算法或者SMO算法思想是一致的，其实他们思想的本质是相同的。

一般的，用Y表示观测随机变量的数据，Z表示隐藏随机变量的数据。Y和Z连在一起称为完全数据（complete-data），观测数据Y又称为不完全数据（imcomplete-data）。假设给定观测数据Y，其概率分布是，其中是需要估计的模型参数，那么不完全数据Y的似然函数是,对数似然函数是；假设Y和Z的联合概率分布是,那么完全数据的对数似然函数是.

EM算法通过迭代求的极大似然估计，每次迭代包含两步：E步，求期望，M步，求极大化。下面来介绍EM算法。

我们还是直接上图吧，（那些公式我是真的不想编辑了….)

EM算法的推导

假设我们有一个样本集{x(1),…,x(m)}，包含m个独立的样本。但每个样本i对应的类别z(i)是未知的（相当于聚类），也即隐含变量。故我们需要估计概率模型p(x,z)的参数θ，但是由于里面包含隐含变量z，所以很难用最大似然求解，但如果z知道了，那我们就很容易求解了。

对于参数估计，我们本质上还是想获得一个使得似然函数最大化的那个参数θ，现在与最大化似然不同的只是似然函数式中多了一个未知的变量z。也就是说我们的目标是：找到合适的θ和z让L(θ)最大。那我们也许会想，你就是多了一个未知的变量而已啊，我们也可以分别对未知的变量θ和z分别求偏导，然后再令其为0，求解出来不也一样吗？

本质上我们是需要最大化（1）式（对（1）式，我们回忆下联合概率密度下某个变量的边缘概率密度函数的求解，注意这里z也是随机变量。对每一个样本i的所有可能类别z求等式右边的联合概率密度函数和，也就得到等式左边为随机变量x的边缘概率密度），也就是似然函数，但是可以看到里面有“和的对数”，求导后形式会非常复杂（自己可以想象下log(f1(x)+ f2(x)+ f3(x)+…)复合函数的求导），所以很难求解得到未知参数z和θ。那OK，我们可否对（1）式做一些改变呢？我们看（2）式，（2）式只是分子分母同乘以一个相等的函数，还是有“和的对数”啊，还是求解不了，那为什么要这么做呢？咱们先不管，看（3）式，发现（3）式变成了“对数的和”，那这样求导就容易了。我们注意点，还发现等号变成了不等号，为什么能这么变呢？这就是Jensen不等式的大显神威的地方。

　jensen不等式

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

Jensen不等式表述如下：
如果f是凸函数，X是随机变量，那么：E[f(X)]>=f(E[X])
特别地，如果f是严格凸函数，当且仅当X是常量时，上式取等号。

如果用图表示会很清晰：

图中，实线f是凸函数，X是随机变量，有0.5的概率是a，有0.5的概率是b。（就像掷硬币一样）。X的期望值就是a和b的中值了，图中可以看到E[f(X)]>=f(E[X])成立。

当f是（严格）凹函数当且仅当-f是（严格）凸函数。

Jensen不等式应用于凹函数时，不等号方向反向。
回到公式（2），因为f(x)=log x为凹函数（其二次导数为-1/x^2<0）。

（2）式中是的期望，（考虑到，是X的函数，则，又，所以就可以得到公式（3）的不等式了：

OK，到这里，现在式（3）就容易地求导了，但是式（2）和式（3）是不等号啊，式（2）的最大值不是式（3）的最大值啊，而我们想得到式（2）的最大值，那怎么办呢？

现在我们就需要一点想象力了，上面的式（2）和式（3）不等式可以写成：似然函数L(θ)>=J(z,Q)，那么我们可以通过不断的最大化这个下界J，来使得L(θ)不断提高，最终达到它的最大值。

见上图，我们固定θ，调整Q(z)使下界J(z,Q)上升至与L(θ)在此点θ处相等（绿色曲线到蓝色曲线），然后固定Q(z)，调整θ使下界J(z,Q)达到最大值（到），然后再固定θ，调整Q(z)……直到收敛到似然函数L(θ)的最大值处的θ*。这里有两个问题：什么时候下界J(z,Q)与L(θ)在此点θ处相等？为什么一定会收敛？

首先第一个问题，在Jensen不等式中说到，当自变量X是常数的时候，等式成立。而在这里，即：

再推导下，由于（因为Q是随机变量z(i)的概率密度函数），则可以得到：分子的和等于c（分子分母都对所有z(i)求和：多个等式分子分母相加比值不变，这个认为每个样例的两个概率比值都是c），则：

至此，我们推出了在固定参数θ后，使下界拉升的Q(z)的计算公式就是后验概率，解决了Q(z)如何选择的问题。这一步就是E步，建立L(θ)的下界。接下来的M步，就是在给定Q(z)后，调整θ，去极大化L(θ)的下界J（在固定Q(z)后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：

EM的算法流程：

初始化分布参数θ；重复以下步骤直到收敛：E步骤：根据参数初始值或上一次迭代的模型参数来计算出隐性变量的后验概率，其实就是隐性变量的期望。作为隐藏变量的现估计值：

M步骤：将似然函数最大化以获得新的参数值：

这个不断的迭代，就可以得到使似然函数L(θ)最大化的参数θ了。那就得回答刚才的第二个问题了，它会收敛吗？

感性的说，因为下界不断提高，所以极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。理性分析的话，就会得到下面的东西：

具体如何证明的，看推导过程参考

接下来我贴几张图片(最后的参考文献部分也给出了链接)，是我关注的一个人写的，有助于大家的理解：

利用坐标上升法理解EM算法

图中的直线式迭代优化的路径，可以看到每一步都会向最优值前进一步，而且前进路线是平行于坐标轴的，因为每一步只优化一个变量。

这犹如在x-y坐标系中找一个曲线的极值，然而曲线函数不能直接求导，因此什么梯度下降方法就不适用了。但固定一个变量后，另外一个可以通过求导得到，因此可以使用坐标上升法，一次固定一个变量，对另外的求极值，最后逐步逼近极值。对应到EM上，E步：固定θ，优化Q；M步：固定Q，优化θ；交替将极值推向最大。

助于理解EM算法的例子：

实例一

这是一个抛硬币的例子，H表示正面向上，T表示反面向上，参数θ表示正面朝上的概率。硬币有两个，A和B，硬币是有偏的。本次实验总共做了5组，每组随机选一个硬币，连续抛10次。如果知道每次抛的是哪个硬币，那么计算参数θ就非常简单了，如下图所示。

如果不知道每次抛的是哪个硬币呢？那么，我们就需要用EM算法，基本步骤为：1、给θA和θB一个初始值；2、（E-step）估计每组实验是硬币A的概率（本组实验是硬币B的概率=1-本组实验是硬币A的概率）。分别计算每组实验中，选择A硬币且正面朝上次数的期望值，选择B硬币且正面朝上次数的期望值；3、（M-step）利用第三步求得的期望值重新计算θA和θB；4、当迭代到一定次数，或者算法收敛到一定精度，结束算法，否则，回到第2步。

稍微解释一下上图的计算过程。初始值θA=0.6,θB=0.5。

图中的0.45是怎么得来的呢？由两个硬币的初始值0.6和0.5，容易得出投掷出5正5反的概率是，, 0.45就是0.449近似而来的，表示第一组实验选择的硬币是A的概率为0.45。图中的2.2H，2.2T是怎么得来的呢？ ，，表示第一组实验选择A硬币且正面朝上次数的期望值是2.2。其他的值依次类推。

Python 实现

现在我们用Python实现上述例子中的第一次迭代过程，也就是对应上图的（1）、（2）、（3）步，如下所示：

# -*- encoding:utf-8 -*-import numpyfrom scipy.stats import binom#模拟实现实例中的一次迭代过程：掷硬币def em_single(priors,observations):    '''    EM算法的一次迭代    :param priors: 初始化参数theta_A,theta_B  [theta_A,theta_B]    :param observations: 观测矩阵：M*N:M代表实验组数，N:代表每组实验的次数    :return:[new_theta_A,new_theta_B]    '''    counts={'A':{'H':0,'T':0},'B':{'H':0,'T':0}}    theta_A=priors[0]    theta_B=priors[1]    # E step:    for observation in observations:        len_observation = len(observation)        num_heads = observation.sum()#正面        num_tails = len_observation-num_heads#反面        #二项分布求解公式,抛硬币是一个二项分布        contribution_A = binom.pmf(num_heads,len_observation,theta_A)        contribution_B = binom.pmf(num_heads,len_observation,theta_B)        #将两个概率正规化，得到数据来自硬币A，B的概率。        weight_A = contribution_A / (contribution_A + contribution_B)        weight_B = contribution_B / (contribution_A + contribution_B)        #更新在当前参数下A，B硬币产生的正反面次数        counts['A']['H'] += weight_A * num_heads        counts['A']['T'] += weight_A * num_tails        counts['B']['H'] += weight_B * num_heads        counts['B']['T'] += weight_B * num_tails    # M step：    new_theta_A = counts['A']['H'] / (counts['A']['H'] + counts['A']['T'])    new_theta_B = counts['B']['H'] / (counts['B']['H'] + counts['B']['T'])    return [new_theta_A,new_theta_B]if __name__ == '__main__':    #采集数据集，硬币投掷结果，1：表示H 正面；0：表示T 反面    observations = numpy.array([[1,0,0,0,1,1,0,1,0,1],                        [1,1,1,1,0,1,1,1,1,1],                        [1,0,1,1,1,1,1,0,1,1],                        [1,0,1,0,0,0,1,1,0,0],                        [0,1,1,1,0,1,1,1,0,1]])    priors=[0.6,0.5]    new_thetas=em_single(priors,observations)    print new_thetas

实验结果：[0.71301223540051617, 0.58133930831366265]和途中的结果是一致的。

下面我们再给出EM算法的主循环：也就是在收敛前迭代的调用单次EM方法。

'''#EM算法的主循环给定循环的两个终止条件：模型参数变化小于阈值；                     循环达到最大次数'''def em(observations,prior,threshold = 1e-6,iterations=10000):    """    EM算法    ：param observations :观测数据    ：param prior：模型初值    ：param tol：迭代结束阈值    ：param iterations：最大迭代次数    ：return：局部最优的模型参数    """    iteration = 0    while iteration < iterations:        new_prior = em_single(prior,observations)        delta_change = numpy.abs(prior[0]-new_prior[0])        if delta_change < threshold:            break        else:            prior = new_prior            iteration +=1    #返回最终的权重，以及迭代次数    return (new_prior,iteration)

实验结果：

最终得到硬币A正面朝上的概率参数为： 0.796788759383最终得到硬币B正面朝上的概率参数为： 0.519583935675最终的迭代次数为： 14

参考链接：
http://blog.csdn.net/u011300443/article/details/46763743

http://blog.csdn.net/abcjennifer/article/details/8170378

http://blog.csdn.net/zouxy09/article/details/8537620/

https://www.zhihu.com/question/27976634/answer/39132183