机器学习之---EM算法

http://blog.csdn.net/pipisorry/article/details/42550815

EM算法有很多的应用，最广泛的就是GMM混合高斯模型、聚类、HMM、基于概率的PLSA模型等等。

本文详细讲述EM算法的由来、EM算法的实现思路、EM算法解决PLSA和LDA的方法。
概述

EM是一种解决存在隐含变量优化问题的有效方法。EM的意思是“Expectation Maximization”，与最大似然估计MLE的关系，EM是解决（不完全数据的）MLE问题的迭代算法 iterative algorithm。

皮皮blog

 

EM算法Intuition示例
概率论中的EM算法：一个理解EM算法的简单示例

Note: 这里的四个概率相加和恒为1。
使用最大似然估计MLE求解
求最大似然函数估计值的一般步骤

（1）写出似然函数；（lz所有事件发生的概率）

（2）对似然函数取对数，并整理；

（3）求导数，令导数为0，得到似然方程；

（4）解似然方程，得到的参数即为所求；

所以我们要想办法引入隐变量将其对数似然函数变成一个二次多项式。
引入隐变量使用EM算法求解

{这里比较简单，隐含变量z的分布很容易推断出来}

Note: z1概率分布的解释：b(y1,p)二项分布，表示进行y1次试验，结果为z1的概率为p=(1-θ)/(2-θ)。

EM迭代求解

 

Note: E步中是使用上一步得到的已知的θ(i)值来求解z的期望，并用z的期望取代z。

EM迭代过程

 

[概率论与数理统计教程 2ed. - 茆诗松]
EM算法解高斯混合模型Gaussian Mixture Models

假设我们需要调查我们学校的男生和女生的身高分布。在校园里随便地活捉了100个男生和100个女生，他们共200个人（也就是200个身高的样本数据）。
高斯模型

你开始喊：“男的左边，女的右边，其他的站中间！”。然后你就先统计抽样得到的100个男生的身高。假设他们的身高是服从高斯分布N(u,∂)的。但是这个分布的均值u和方差∂2我们不知道，这两个参数就是我们要估计的。记作θ=[u, ∂]T。用数学的语言来说就是：在学校那么多男生（身高）中，我们独立地按照概率密度p(x|θ)（高斯分布N(u,∂)的形式）抽取100了个（身高），组成样本集X，我们想通过样本集X来估计出未知参数θ。抽到的样本集是X={x1,x2,…,xN}，其中xi表示抽到的第i个人的身高，这里抽到的样本个数N就是100。用MLE解决的公式：

通过抽取得到的那100个男生的身高和已知其身高服从高斯分布，我们通过最大化其似然函数，就可以得到了对应高斯分布的参数θ=[u, ∂]T了。对于女生的身高分布也可以用同样的方法得到。
高斯混合模型      

再回到例子本身，如果没有“男的左边，女的右边，其他的站中间！”这个步骤，或者说我抽到这200个人都戴了面具。那现在这200个人已经混到一起了，这时候，你从这200个人（的身高）里面随便给我指一个人（的身高），我都无法确定这个人（的身高）是男生（的身高）还是女生（的身高）。也就是说你不知道抽取的那200个人里面的每一个人到底是从男生的那个身高分布里面抽取的，还是女生的那个身高分布抽取的。用数学的语言就是，抽取得到的每个样本都不知道是从哪个分布抽取的。这个时候，对于每一个样本或者你抽取到的人，就有两个东西需要猜测或者估计的了，一是这个人是男的还是女的？二是男生和女生对应的身高的高斯分布的参数是多少？   只有当我们知道了哪些人属于同一个高斯分布的时候，我们才能够对这个分布的参数作出靠谱的预测。只有当我们对这两个分布的参数作出了准确的估计的时候，才能知道到底哪些人属于第一个分布，那些人属于第二个分布。

       为了解决这个你依赖我，我依赖你的循环依赖问题，总得有一方要先打破僵局，先随便整一个值出来，然后我再根据你的变化调整我的变化，然后如此迭代着不断互相推导，最终就会收敛到一个解。

EM解高斯混合模型

      在我们上面这个问题里面，我们是先随便猜一下男生（身高）的正态分布的参数：如均值和方差是多少。例如男生的均值是1米7，方差是0.1米（当然了，刚开始肯定没那么准），然后计算出每个人更可能属于第一个还是第二个正态分布中的（lz某个身高值代入男生和女生的高斯分布中，选择概率大的作为其所在分布），这个是属于Expectation一步。

      有了每个人的归属，已经大概地按上面的方法将这200个人分为男生和女生两部分，我们就可以根据之前说的最大似然那样，通过这些被大概分为男生的n个人来重新估计第一个分布的参数，女生的那个分布同样方法重新估计。这个是Maximization。

      然后，当我们更新了这两个分布的时候，每一个属于这两个分布的概率又变了，那么我们就再需要调整E步……如此往复，直到参数基本不再发生变化为止。

      这里把每个人（样本）的完整描述看做是三元组yi={xi,zi1,zi2}，其中，xi是第i个样本的观测值，对应身高，可观测值。zi1和zi2表示男生和女生这两个高斯分布中哪个被用来产生值xi，就是说这两个值标记这个人到底是男生还是女生（的身高分布产生的）。这两个值我们是不知道的，是隐含变量。确切的说，zij在xi由第j个高斯分布产生时值为1，否则为0。例如一个样本的观测值为1.8，然后他来自男生的那个高斯分布，那么我们可以将这个样本表示为{1.8, 1, 0}。如果zi1和zi2的值已知，也就是说每个人我已经标记为男生或者女生了，那么我们就可以利用上面说的最大似然算法来估计他们各自高斯分布的参数。但是它们未知，因此我们只能用EM算法。

       假设已经知道这个隐含变量了，那么直接按上面说的最大似然估计求解那个分布的参数就好了。我们可以先给这个分布弄一个初始值，然后求这个隐含变量的期望，当成是这个隐含变量的已知值，那么现在就可以用最大似然求解那个分布的参数了吧，那假设这个参数比之前的那个随机的参数要好，它更能表达真实的分布，那么我们再通过这个参数确定的分布去求这个隐含变量的期望，然后再最大化，得到另一个更优的参数，……迭代，就能得到一个皆大欢喜的结果了。

       迭代时咋知道新的参数的估计就比原来的好啊？这些问题在数学上是可以稳当的证明的。

[Mitchell的Machine Learning书中一个EM应用的例子]

[从最大似然到EM算法浅解]

皮皮blog

 

EM算法
Jensen不等式

凸函数：优化理论中，设f是定义域为实数的函数，如果对于所有的实数x，clip_image002，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的（clip_image004），那么f是凸函数。如果clip_image006或者clip_image008，那么称f是严格凸函数。 当f是（严格）凹函数当且仅当-f是（严格）凸函数。比如clip_image037是凹函数。

Note: lz发现好多地方凸函数定义可能是相反的，国内外定义好像也不一样。这里只要记得，正常的碗就是凸函数！

Jensen不等式：如果f是凸函数，X是随机变量，那么clip_image010。特别地，如果f是严格凸函数，那么clip_image012当且仅当clip_image014，也就是说X是常量。Jensen不等式应用于凹函数时，不等号方向反向，也就是clip_image021。

这里我们将clip_image016简写为clip_image018。
      clip_image019

      图中，实线f是凸函数，X是随机变量，有0.5的概率是a，有0.5的概率是b。（就像掷硬币一样）。X的期望值就是a和b的中值了，图中可以看到clip_image010[1]成立。
EM算法
EM问题描述

      给定的训练样本是clip_image023，样例间独立，我们想找到每个样例隐含的类别z，能使得p(x,z)最大。

p(x,z)的最大似然估计如下：

      clip_image024

      第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求θ一般比较困难，因为有隐藏变量z存在，但是一般确定了z后，求解就容易了。竟然不能直接最大化clip_image028，我们可以不断地建立clip_image030的下界（E步），然后优化下界（M步）[E步固定clip_image026[8]，优化clip_image107，M步固定clip_image107[1]优化clip_image026[9]]。

      对于每一个样例i，让clip_image032表示该样例隐含变量z的某种分布（这里隐含变量z的分布概率是未知的），clip_image032[1]满足的条件是clip_image034。（如果z是连续性的，那么clip_image032[2]是概率密度函数，需要将求和符号换做积分符号）。比如要将班上学生聚类，假设隐藏变量z是身高，那么就是连续的高斯分布。如果按照隐藏变量是男女，那么就是伯努利分布了。
建立clip_image028[5]的下界（E步）

可以由前面阐述的内容得到下面的公式：

      clip_image035

上述公式推导：

    （1）到（2）比较直接，就是分子分母同乘以一个相等的函数。（2）到（3）利用了Jensen不等式，考虑到clip_image037是凹函数（二阶导数小于0），而且

      clip_image038就是clip_image039(Zi的函数)的期望

 Lazy Statistician规则:

    设Y是随机变量X的函数clip_image041（g是连续函数），那么

      （1） X是离散型随机变量，它的分布律为clip_image043，k=1,2,…。若clip_image045绝对收敛，则有

      clip_image047

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

      clip_image053

      对应于上述问题，Y是clip_image039[1]，X是clip_image055，clip_image057是clip_image059，g是clip_image055[1]到clip_image039[2]的映射。这样解释了式子（2）中的期望

根据凹函数时的Jensen不等式：

      clip_image060可以得到（3）。

对于clip_image032[3]的选择，有多种可能，那种更好的？

      假设clip_image026[1]已经给定，那么clip_image028[2]的值就决定于clip_image057[1]和clip_image062了。我们可以通过调整这两个概率使下界不断上升，以逼近clip_image028[3]的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率能够等价于clip_image028[4]了。按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：

clip_image063，c为常数，不依赖于clip_image065。对此式子做进一步推导，我们知道clip_image067，那么也就有clip_image069，那么有下式：

      clip_image070

      至此，我们推出了在固定其他参数clip_image026[2]后，clip_image072的计算公式就是后验概率{里面主要包含要求的参数，要初始化}，解决了clip_image072[1]如何选择的问题（保证了在给定clip_image077[2]时，Jensen不等式中的等式成立，取得下界）。
优化下界（M步）

接下来的M步，就是在给定clip_image072[2]后，调整clip_image026[3]，去极大化clip_image028[6]的下界（在固定clip_image072[3]后，下界还可以调整的更大）。
EM算法的步骤

Note: 通过男女混合高斯模型来记！

初始化分布参数θ；

循环重复直到收敛 {

      （E步）对于每一个i，计算

                   #E步骤：根据参数初始值或上一次迭代的模型参数来计算出隐性变量的后验概率，其实就是隐性变量的期望。作为隐藏变量的现估计值：

                  clip_image074

      （M步）计算

                   #将似然函数最大化以获得新的参数值(log式中分母可省略)

                  clip_image075

那么究竟怎么确保EM收敛？

假定clip_image077和clip_image079是EM第t次和t+1次迭代后的结果。

如果我们证明clip_image081，也就是说极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。下面来证明，选定clip_image077[1]后，我们得到E步

      clip_image083

      这一步保证了在给定clip_image077[2]时，Jensen不等式中的等式成立，也就是

      clip_image084

      然后进行M步，固定clip_image086，并将clip_image088视作变量，对上面的clip_image090求导后，得到clip_image092，这样经过一些推导会有以下式子成立：

      clip_image093

      解释第（4）步，得到clip_image092[1]时，只是最大化clip_image090[1]，也就是clip_image095的下界，而没有使等式成立，等式成立只有是在固定clip_image026[4]，并按E步得到clip_image097时才能成立。

      况且根据我们前面得到的下式，对于所有的clip_image097[1]和clip_image026[5]都成立

      clip_image098

      第（5）步利用了M步的定义，M步就是将clip_image088[1]调整到clip_image100，使得下界最大化。因此（5）成立。（自证明？）

      （6）是之前的等式结果

      这样就证明了clip_image102会单调增加。一种收敛方法是clip_image102[1]不再变化，还有一种就是变化幅度很小。

      再次解释一下（4）、（5）、（6）。首先（4）对所有的参数都满足，而其等式成立条件只是在固定clip_image026[6]，并调整好Q时成立，而第（4）步只是固定Q，调整clip_image026[7]，不能保证等式一定成立。（4）到（5）就是M步的定义，（5）到（6）是前面E步所保证等式成立条件。也就是说E步会将下界拉到与clip_image102[2]一个特定值（这里clip_image088[2]）一样的高度，而此时发现下界仍然可以上升，因此经过M步后，下界又被拉升，但达不到与clip_image102[3]另外一个特定值一样的高度，之后E步又将下界拉到与这个特定值一样的高度，重复下去，直到最大值。

      如果我们定义

      clip_image103

      从前面的推导中我们知道clip_image105，EM可以看作是J的坐标上升法，E步固定clip_image026[8]，优化clip_image107，M步固定clip_image107[1]优化clip_image026[9]。
重新审视混合高斯模型

      我们已经知道了EM的精髓和推导过程，再次审视一下混合高斯模型。之前提到的混合高斯模型的参数clip_image109和clip_image111计算公式都是根据很多假定得出的，有些没有说明来由。为了简单，这里在M步只给出clip_image113和clip_image115的推导方法。

E步很简单，按照一般EM公式得到：

      clip_image116

      简单解释就是每个样例i的隐含类别clip_image055[2]为j的概率可以通过后验概率计算得到。

      在M步中，我们需要在固定clip_image072[4]后最大化最大似然估计，也就是

      clip_image118

      这是将clip_image120的k种情况展开后的样子，未知参数clip_image122和clip_image124。

      固定clip_image126和clip_image128，对clip_image130求导得

      clip_image131

      等于0时，得到之前模型中的clip_image115[1]的更新公式

      clip_image132

      然后推导clip_image126[1]的更新公式。看之前得到的

      clip_image133

      在clip_image113[1]和clip_image115[2]确定后，分子上面的一串都是常数了，实际上需要优化的公式是：

      clip_image134

      需要知道的是，clip_image126[2]还需要满足一定的约束条件就是clip_image136。

      这个优化问题我们很熟悉了，直接构造拉格朗日乘子。

      clip_image137

      还有一点就是clip_image139，但这一点会在得到的公式里自动满足。

      求导得，

      clip_image141

      等于0，得到

      clip_image142

      也就是说clip_image143再次使用clip_image136[1]，得到

      clip_image144

      这样就神奇地得到了clip_image146。

      那么就顺势得到M步中clip_image126[3]的更新公式：

      clip_image147

      clip_image111[1]的推导也类似，不过稍微复杂一些，毕竟是矩阵。结果在之前的混合高斯模型中已经给出。
总结

      如果将样本看作观察值，潜在类别看作是隐藏变量，那么聚类问题也就是参数估计问题，只不过聚类问题中参数分为隐含类别变量和其他参数，这犹如在x-y坐标系中找一个曲线的极值，然而曲线函数不能直接求导，因此什么梯度下降方法就不适用了。但固定一个变量后，另外一个可以通过求导得到，因此可以使用坐标上升法，一次固定一个变量，对另外的求极值，最后逐步逼近极值。对应到EM上，E步估计隐含变量，M步估计其他参数，交替将极值推向最大。EM中还有“硬”指定和“软”指定的概念，“软”指定看似更为合理，但计算量要大，“硬”指定在某些场合如K-means中更为实用（要是保持一个样本点到其他所有中心的概率，就会很麻烦）。

      另外，EM的收敛性证明方法确实很牛，能够利用log的凹函数性质，还能够想到利用创造下界，拉平函数下界，优化下界的方法来逐步逼近极大值。而且每一步迭代都能保证是单调的。最重要的是证明的数学公式非常精妙，硬是分子分母都乘以z的概率变成期望来套上Jensen不等式。
EM算法另一种理解

坐标上升法（Coordinate ascent）：

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

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

lz NG的协同过滤参数更新算法也有类似EM的算法，再研究。

from:http://blog.csdn.net/pipisorry/article/details/42550815

ref:http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html

Tom Mitchell的Machine Learning书

M. Bishop - PatternRecognitionAndMachineLearning: EM算法原理
---------------------  
作者：-柚子皮-  
来源：CSDN  
原文：https://blog.csdn.net/pipisorry/article/details/42550815  
版权声明：本文为博主原创文章，转载请附上博文链接！

 

 

　今天不太想学习，炒个冷饭，讲讲机器学习十大算法里有名的EM算法，文章里面有些个人理解，如有错漏，还请读者不吝赐教。

　　众所周知，极大似然估计是一种应用很广泛的参数估计方法。例如我手头有一些东北人的身高的数据，又知道身高的概率模型是高斯分布，那么利用极大化似然函数的方法可以估计出高斯分布的两个参数，均值和方差。这个方法基本上所有概率课本上都会讲，我这就不多说了，不清楚的请百度。

　　然而现在我面临的是这种情况，我手上的数据是四川人和东北人的身高合集，然而对于其中具体的每一个数据，并没有标定出它来自“东北人”还是“四川人”，我想如果把这个数据集的概率密度画出来，大约是这个样子：

　　

　　好了不要吐槽了，能画成这个样子我已经很用心了= =

　　其实这个双峰的概率密度函数是有模型的，称作高斯混合模型（GMM），写作：

　　

　　话说往博客上加公式真是费劲= =这模型很好理解，就是k个高斯模型加权组成，α是各高斯分布的权重，Θ是参数。对GMM模型的参数估计，就要用EM算法。更一般的讲，EM算法适用于带有隐变量的概率模型的估计，什么是隐变量呢？就是观测不到的变量，对于上面四川人和东北人的例子，对每一个身高而言，它来自四川还是东北，就是一个隐变量。

　　为什么要用EM，我们来具体考虑一下上面这个问题。如果使用极大似然估计——这是我们最开始最单纯的想法，那么我们需要极大化的似然函数应该是这个：

　　

　　然而我们并不知道p(x;θ)的表达式，有同学说我知道啊，不就是上面那个混个高斯模型？不就是参数多一点麽。

　　仔细想想，GMM里的θ可是由四川人和东北人两部分组成哟，假如你要估计四川人的身高均值，直接用GMM做似然函数，会把四川人和东北人全考虑进去，显然不合适。

　　另一个想法是考虑隐变量，如果我们已经知道哪些样本来自四川，哪些样本来自东北，那就好了。用Z=0或Z=1标记样本来自哪个总体，则Z就是隐变量，需要最大化的似然函数就变为：

　　然而并没有卵用，因为隐变量确实不知道。要估计一个样本是来自四川还是东北，我们就要有模型参数，要估计模型参数，我们首先要知道一个样本是来自四川或东北的可能性...

　　到底是鸡生蛋，还是蛋生鸡？

　　不闹了，我们的方法是假设。首先假设一个模型参数θ，然后每个样本来自四川/东北的概率p(zi)就能算出来了，p(xi,zi)=p(xi|zi)p(zi)，而x|z=0服从四川人分布，x|z=1服从东北人分布，所以似然函数可以写成含有θ的函数，极大化它我们可以得到一个新的θ。新的θ因为考虑了样本来自哪个分布，会比原来的更能反应数据规律。有了这个更好的θ我们再对每个样本重新计算它来自四川和东北的概率，用更好的θ算出来的概率会更准确，有了更准确的信息，我们可以继续像上面一样估计θ，自然而然这次得到的θ会比上一次更棒，如此蒸蒸日上，直到收敛（参数变动不明显了），理论上，EM算法就说完了。

　　然而事情并没有这么简单，上面的思想理论上可行，实践起来不成。主要是因为似然函数有“和的log”这一项，log里面是一个和的形式，一求导这画面不要太美，直接强来你要面对 “两个正态分布的概率密度函数相加”做分母，“两个正态分布分别求导再相加”做分子的分数形式。m个这玩意加起来令它等于0，要求出关于θ的解析解，你对自己的数学水平想的不要太高。

　　怎么办？先介绍一个不等式，叫Jensen不等式，是这样说的：

　　X是一个随机变量，f(X)是一个凸函数（二阶导数大或等于0），那么有：

　　

 

　　当且仅当X是常数的时候等号成立

　　如果f（X）是凹函数，不等号反向

　　关于这个不等式，我既不打算证明，也不打算说明，希望你承认它正确就好。

　　半路杀出一个Jensen不等式，要用它解决上面的困境也是应有之义，不然说它做什么。直接最大化似然函数做不到，那么如果我们能找到似然函数的一个紧的下界一直优化它，并保证每次迭代能够使总的似然函数一直增大，其实也是一样的。怎么说？画个图你就明白了：

　　

　　图画的不好，多见谅。横坐标是参数，纵坐标是似然函数，首先我们初始化一个θ1，根据它求似然函数一个紧的下界，也就是图中第一条黑短线，黑短线上的值虽然都小于似然函数的值，但至少有一点可以满足等号（所以称为紧下界），最大化小黑短线我们就hit到至少与似然函数刚好相等的位置，对应的横坐标就是我们的新的θ2，如此进行，只要保证随着θ的更新，每次最大化的小黑短线值都比上次的更大，那么算法收敛，最后就能最大化到似然函数的极大值处。

　　构造这个小黑短线，就要靠Jensen不等式。注意我们这里的log函数是个凹函数，所以我们使用的Jensen不等式的凹函数版本。根据Jensen函数，需要把log里面的东西写成一个数学期望的形式，注意到log里的和是关于隐变量Z的和，于是自然而然，这个数学期望一定是和Z有关，如果设Q(z)是Z的分布函数，那么可以这样构造：

　　这几句公式比较多，我不一一敲了，直接把我PPT里的内容截图过来：

　　

　　所以log里其实构造了一个随机变量Y，Y是Z的函数，Y取p/Q的值的概率是Q，这点说的很清楚了。

　　构造好数学期望，下一步根据Jensen不等式进行放缩：

　　

　　有了这一步，我们看一下整个式子：

　　

　　也就是说我们找到了似然函数的一个下界，那么优化它是否就可以呢？不是的，上面说了必须保证这个下界是紧的，也就是至少有点能使等号成立。由Jensen不等式，等式成立的条件是随机变量是常数，具体到这里，就是：

　　　　

　　又因为Q（z）是z的分布函数，所以：

　　把C乘过去，可得C就是p(xi,z)对z求和，所以我们终于知道了：

　　

　　得到Q(z)，大功告成，Q(z)就是p(zi|xi)，或者写成p(zi)，都是一回事，代表第i个数据是来自zi的概率。

　　于是EM算法出炉，它是这样做的：

　　首先，初始化参数θ

　　（1）E-Step：根据参数θ计算每个样本属于zi的概率，即这个身高来自四川或东北的概率，这个概率就是Q

　　（2）M-Step：根据计算得到的Q，求出含有θ的似然函数的下界并最大化它，得到新的参数θ

　　重复（1）和（2）直到收敛，可以看到，从思想上来说，和最开始没什么两样，只不过直接最大化似然函数不好做，曲线救国而已。

　　至于为什么这样的迭代会保证似然函数单调不减，即EM算法的收敛性证明，我就先不写了，以后有时间再考虑补。需要额外说明的是，EM算法在一般情况是收敛的，但是不保证收敛到全局最优，即有可能进入局部的最优。EM算法在混合高斯模型，隐马尔科夫模型中都有应用，是著名的数据挖掘十大算法之一。

　　就酱~有什么错漏和不同见解欢迎留言评论，我的推导和思路整理和网上的其他并不是完全相同，见仁见智吧~

 

注：

1.利用联合概率求最大后验概率，可以用似然估计，但是当式子较为复杂的时候可以使用EM算法，

2.EM算法是一种求下界的逼近算法，不断的使下界值最大，即可求出最大后验概率，

3.利用来jensen不等式，和概率里面的lazy规则