EM算法理論及其應用

1、2009年11月第15卷第4期安慶師范學院學報(自然科學版Journal of Anqing Teachers College(Natural Science EditionNov.2009楊基棟(華東師范大學金融與統(tǒng)計學院,上海200062摘要:EM算法是一種迭代算法,主要用來計算后驗分布的眾數或極大似然估計,廣泛地應用于缺損數據、截尾數據、成群數據、帶有討厭參數的數據等所謂的不完全數據的統(tǒng)計推斷問題。在介紹EM算法的基礎上,針對EM算法收斂速度慢的缺陷,具體討論了加速EM算法:EMB算法和MEMB算法;針對EM算法計算的局限性,給出了EM算法的推廣:GEM和MCEM算法。最后給出了EM的實

2、值實例,結果精確。關鍵詞:EM算法;極大似然估計;GEM算法;MCEM算法;EMB算法;MEMB算法中圖分類號:O212.8文獻標識碼:A文章編號:1007-4260(200904-0030-060引言在統(tǒng)計領域里,統(tǒng)計計算技術近年來發(fā)展很快,它使許多統(tǒng)計方法,尤其是Bayes統(tǒng)計得到廣泛的運用。Bayes計算方法有很多,大體上可分為兩大類:一類是直接應用于后驗分布以得到后驗均值或后驗眾數的估計,以及這種估計的漸進方差或其近似;另一類算法可以總稱為數據添加算法,這是近年發(fā)展很快而且應用很廣的一種算法,它是在觀測數據的基礎上加一些“潛在數據”,從而簡化計算并完成一系列簡單的極大化或模擬,該“潛在

3、數據”可以是“缺損數據”或未知參數。其原理可以表述如下:設我們能觀測到的數據是Y,關于Y的后驗分布p(|Y很復雜,難以直接進行各種統(tǒng)計計算,假如我們能假定一些沒有能觀測到的潛在數據Z為已知(譬如,Y為某變量的截尾觀測值,則Z為該變量的真值,則可能得到一個關于的簡單的添加后驗分布p(|Y,Z,利用p(|Y,Z的簡單性我們可以對Z的假定作檢查和改進,如此進行,我們就將一個復雜的極大化或抽樣問題轉變?yōu)橐幌盗泻唵蔚臉O大化或抽樣。EM算法就是一種常用的數據添加算法。1EM算法及其理論先考慮一個簡單情形。設某元件的失效時間Y關于某個變量x有直線回歸關系,假設在一次試驗中得到一批觀測數據,見右圖,“

4、5;”表示該種元件的失 效時間對應的值,“”對應元件的(右截尾時間(比實際失效時間要小。如果直線斜率和截距的估計值已知,則我們可以在真實數據不小于截尾數據的前提下,將各個被截尾的失效時間估計出來(譬如,若E(Y|x>Z,則用E(Y|x作為真實數據,否則取Z作為真實數據,從而得到所謂的“完全數據”,由此“完全數據”,我們可以對直線斜率和截距進行估計(如極大似然估計,估計出新的斜率和截距后,在重新估計各個被截尾的失效時間,得到新的完全數據,如此重復,我們將一個復雜的估計時間替換成一系列簡單的估計問題。將之一般化,我們可以給出EM算法。EM算法是一種迭代方法,最初由Demp ster等于197

5、7年首次提出,主要用來計算后驗分布的眾數或極大似然估計。近十年來引起了統(tǒng)計學家們的極大興趣,在統(tǒng)計領域得到廣泛應用。這種方法可以廣3收稿日期:2009-05-21作者簡介:楊基棟,女,安徽安慶人,華東師范大學金融與統(tǒng)計學院助理工程師。泛的應用于缺損數據,截尾數據,成群數據,帶有討厭參數的數據等所謂的不完全數據。它的每一次迭代有兩步組成:E 步(求期望和M 步(極大化。一般的,以p (|Y 表示的基于觀測數據的后驗分布密度函數,稱為觀測后驗分布,p (|Y ,Z 表示添加數據Z 后得到的關于的后驗分布密度函數,稱為添加后驗分布,p (Z |,Y 表示在給定和觀測數據Y 下潛在數據Z 的條件分布密

6、度函數。我們的目的是計算觀測后驗分布p (|Y 的眾數,于是,EM 算法如下進行。記(i 為第i +1次迭代開始時后驗眾數的估計值,則第i +1次迭代的兩步為E 步:將p (|Y ,Z lo g p (|Y ,Z 后關于Z 的條件分布求期望,從而把Z 積掉,即Q (|(i ,Y E z lo g p (|Y ,Z |(i ,Y =log p (|Y ,Z p (Z |(i ,Y d Z (1M 步:將Q (|(i ,Y 極大化,即找一個點(i+1,使Q (i+1|(i ,Y =max Q (|(i ,Y (2如此形成了一次迭代(i (i+1。將上述E 步和M 步進行迭代直至(i+1-(i 或Q

7、 (i+1|(i ,Y -Q (i |(i ,Y 充分小時停止。例1假設一次試驗可能有四個結果,其發(fā)生的概率分別為125,18,20,34。此處觀測數據為Y =(y 1,y 2,y 3,y 4=(125,18,20,34,取的先驗分布(為(0,1上均勻分布,則的觀測后驗分布為p (|Y (p (Y |=(12+4y 114(1-y 214(1-y 3(4y 4(2+y 1(1-y 2+y 3y 4(3假設第一種結果可以分解成兩部分,其發(fā)生概率分別為12和4,令Z 和y 1-Z 分別表示試驗中結果落入這兩部分的次數(Z 是不能觀測的潛在數據,則的添加后驗分布為p (|Y ,Z (p (Y ,Z

8、|=(12z (4y 1-z 14(1-y 214(1-y 3(4y 4y 1-z +y 4(1-y 2+y 3(4用(3求的后驗眾數比較麻煩,而用(4求后驗眾數則十分簡單。在第i +1次迭代中,假設有估計值(i ,則可通過E 步和M 步得到一個新的估計,在E 步中有:Q (|(i ,Y =E z (y 1-Z +y 4log +(y 2+y 3lo g (1-|(i ,Y =y 1-E z (Z |(i ,Y +y 4lo g+(y 2+y 3log (1-因在(i 和Y 給定下,Z b (y 1,2/(i +2,故E z (Z |(i ,Y =2y 1/(i +2在M 步中,我們將Q (|

9、(i ,Y 對求導并令其為0,有(i+1=y 1+y 4-E z (Z |(i ,Y y 1+y 2+y 3+y 4-E z (Z |(i ,Y =(i y 1+(i +2y 4(i y 1+(i +2(y 2+y 3+y 4=159(i +68197(i +144(5(5式給出了有EM 算法得到的迭代公式。由此公式進行迭代,可得到觀測后驗分布的眾數。由(5式不用迭代也可求出的估計。事實上,它是一個迭代公式,若收斂到=(159+68/(197+144,解之有=0.626821497。2EM 算法的收斂性2.1EM 算法的收斂性定理EM 算法的最大優(yōu)點是簡單和穩(wěn)定。EM 算法的主要目的是提供一個

10、簡單的迭代算法來計算極大似然估計,人們自然會問,如此建立的EM 算法是否達到預期的要求,就是說,由EM 算法得到的估計序列是否收斂,如果收斂,其結果是p (|Y 的最大值或局部最大值。為此我們給出以下兩個定理。記EM 算法得到的估計序列為(i ,i =1,2,(|Y =lo g p (|Y 。定理1EM 算法在每一次迭代后均提高(觀測后驗密度函數值,即p (i+1|Y p (i |Y (6證明由全概率公式p (,Z |Y =p (Z |,Y p (|Y =p (|Y ,Z p (Z |Y 將上式后面兩項取對數有l(wèi)o g p (|Y =log p (|Y ,Z -log p (Z |,Y +lo

11、g p (Z |Y (7設現有估計(i ,將上式對Z 關于p (Z |(i ,Y 求期望,有l(wèi)og p (|Y =lo g p (|Y ,Z -log p (Z |,Y +log p (Z |Y p (Z |(i ,Y d Z Q (|(i ,Y -H (|(i ,Y +K (i ,Y (8其中Q (|(i ,Y 已在(1中定義,H (|(i ,Y =log p (Z |,Y p (Z |(i ,Y d Z ,K (i ,Y =log p (Z |Y p (Z |(i ,Y d Z13第4期楊基棟:EM 算法理論及其應用分別在(8中取為(i 和(i+1并相減,有l(wèi)og p (i+1|Y -lo

12、 g p (i |Y =Q (i+1|(i ,Y -Q (i |(i ,Y -H (i+1|(i ,Y -H (i |(i ,Y 由J ensen 不等式,EZ|(i ,y log (p (Z |(i ,Y p (Z |(i ,Y log E Z|(i ,y (p (Z |(i+1,Y p (Z |(i ,Y =0,故H (i+1|(i ,Y -H (i |(i ,Y 0,而(i+1是使Q (|(i ,Y 達到最大的,顯然Q (i+1|(i ,Y -Q (i |(i ,Y 0。這就說明,用EM 算法求出的(i 的確使其對數似然L (i |Y 關于i 單增,但只有在嚴格單增的情況下,才能保證EM

13、 算法求出的(×為的極大似然估計。這是從純數學的,實際應用有一定的困難。定理2(1如果p (|Y 有上界,則L (i |Y 收斂到某個L ×。(2如果Q (|關于和都連續(xù),則在關于L 的很一般的條件下,由EM 算法得到的估計序列(i 的收斂值×是L 的穩(wěn)定點。證明(1的證明由單調收斂定理立得;(2的證明見文獻4。關于定理2,我們指出,定理的條件在大多數場合是滿足的,定理的收斂性結論是針對(對數后驗密度函數數值給出的,而后驗密度函數值序列的收斂性比估計序列本身的收斂性更具意義。另外,在定理2的條件下,EM 算法的結果只能保證收斂到后驗密度函數的穩(wěn)定點,并不能保證收斂

14、到極大值點,事實上,任何一種算法都很難保證其結果為極大值點。較可行的辦法是選取幾個不同的初值進行迭代,然后在諸估計間加以選擇,這可以減輕初值選取對結果的影響。2.2EM 算法收斂速度的探討EM 算法是一種求參數極大似然估計的迭代算法,在處理不完全數據中有重要應用。EM 算法實現簡單,數值計算穩(wěn)定,存儲量小,并具有良好的全局收斂性。但是,EM 算法收斂速度相當慢,只是次線性的收斂速度,這個缺點防礙了EM 算法的應用?，F已提出了多種加速EM 算法收斂的方法,其中使用非線性規(guī)劃中的Broyden 對稱秩1校正公式,提出了一種加速EM 算法收斂的方法。與其他加速收斂方法相比,本方法簡便易行,不必對不完

15、全數據的似然函數一維搜索,在收斂性方面,與EM 算法一樣具有全局收斂性。數值試驗結果表明,本方法的收斂速度比EM 算法快的多。一般地,在EM 算法的M 步求(i+1時,可求解方程組5Q (|(i 5=0,得到參數的迭代公式(i+1=G (i 。這種迭代公式通常使EM 算法的收斂速度很慢,為加速收斂,可考慮使用其他方法求解。根據著名的Fisher 公式5Q (|(i 5|=(i =g (i ,這里g (i 是不完全數據對數似然函數L (在(i 處的梯度g (i =L (|=(i ,于是問題轉化為求(×,使g (×=0,從而可以使用非線性規(guī)劃中的有效方法求解,達到加速收斂的目的

16、。在非線性規(guī)劃的求解方法中,Broyden 對稱秩1校正公式具有二次終止性和超線性的收斂速度,使用它可望改善EM 步長的緩慢變化,以加速EM 算法的收斂。但是,如果不對L (進行一維搜索,Broyden 對稱秩1校正公式僅具有局部收斂性,因此考慮將Broyden 對稱秩1校正公式的二次終止性和EM 算法良好的全局收斂性相結合,形成一種能夠加速EM 算法收斂的新方法EM B 算法。在迭代初期,使用EM 算法,而當接近最優(yōu)解時,EM 步長的變化極為緩慢,這時使用Broyden 對稱秩1校正公式進行校正。算法1(EMB 算法1取初始值(0,常數(0,1,H (0=I ,r (0=g (0,k =0。

17、2若g (k =0,停止計算,(k 即為所求;否則轉3。3若g (k <r (k ,則取r (k+1滿足g (k ÷r (k+1<r (k ,置(k+1=(k -H (k g (k ,轉5,否則轉4。4置r (k+1=r (k ,取EM 步長(k ,置(k+1=(k +(k 。5對H (k 使用Broyden 對稱秩1校正公式,置H (k+1=H (k +(k -H (k (k (k -H (k (k T (k -H (k (k T (k ,其中k =g (k+1-g (k ,(k =(k+1-(k 。6置k =k +1,轉2。關于EMB 算法的收斂性有如下結論:定理3在

18、EM 算法收斂的條件下,EMB 算法是收斂的。證明分三種情形討論23安慶師范學院學報(自然科學版2009年1若在EMB 算法中,每次迭代都取EM 步長(k 計算(k+1,即EMB 算法即為EM 算法。2若在EMB 算法中,每次迭代都按公式(k+1=(k -H (k g (k 計算(k+1,則由g (k ÷r (k+1<r (k 可得r (k+1<r (k <2r (k-1<<k+1r (0,因<1,故k 時,記(k (×,則g (×=0,因此EMB 算法收斂。3一般情形時,若設EMB 算法局部收斂,即存在某個>0及k 1,

19、當k >k 1時,g (k >,而由非負單調遞減數列r (k 的構造可知,存在子列r (k ,r (k 0。因此存在k 2,當k >max k 1,k 2時,r (k <,于是g (k <r (k <(t 1,矛盾。因此EMB 算法收斂。從數值試驗結果來看,EMB 算法能顯著減少迭代次數,但有些時候效果不好,為此考慮在EMB 算法中加進一個延遲參數,目的是將校正盡量靠后,當接近最優(yōu)解時,再用Broyden 對稱秩1校正公式。算法2(M EMB 算法1取初始值(0,常數(0,1,非負整數,r (0=g (0,k =0。2若g (k =0,停止計算,(k 即為所

20、求;否則轉3。3若k >,轉4,否則轉5。4若g (k <r (k ,則取r (k+1滿足g (k ÷r (k+1<r (k ,置(k+1=(k -H (k g (k ,轉6,否則轉5。5置r (k+1=r (k ,取EM 步長(k ,置(k+1=(k +(k 。6對H (k 使用Broyden 對稱秩1校正公式,置H (k+1=H (k +(k -H (k (k (k -H (k (k T (k -H (k (k T (k ,其中k =g (k+1-g (k ,(k =(k+1-(k 。7置k =k +1,轉2。2.3EM 算法的推廣算法EM 算法得到廣泛應用的一

21、個重要原因是在M 步中,求極大化的方法與完全數據下求極大化的方法完全一樣。在許多場合,這樣的極大化有顯式表示,然而并不總是這樣,有時要找Q (|(i ,Y 達到最大的是困難的,一個比較簡單的方法是找一個(i+1,使Q (i+1|(i ,Y >Q (i |(i ,Y (9由(1和(9組成的EM 算法稱為廣義EM 算法(GEM 算法。我們可以利用無約束最優(yōu)化方法中的一些方法,像最速下降法,Newto n 法,共軛方向法和共軛梯度法,擬Newto n 法,Powell 方向加速法,來構造廣義EM 算法(GEM 算法,既可避免EM 算法M 步中表達式的局限,又可加速EM 算法,使EM 算法的收斂

22、速度得到很大的提高。在EM 算法的M 步求(i+1時,根據著名的Fisher 公式5Q (|(i 5|=(i =g (i ,這里g (i 是不完全數據對數似然函數L (在(i 處的梯度,于是問題轉化為求(3,使g (3=0,從而可以使用非線性規(guī)劃中的有效方法求解,達到加速收斂的目的。像上一章介紹的EMB 算法和M EMB 算法,就是兩種廣義的EM 算法,及GEM 算法。還有Laird 等利用數值分析中的外推Ait ken 加速來加速EM 算法,Lo uis 在Newto n 法的基礎上關于EM 算法提出了一種加速的方法,Amshidian 等提出了EM 算法的廣義共軛梯度加速方法等。而對于EM

23、 算法的E 步,有時要獲得期望的顯式表示是不可能的,即使近似計算也很困難,這時用Mo nte Carlo 方法來完成,就是所謂的Mo nte Carlo EM (MCEM 方法,它將E 步改為:(E1由p (Z |(i ,Y 隨機的抽取m 個隨機數Z 1,Z m ;(E2計算Q (|(i ,Y =1m m j =1log p (|Z j ,Y 。由大數定理,只要m 足夠大,Q (|(i ,Y 與(|(i ,Y 很接近,從而我們可以在M 步中對Q (|(i ,Y 求極大化。在MCEM 算法中有兩點是主要考慮的:一是m 的大小的確定,從精確角度來講,m 自然越大越好,但過大的m 使得計算的效率太低

24、,一般在開始時m 不需要很大;另一點是收斂性的判斷,因在E 步中是采用的Mo nte Carlo 方法,若要求這樣得到的(i 收斂到一點顯然是不現實的。在MCEM 中,收斂性的判別往往可借助圖形來進行。若經過若干次迭代后,迭代值圍繞直線=3小幅波動,則可以認為算法收33第4期楊基棟:EM 算法理論及其應用斂了,此時,為增加估計精度,可增加m 的值再運行一段時間,就可停止。MCEM 算法比較靈活,但是需要仔細選擇模擬容量和確保正確的收斂性準則。我們可以通過增加迭代次數來提高模擬容量,例如,在McColloch (1994方法中,模擬容量隨著迭代次數線性增加。除此以外,由于蒙特卡羅誤差,該EM 算

25、法不具有單調性,難以估計其收斂性。在McColloch 方法中,MCEM 算法通過規(guī)定迭代的次數來結束算法。Shi 和Lee (2000用臨時抽樣的方法直觀的觀察MCEM 算法的收斂性。3EM 算法的應用實例EM 算法可以應用于醫(yī)學研究中,尤其是臨床醫(yī)學中十分常見的一種數據觀測形式為重復觀測,其特點是在同一實驗單位上進行多次重復觀測,這個過程由于各種原因經常導致實驗觀測數據缺失,如動物的意外死亡,記錄儀器發(fā)生故障,被調查者拒絕回答相關調查項目等,然而,目前絕大多數重復觀測數據統(tǒng)計分析方法都有一個基本假定,即每個實驗單位具有完全的重復觀測數據(通常假定服從多元正態(tài)分布。因此有必要對缺失數據進行適

26、當的處理,否則數據分析時,通常的做法是刪去具有缺失的部分觀察記錄而不考慮記錄數據所蘊涵的信息,從而造成信息的損失以及分析結果的偏性。缺失數據的處理方式在很大程度上依賴醫(yī)學數據觀測的實際背景,例如,常規(guī)的實驗設計得到的獨立結構數據方差分析前的缺項估計的原理,都是解誤差離均差平方和對缺項變量的導數等于0的方程,從而獲得缺項的估計值。而醫(yī)學重復觀測設計得到的數據有別于獨立結構數據,其數據間存在一定的相關性,因此需要進行專門的缺項估計方法的研究。為此。我們在軟件編程的基礎上,應用EM 算法進行了這方面的數據處理。目前,可用于進行缺項估計的方法很多,但EM 算法較為優(yōu)秀,它是利用所有的資料信息來進行缺項

27、估計,以“補缺”的方式將含有缺失值的“不完全資料”轉化為“完全資料”,因此,對EM 算法得到的“完全資料”可以應用所有的重復觀測資料分析方法進行處理分析。EM 方法補缺使估計系數的方差明顯變小,從而提高了生長曲線系數的估計精度。當醫(yī)院科研數據量較大時,為方便缺失值的估計,就需用MA TL AB 進行編程,另外,在實際工作中為了了解多個指標間的關系及變化規(guī)律,常常需要同時觀測個體的多個反應指標,此時,雖然數據不是重復觀測數據,但如果缺失數據隨機發(fā)生,即缺失數據發(fā)生的可能性不被該數據值的大小所影響,則同樣可用EM 算法進行估計,因此EM 算法不但適用于重復觀測缺失數據的處理,而且適用于一般情況下多變