混合Bernoulli分布參數估計的EM算法研究

張寶龍，魏立力

（寧夏大學 數學計算機學院，寧夏 銀川 750021）

成敗型隨機試驗在統計學上稱為伯努利試驗（Bernoulli trial）.很多實際問題都可以歸結為伯努利試驗.比如在醫學領域考察對病人治療結果的有效與無效、某種化驗結果的陽性與陰性、接觸某傳染源的感染與未感染等；在系統可靠性理論中元件工作正常與失效；決定人類的某一特別屬性（比如是否為左撇子）的一對基因的顯性表現與隱性表現；某陪審團的陪審員對被告人的投票結果為有罪和無罪等等.伯努利試驗必須滿足兩個基本條件：每次試驗的結果獨立且只有“成功”與“失敗”，每次試驗中“成功”的概率保持不變.

伯努利試驗的一種推廣是假設每次試驗相互獨立，但其成功概率允許不盡相同.這樣的情形可以用一個混合Bernoulli分布來描述：

其中所有πj≥0,且…,pl)為未知參數.如果希望考察某種藥物對患者的療效（有效或無效），則該模型非常適用，因為我們很難保證同種藥物對不同患者的療效完全相同.也就是說，我們預期對于眾多患者的療效可以分成l個不同的類別.

現在設y=(y1,y2,…,yn)來自于混合Bernoulli分布(1.1)的樣本，我們的目的是求未知參數θ=(π1,π2,…,πl-1,p1,p2,…,pl)的極大似然估計.為此先考查其對數似然函數：

不難看出，直接求(1.2)式的最大值點是很困難的，我們下面將推導該問題的EM算法.

EM算法是一種迭代計算，其每次迭代由兩步組成：E步（求條件期望）和M步（極大化），這正是該算法名稱的由來.該算法最初由Dempster，Laird和Rubin提出[1]，主要用來求后驗分布的眾數（極大似然估計），廣泛應用于刪失數據，截尾數據，成群數據等.其基本思想是在給出缺失數據初值的條件下，估計出模型參數的值；然后再根據參數值估計出缺失數據的值.根據估計出的缺失數據的值再對參數值進行更新，如此反復迭代，直至收斂，迭代結束.

EM算法提出之后，很快引起國內外眾多學者的關注，文獻[2]很好地總結了EM算法及其推廣算法的很多成果.文獻[3]詳細介紹了有限混合模型及其應用.文獻[4]介紹了有限混合模型及其應用的研究進展.本文基于EM算法研究了有限混合Bernoulli分布模型的參數估計，并利用R軟件進行了數值模擬.

1 EM 算法簡介

一般而言，形式上[1]我們有兩個樣本空間X,Y，以及X到Y的一個多對一映射xay(x).其中X中x=(x1,x2,…,xn)不能直接觀測到，只能通過y間接的觀測到，x被稱為“完全數據”.Y里的y=(y1,y2,…,yn)是能夠觀測到的數據，即“不完全數據”.

設參數θ∈Θ,x的密度函數為fc(x|θ),則y的密度函數為

其中X(y)={x:y(x)=y}

我們的目的是求參數θ 的極大似然估計，對坌θ∈Θ,使得g(y|θ)極大化.令x=(y,z)表示y的完全數據，其中z=(z1,z2,…,zn)表示不可觀測數據或缺失數據，即將yi,i=1,…,n用缺值擴張為xi=(yi,zi).由于在具體問題中，極大化不完全數據的密度函數g(y|θ)往往要比極大化完全數據的密度函數fc(x|θ)困難很多.EM算法就是試圖通過極大化后者（或者對數似然函數）而達到極大化g(y|θ)的目的.但是由于x不能被觀測到，從而logfc(x|θ)未知，我們的策略是用logfc(x|θ)在給定y和θ(k)（經過k次迭代后的值）下的條件期望來代替.

具體地說，設θ(0)是θ 的初值.則在第一次迭代中，E步需要計算

M步則需要關于θ 最大化Q(θ;θ(0)).也就是求θ(1)使得

重復執行E步和M步，但是這次用θ(1)的當前值來代替θ(0)，然后極大化得到θ(2)作為下一次迭代的初值.則在k+1次迭代時，E步和M步可以被定義如下:

簡單和穩定是EM算法的最大優點，EM算法所得到的估計序列具有良好的收斂性，且收斂到g(y|θ)的最大值.具體而言，記估計序列為θ(k),k=1,2,…,L(θ|y)=logg(y|θ),則在關于L的很一般的條件下，θ(k)的收斂值θ*是L的穩定點[1].

2 有限混合Bernoulli 分布模型參數估計的EM 算法

設離散型隨機變量X服從有限混合Bernoulli分布，其概率函數為(1.1)式，設)是來自于該分布的簡單隨機樣本，其似然函數為(1.2)式，我們的目的是求未知參數θ=(π1,π2,…,πl-1，p1,p2,…,pl)的極大似然估計.下面我們討論相應的EM算法.

引入潛在變量z=(z1,z2,…,zn)，其中z=(zi1,zi2,…,zil)，且z1,z2,…,zn相互獨立，zij是取值為0或1的指示變量，zij=1表示yi來自于第j個分支密度，且

于是有yi|zij=1:B(1,pj),j=1,2,…,l.這樣可以得到完全樣本為xi=(zi,yi),i=1,2,…,n,其似然函數為

對上述似然函數取對數并去掉與π1,π2,K,πl-1,p1,p2,K,pl無關的量，得

在E步中，令

易知

在M步中，解

得到

基于以上導出的迭代公式，用R軟件進行參數估計的算法如下：

1.產生n個來自有限混合Bernoulli分布的隨機變量；

3.令k=1，將2中的初值代入迭代公式(1.6)中，得到θ(k)；

4.k→k+1，將θ(k)代入迭代公式(1.6)中，得到θ(k+1)；

5.重復4，直至||θ(k+1)-θ(k)||充分小時停止.

以二階混合二項分布為例：模擬1中隨機變量Y的算法為

1.產生服從參數為(1,p1)的二項隨機變量Y1；

2.產生服從參數為(1,p2)的二項隨機變量Y2；

3.產生隨機數U；

4.若U<π，令Y=Y1；否則令Y=Y2；

3 數值模擬

下面用R統計軟件進行隨機模擬.

第一種情況：建立混合模型0.4B(1,0.3)+0.6B(1,0.5)，分別產生500個和1000個來自該混合模型的隨機數，選取三組不同的初值為

表1 混合模型0.4B(1,0.3)+0.6B(1,0.5)參數估計結果

表2 混合模型0.4B(1,0.3)+0.6B(1,0.5)參數估計結果

第二種情況：建立混合模型0.4B(1,0.3)+0.6B(1,0.5)，分布產生100個、500個、1000個來自該混合模型的隨機數，選取初值為π(0)=0.4,p1(0)=0.1,p2(0)=0.2.進行數值模擬，EM算法參數估計結果見表3.

從表1和表2可以明顯看出,隨著初值逐漸接近真值時，估計值亦趨于真值.當估計值變化不大時，說明估計值收斂到穩定點.由表3可以看出，隨著樣本容量的增加，參數的估計值逐漸接近于真值.同樣，當估計值變化不大時，說明估計值收斂到穩定點.

〔1〕Dempster A P，Laird N.Maximum Likelihood from Incomplete Data via EM Algorithm[J].J.Royal Statistical Society，Series B，1977，39:1-38.

〔2〕Gelffrey J.McLachlan.The EM Algorithm and Extensions(Second Edition)[M].New York: Wiley & Sons，Inc，2008.

〔3〕McLachlan G，Peel D.Finite Mixture Models[M].New York:Wiley&Sons，Inc，2000.

〔4〕孫蘭.有限混合模型及其應用的研究進展[D].長春:東北師范大學，2006.

〔5〕魏立力，馬江洪，顏榮芳.概率統計引論[M].北京:科學出版社，2012.