EM算法学习(一)
EM算法是英文expectation-maximization算法的英文简写,翻译过来就是期望最大化算法,其实是一种根据求参的极大似然估计的一种迭代的优化策略,EM算法可以广泛估计是因为他可以从非完整的数据集中对于参数进行极大似然的估计,这样的方法对于处理残缺数据,截尾数据和一些带有噪声的数据来说是很有效的.
在写这篇文章之前,我看了很多篇博客,学习了很多的知识,也参照了很多的资料,希望可以从EM算法的迭代优化理论和一般的步骤中出发,然后能够举一个例子来使我们理解这个EM算法,然后在对其收敛性进行证明,目的是为了说明EM算法每一次迭代都是能够提高似然函数值然后收敛到一个稳定的点,再引出EM算法的收敛速度.
大概通过上述部分,我们可以得到基于其简单,收敛,稳定上升的优势,但是也会产生一些缺点,比如收敛速度过慢的加速方法等,在第二篇文章中将会介绍这个处理缺点的方法,然后会写一些关于EM算法的重要应用,包括EM算法在二元正态分布上的参数估计的应用,混合高斯分布参数估计方面的应用,以及EM算法在隐马尔科夫模型上参数的应用(一种EM算法的特殊情形),希望通过这一系列的文章可以让大家理解好EM算法的明显优势以及原理,让我们开始吧!
背景:
极大似然估计和贝叶斯统计其实是作为现在的统计领域中非常热门的领域了,其实来说他们的计算过程是有一定的相似成分的,比如极大似然函数估计在计算的方法上跟贝叶斯的后验概率的计算是非常相似的,学过统计学习的我们知道,贝叶斯是分为两种的大类的,一种是拥有显式的后验分布,这样的一般用于简单的似然函数,另外一种是数据添加的算法,有些时候我们的数据可能会存在缺失或者是似然函数不是显性的,数据添加类在这时候就可以很好的应用,他可以将已经观测到的数据基础上加上一些”潜在数据”,从而使得变得更简单,完成极大化的工作,然后我们常用的一种数据添加法其实就是我们今天介绍的EM算法.
EM算法是一种迭代的优化策略,他的计算方法是分为期望步(E步)和极大步(M步)的,所以这个算法的名字是这样来的,EM算法受到了缺失算法的影响,最初就是为了解决上边提到的数据缺失的问题,基本的思想就是首先根据已经观测出来的数据估计出模型参数的值,然后再根据上一步估计出的参数值来估计缺失数据的值,然后再根据估计中缺失的数据加上之前的已经观测到的数据重新在对参数值进行估计,然后反复的进行迭代,直到最后收敛,迭代结束.
而现在EM算法发展了几十年了,在当时的数据快速增长得那个时代,那时候处理数据很困难,经常会出现数据缺失或者不可用的情况,当时无非就是用用神经网络拟合,添补法,卡尔曼滤波法等等,但是最后还是EM脱颖而出,最主要还是他的算法步骤简单,稳定上升可以很可靠的找到最优的收敛值,但是运用这种思想,我们拓展到了简化问题策略,有时候缺失数据并非真的缺少了,这时候EM引入恰当的数据添加技术,这样的数据被称为”潜在数据”,复杂问题通过引入潜在数据,能够有效的解决我们的问题
“潜在数据”可以解释为数据本身并不存在缺失变量,但观察数据比较难以处理,如果添加上额外的变量,处理起来会变得比较简 单。假设X是已知的观测数据,想象由随机变量X生成的观察数据连同来 自随机变量y的缺失或未观测数据,得到Z=(X,Y)为完全数据。通过给定观察数据X,我们希望最大化似然的函数L(0/x).由于数据缺失或者其他原因导
致采用该似然函数会难以处理,而采用Z|0和Y|(x,0)的密度则比较容易处 理。EM算法通过采用这些较容易的密度p(0|z),从而避开考虑了P(0|X).但是这在贝叶斯应用中,对于后验分布的P都是随机变量.
但是不可避免EM算法也有一些缺点:
1:在缺失数据较多的情形,收敛的速度较慢.
2:对于某些情况下,要计算算法中的M步,即完成对似然函数的估计是非常困难的
3:在某些情况下是要获得EM算法中的E步的期望显式是非常困难或者不可能的
算法原理和步骤:
现在我们假设X是观测数据,Y是潜在数据,EM算法迭代是为了寻求关于0最大化函数L(0|X),设0(k)是在进行K次迭代以后估计得到的最大值点,K属于(0,1,2......),定义Q(0|0(k))
是在观测数据X ={x1,x2….xn}的条件下完全数据的联合对数函数似然的期望,既可以获得以下式子:
![60f08a6e7be42c9c24d1a050eaa8c684.png](https://i-blog.csdnimg.cn/blog_migrate/49ec133f20bf7325825c8abbf93f84db.jpeg)
通过上式我们可以发现,一旦我们的样本点X给定,那么Y就是Z的唯一的随机部分,通过对于Y求条件期望,就又可以把Y给积分掉了,这样使得Q(0|0(k))就完全成为了一个关于0的函数,这样就可以求使得Q(0|0(k))最大的0,并且记为0(k+1),以供下一次迭代使用.
EM算法从0(0)开始,然后在这两部之间进行交替,E表示期望,M表示最大化,该算法可以概括如下:
1:E步,在给定的观测数据X和已经知道的参数条件下,求缺失数据Y的数学期望,即计算上面的提到的对数似然函数的条件期望Q(0|0(k)).
2:M步,就像不存在缺失数据Y一样(在填充缺失数据后),针对完全数据下的对数似然函数的期望进行最大化处理,即求关于0似然函数Q(0|0(k))的最大化.设0(k+1)等于最大值点,更新0(k):
![7f299e0a5264f79e15b1bd79f1725f3f.png](https://i-blog.csdnimg.cn/blog_migrate/1bfedc53a9cfddedbd62098f5513d977.jpeg)
3:返回E步,知道满足某停止规则为止.一般依赖于
![22239d7d25ff974fa84cdc4ba2add0e2.png](https://i-blog.csdnimg.cn/blog_migrate/dc1aab5da46cd5c5b6c4808e44aa771f.jpeg)
现在使用韩佳伟的数据分析中的例子来详细的说明EM算法:
1:现在假设进行一个实验会出现现在四种结果,每种结果发生的概率分别为:
![8dd6129417ba9702125754ca2aea9dfa.png](https://i-blog.csdnimg.cn/blog_migrate/c1ee2c817d7694d384012c4100237859.jpeg)
实验进行了197次,四种结果发生了125,18,20,34次,这个时候的X={x1,x2,x3,x4}={125,18,20,34}
为了估计参数,我们可以取0的先验分布p(0)为U(0,1),由贝叶斯公式可知,0的后验分布为:
![ee1b4c36fdabd9dfd8a492f1a8e5f17c.png](https://i-blog.csdnimg.cn/blog_migrate/65108d7a318f1cdac9e0b72ec117f63b.jpeg)
把第一种结果分成发生概率分别为1/2和0/4的两个部分,令Y和x1-Y分别表示为这两部分试验成功的次数(Y为缺失数据),则0的添加的后验分布为:
![4db67f9e9bd579636b020759a8bd7df5.png](https://i-blog.csdnimg.cn/blog_migrate/465791a09f44a17a89966b5911b1ac02.jpeg)
这时候就该轮到了EM算法添加数据了,直接求0的极大似然估计也是比较麻烦的,现在使用EM算法后迭代最后一个后验分布函数就简单多了,(在上面的计算过程中,最下边的那个符号,他表示的是符号两端的式子成比例,并且这个比例跟0无关,这个比例不会影响到EM迭代算法估计的结果,因为后面的极大化可以进行约去).
现在我们假设在第i+1次的迭代中,有估计值0(i),则可以通过EM算法中的E步和M步得到一个新的估计,在E步中:
![e0c4d669d6ddd494a449efb2fe249513.png](https://i-blog.csdnimg.cn/blog_migrate/f0802551c6ac39d0c74c10e2b14895be.jpeg)
因为在x和0(i)给定的情况下,Y服从二项式分布,因此可以得到E(Y|X,0(i))=2x(1)/[2+0(i)],于是便有:
![a603d2af638f5533bc2dfc884abbc919.png](https://i-blog.csdnimg.cn/blog_migrate/35ef37688227ed560e2d38b960eef672.jpeg)
在M步中,对上边的式子中的0进行求导并且使结果为0,可以得到迭代的公式为:
![063028a5d2ec86c0a64b452348a51fee.png](https://i-blog.csdnimg.cn/blog_migrate/873c5268c3d2e54e1b1f47c35eb33832.jpeg)
这个时候从0(0)=0.5开始,经过EM四次迭代以后收敛到0.6268,根据Matlab进行作图,收敛曲线如下所示:
![64a6bd37fb2beee9e146c88c3c64dd7b.png](https://i-blog.csdnimg.cn/blog_migrate/964e43befd71e97656fbf806aadbb59c.jpeg)
另外韩佳伟那本书使用的R,我这里使用的是MATLAB,效果类似,都可以实现.
EM算法收敛性证明:
算法简单、收敛稳定是EM算法的最大优点,EM算法的简便性在上文中已经得到充足的认识,现在需要对其估计得到的序列是不是达到了预期 的收敛要求,即EM算法估计的结果是不是能稳定收敛,以及收敛结果是不 是p(0IX)的最大值或局部最大值来进行验证,本文下面需要来证明EM算法的收敛性。
第一步:吉布斯不等式:
![d0a396331970d6d6fb32ff91e03754aa.png](https://i-blog.csdnimg.cn/blog_migrate/b0560fa95ac09e116817a99335c40040.jpeg)
当且仅当x=1时,等号成立
{PS:其实吉布斯不等式是詹森不等式的一个特例,因为log(x)是一个凸函数,作图的话就能够发现在某一个时刻,吉布斯不等式和詹森不等式是有交点的在曲线上}
第二步:
现在先引出几个定理:如果f(x),g(x)具有相同的支集的概率密度函数,现在令:
![1678b313add3979c67202dffacefd680.png](https://i-blog.csdnimg.cn/blog_migrate/72754b0ade0e778308026d0a5a6aebe5.jpeg)
则:
![463c9a7d57d68455a7b58a281bf16cbe.png](https://i-blog.csdnimg.cn/blog_migrate/2408bc1984037a156ee351c4bbe3bb3a.jpeg)
证明:
![68b37431c6df0dd6df9a1f1f94e75ca5.png](https://i-blog.csdnimg.cn/blog_migrate/8c939a0917db2c6a7eaf4a34ff22984c.jpeg)
当且仅当f(x)=g(x)时,等号成立,此时log(g(x)/f(x))=f(x)/g(x)-1
第三步:
记录EM算法的估计序列为{0(k)},证明EM算法每个最大化步都提高了对于观测数据的对数似然函数L(0|X)的值,即:
![6a3baead531778d860fe9a660191d1ca.png](https://i-blog.csdnimg.cn/blog_migrate/3615499959f3aa7d7e5c3450602cea2e.jpeg)
记:
![726b91affea990019eb4eb7f1d1e924a.png](https://i-blog.csdnimg.cn/blog_migrate/feecfbf143f3ee62f8e640406bf4dfbc.jpeg)
证明:已知:
![5b0b4b3c75e5617af39ac0e440aa7b8c.png](https://i-blog.csdnimg.cn/blog_migrate/3ad95d7bba3d8ee32c1602dec5c9623e.jpeg)
然后根据贝叶斯统计先验后验的函数关系式得到:
![b6f4f7175ee6139f6f0dc23721eccbe0.png](https://i-blog.csdnimg.cn/blog_migrate/ba8374c1a484e1ed0b3cde39d4d09de4.jpeg)
可以得知,在大样本以及0均匀分布的前提下,其中p(x|0)是已经观测到的X的密度,而p(y|x,0)是给定已经观测到的数据后缺失数据的密度:
![82afa531d4afdab740206c26252b4505.png](https://i-blog.csdnimg.cn/blog_migrate/3811326b1c5480076e49ca891dc07c57.jpeg)
记:
![d0297b1dec7c6c5f1862db5406ddca06.png](https://i-blog.csdnimg.cn/blog_migrate/d9bd4fa42322e245bb12e62f577d5d23.jpeg)
则:
![1f7eccdea4d37d3d25d37b8e02ede14c.png](https://i-blog.csdnimg.cn/blog_migrate/3c2f5f9be44591a9941845929144d42b.jpeg)
从而:
![a251aaff5f5311a3ec159ad7dfebb295.png](https://i-blog.csdnimg.cn/blog_migrate/3dba81a455b6136d3e5243fe6b8ffd24.jpeg)
进而得到:
![864697a94fa9bdcfb5bd8862417f13ce.png](https://i-blog.csdnimg.cn/blog_migrate/f31311afda6d7965653742294279176f.jpeg)
所以:
![d37e0f9d414875a7827e98a3440a5742.png](https://i-blog.csdnimg.cn/blog_migrate/1acacfdc77f7fa3b5e88e96eebcf80ff.jpeg)
由M步一直求极大化过程:
![7b53e28c19e293c4abffb8a38c10cdd5.png](https://i-blog.csdnimg.cn/blog_migrate/2357ae07bd8adbaab5eaeb2883b95ea0.jpeg)
可知Q一直在增大,有第二步的定理1得到:
![46bf5d952e9c423b5262ac627881c709.png](https://i-blog.csdnimg.cn/blog_migrate/31f2530458b9a20cd413a925104f2634.jpeg)
这里可以知道在EM算法在E步和M步的交替运算下,似然函数是不断地变大的,我们就可以认为估计的参数序列最终会收敛到极大似然估计.
如果在每次迭代中,都是通过求似然函数的极大似然估计,选择最大化的0(k+1)来代替0,这样就构成了EM算法,大部分时候极大似然函数都是有显式表达式,但是不是总是这样,所以有时候会有广义EM算法(GEM),增大Q的哪一步也增大了对数似然函数.
这里我不加证明的给出,得到的收敛性结论主要是针对对数似然函数值给出的,而不是针对的估计序列的收敛性;而且在一般的情况下,我们用EM算法得到的估计值0(k),只能保证收敛到似然函数的一个稳定点,并不能其保证收敛到全局最大值点或者局部最大值点。
参考资料:
The EM algorithm 吴恩达
http://blog.csdn.net/zouxy09/article/details/8537620/
http://www.cnblogs.com/openeim/p/3921835.html
以及更多的人的帮助,谢谢他们,接下来将会进行后续的对于加速EM算法收敛的方法,感谢阅读参考文献:K码农-http://kmanong.top/kmn/qxw/form/home?top_cate=28