EM算法

本文轉載自查看原文 2017-03-17 13:13 2028 23-模式識別/ Jensen不等式/ EM算法/ 混合分布

作者：桂。

時間：2017-03-17 12:18:15

鏈接：http://www.cnblogs.com/xingshansi/p/6557665.html

前言

本文為曲線與分布擬合一文的補充，由於混合分布的推導中用到EM算法，在此梳理一下。全文主要包括：

　　1）EM算法背景介紹；

　　2）EM算法原理推導；

　　3）EM應用

內容多有借鑒他人，最后一並附上鏈接。

一、背景介紹

　　A-硬幣第一拋

本文以最大似然估計（MLE）為例（MAP等同樣有效）。盒子里僅有硬幣A，假設其正面出現的概率為$p$，並將出現正面記作1，反面記作0。進行$N$次獨立重復試驗（取$N=10$）,

得到結果為：$y = ${1，1，0，1，0，0，1，1，1，1}。利用最大似然估計：

$J\left( p \right) = \mathop \prod \limits_{i = 1}^N {p^{{y_i}}}{\left( {1 - p} \right)^{1 - {y_i}}} \Rightarrow \sum\limits_{i = 1}^N {{y_i}\ln p + \left( {1 - {y_i}} \right)} \ln \left( {1 - p} \right)$

將$y_i$的結果值代入，得：

$L\left( p \right) = \ln J\left( p \right) = 7\ln p + 3\ln \left( {1 - p} \right)$

求解：$p = \frac{7}{10}$。

題外話：

MLE：所見即所得，$y_i = 1$個數為7，事實上統計概率$7/10$就是$p$的最佳估計。

MLE與熵（Entropy）同樣存在聯系，重寫代價函數：

$L\left( p \right) = \ln J\left( p \right) = 7\ln p + 3\ln \left( {1 - p} \right)$

微調：

$L\left( p \right) = \ln J\left( p \right) = 7/10\ln p + 3/10\ln \left( {1 - p} \right) = q\ln p + \left(1-q \right)\ln \left( {1 - p} \right)$

若要熵最大，有$q = p$，$p$為待估計的概率，$q$為統計概率0.7。

　　B-硬幣第二拋

回到拋硬幣的問題，接着拋硬幣：盒子里僅有硬幣B，假設其正面出現的概率為$s$，並將出現正面記作1，反面記作0，進行$N$次獨立重復試驗（取$N=10$）,

得到結果為：$y = ${1，0，0，1，0，0，1，0，0，1}。利用最大似然估計，同樣可以估計出：$s = 0.4$.

從這兩次拋硬幣來看，最大似然估計MLE都讓問題得到解決。下面進行第三拋，MLE能否勝任呢？接着往下看。

　　C-硬幣第三拋

可能是因為趕時間，這次拿出硬幣就拋了，卻沒有留意硬幣A、B都在盒子里！將出現正面記作1，反面記作0，進行$N$次獨立重復試驗（取$N=10$）,得到觀測結果為：$y = ${1，1，0，1，0，0，1，0，1，1}，下面問題來了：根據觀測數據$y$，分別估計A、B硬幣正面的出現概率$p、s$?連哪一次是A、哪一次是B都不知道，怎么估計概率呢？MLE仍然不死心：假設A、B兩枚硬幣，拿出A的概率記為$\pi$，拿出B的概率為$1-\pi$，得到准則函數：

$J\left( {\pi,p,s} \right) = \mathop \prod \limits_{i = 1}^N \left[ {\pi {p^{{y_i}}}{{\left( {1 - p} \right)}^{1 - {y_i}}} + \left( {1 - \pi } \right){s^{{y_i}}}{{\left( {1 - s} \right)}^{1 - {y_i}}}} \right]$

這個問題涉及三個未知參數$\theta = (\pi,p,s)$，MLE強行優化准則：

$\hat \theta = \arg \mathop {\max }\limits_\theta \log J\left( {\pi ,p,s} \right)$

再求偏導看看？MLE略顯沮喪，EM走過來拍了拍他的肩膀：困難總會有解決辦法，不是嗎？

二、原理推導

　　A-算法步驟

硬幣第三拋的問題可以簡化為：

輸入：觀測變量數據Y，隱變量數據Z（可以將從盒中選A/B的概率看作是拋硬幣Z決定，例如選硬幣A的概率為$\pi$，選硬幣B的概率為$1-\pi$，則等價為：拋硬幣Z，正面則選擇拋A，反面則選擇拋B，且Z的正面概率為$\pi$，反面概率為$1-\pi$，由於觀測不到Z，硬幣Z就是隱變量數據。）聯合分布P（Y,Z|$\theta$），條件分布P(Z,|Y,$\theta$);

輸出：模型參數$\theta$.

對於輸入模型，求解似然函數：

對於只有觀測變量的情形（記為：完全數據情形），准則函數L借助MLE可解；但隱變量$Z$的存在（記為：缺失數據情形），$L(\theta)$求解沒有閉式解，但如果$P(Z|\theta)$變為一個常數呢？即可將缺失數據情形轉化為完全數據情形，從而借助MLE求解。按下面的步驟求解：

山寨版EM算法——

步驟1：給定一個初始值$P(Z|\theta)$；

步驟2：利用MLE：對$L(\theta)$求偏導，借助梯度下降求解$\theta$；

步驟3：再將求解的$\theta$看作常數，對$L(\theta)$關於$P(Z|\theta)$求解，得出新的$P(Z|\theta)$估計；

步驟4：重復步驟1，直到滿足收斂條件。

山寨版EM算法，利用了循環迭代的算法優化$L(\theta)$，對數$log（.）$形式的求導造成分母有求和/積分項，可見即使轉化為完全數據情形，求解也非常困難。

觀察隱變量特性：通常$Z$一種取值決定了一組$\theta$，如$Z$取正面則觀測數據影響$A$的概率求解，$Z$取反面則觀測數據影響$B$的概率求解，如圖所示：

如果將對數{關於Z求和}的形式，轉換為關於Z求和{對數}的形式，則針對每一種可能的Z求偏導便可以避免上述求解的麻煩，這是EM的關鍵所在。EM算法最大的特色不在於它可以對含有隱變量的問題進行估計，而在於：提供了一種近似$L(\theta)$的准則函數$Q$，實現了求和與對數的轉化，並證明了准則函數$Q$的有效性。

給出求解步驟：

正規版EM算法——

步驟1：選擇參數的初值${\theta ^{\left( 0 \right)}}$，開始迭代；

步驟2：E步（求$Q\left( \theta ,\theta ^\left( i \right) \right)$）：記${\theta ^{\left( i \right)}}$為第i次迭代參數$\theta$的估值，在第i+1次迭代的E步，計算:

其中，${P\left( {Z|Y,{\theta ^{\left( i \right)}}} \right)}$是給定觀測數據Y和當前參數估計$\theta ^\left( i \right)$下隱變量Z的條件概率分布。

事實上，E步主要求解隱變量條件概率密度$P\left( {Z|Y,{\theta ^{\left( i \right)}}} \right)$，這一步實現了缺失數據—>完全數據的轉化，進而構造准則函數$Q$，為MLE求解參數（即M步）作准備；

步驟3：M步（在隱變量條件概率密度給定的前提下，利用MLE實現參數估計）求使$Q\left( {\theta ,{\theta ^{\left( i \right)}}} \right)$最大化的$\theta$，確定第$i+1$次迭代的參數的估計值${{\theta ^{\left( {i + 1} \right)}}}$:

步驟4：重復步驟2、3，直到滿足收斂條件。

EM算法的精華在於Q函數的給出，關於Q的來龍去脈，將在本文后半部分給出。

　　B-遺留問題求解

再回過頭來看看MLE感到沮喪的問題——硬幣第三拋。有了EM算法，我們來一步步分解這個問題。

記：觀測數據為$Y$={$Y_1,Y_2,...Y_N$}，對應隱變量為$Z$={$Z_1,Z_2,...Z_N$}。

E-Step:

1)將缺失數據，轉化為完全數據

假設已經得到第$i$次迭代的參數$\theta ^{\left( i \right)}$，$\theta$ = {$\pi,p,s$}.

E步的核心就是缺失數據到完全數據的轉化，即計算條件概率$P\left( {Z|Y,{\theta ^{\left( i \right)}}} \right)$。

對於任意$Z_j$以及對應的$Y_j$，條件概率$P\left( {Z_j|Y_j,{\theta ^{\left( i \right)}}} \right)$有兩種取值：

$P\left( {{Z_j} = 1|{Y_j},{\theta ^{\left( i \right)}}} \right)$

$P\left( {{Z_j} = 0|{Y_j},{\theta ^{\left( i \right)}}} \right)$

二者概率之和為1，故只計算其中一個即可，假設計算$P\left( {{Z_j} = 1|{Y_j},{\theta ^{\left( i \right)}}} \right)$。利用條件概率公式：

其中：

${P\left( {{Z_j} = 1|{\theta ^{\left( i \right)}}} \right)}$表示硬幣Z第$j$次拋正面的概率，對應概率值為$\pi^{(i)}$；

${P\left( {{Y_j}|{Z_j} = 1,{\theta ^{\left( i \right)}}} \right)}$表示硬幣$Z$第$j$次拋正面時對應的$Y$的概率。Z拋正面對應選擇硬幣$A$，此時概率有兩種可能：

$Y_j = 0$	$Y_j = 1$
$p^{(i)}$	$1-p^{(i)}$

${P\left( {{Y_j}|{Z_j} = 1,{\theta ^{\left( i \right)}}} \right)}$即為表格中的概率分布，簡化成一個表達式：

分母為分子所有組合的求和。

硬幣有兩面，如果是3/4/..更多可能呢？事實上，每一個觀測點可以表示為$P\left( {{Z_j} \in {\Upsilon _k}|{Y_j},{\theta ^{\left( i \right)}}} \right)$，${{Z_j} \in {\Upsilon _k}}$表示第$j$個觀測點來自概率模型$k$，對應到這里就是：${\Upsilon _1}$對應硬幣正面，${\Upsilon _2}$對應硬幣反面，至於正面/反面記作0還是1或其他值，則無影響。

至此，完成了條件概率的求解，缺失數據變為完全數據。按照前文的分析，得到完全數據，對原始准則函數$L$求偏導亦可解。

但EM的精髓在於：將 $log[求和f(或積分)]$的形式轉化為$求和(或積分)log(f)$的形式，簡化求解。為此，我們需要寫出新的准則函數Q。需要知道：L也好，Q也好，都是基於完全數據的求解，所以說條件概率的求解是E-step的核心。

2）構造准則函數Q

給出准則函數Q的固定形式：

展開為具體形式：

條件概率已經求出，又${\log P\left( {{Y_j},{Z_j} = 0|\theta } \right)}$同${\log P\left( {{Y_j},{Z_j} = 1|\theta } \right)}$一樣，在求解條件概率的中間過程已經求出（上文有對應公式，不過注意$log(.)$的參數是$\theta$，而不是$\theta^{(i)}$），再次印證一點：E-Step重要的是條件概率的求解，一旦求解得出，就會有兩件事發生：1）缺失數據變為完全數據；2）Q函數可直接給出；

硬幣第三拋對應的Q函數：

其中$\mu _j^{(i + 1)}$:

M-Step：

1）利用MLE求解參數

有了Q函數，利用MLE分別求偏導：

至此，EM算法完成求解，MLE笑從雙臉生。

　　C-算法推導

1）凸函數

粗略翻譯：

$f$是定義域為實數的函數，如果$x$是實數，滿足$f^{''} \ge 0$；如果$x$為向量，滿足Hessian矩陣$H \ge 0$；這兩種情況下$f$為凸函數，當不取等號時，$f$為嚴格凸函數。

2）Jensen不等式

粗略翻譯：

$f$是凸函數，且$X$為隨機變量，則：

$E\left( {f\left( X \right)} \right) \ge f\left( {E\left( X \right)} \right)$

若$f$是嚴格凸函數，當且僅當$P(X = E(X)) = 1$時上式等號成立（此時，$X$為一常數）。

示意圖：

下面我們證明一下Jensen不等式：

假設$f$為凸函數，且區間為$I$，有${x_1},{x_2},...,{x_n} \in I$，對應概率密度為${p_1},{p_2},...,{p_n}$，且$\sum\limits_{i = 1}^n {{p_i}} = 1$。利用數學歸納法加以證明。

$n=1$沒有討論的意義；當$n=2$，對應凸函數定義，故結論顯然成立；

假設不等式對$n=K$有效，當$n=K+1$時，

至此，Jensen不等式得證。