EM算法1-原理

EM算法用於含有隱含變量的機率模型參數的極大似然估計。什麼是隱含變量的機率模型呢？舉個例子，假設有3枚硬幣，分別記爲A,B,C，它們正面出現的機率分別爲r,p,q。每次實驗先擲硬幣A，若是出現的是正面就投B，若是出現的反面就投C，出現正面記爲1，出現反面記爲0。獨立10次實驗，觀測結果以下：1101001011。若是隻有這個結果，而不知道過程，問如何估計r,q,p?也就是說，咱們能看到每次的觀測結果，但是這個結果是B產生的仍是C產生的咱們不知道，也就是A的結果咱們不知道，這個就是所謂的隱含變量。若是把觀測變量用Y來表示，隱含變量（A的結果）用Z表示，那麼觀測數據的似然函數爲：

\(P(Y|\theta)=\prod_i{rp^{y_i}(1-p)^{1-y_i}+(1-r)q^{y_i}(1-q)^{1-y_i}}\)

把上面的模型泛化一下能夠歸納爲，有觀測數據{\(x_1,x_2,...x_m\)}，由一個具備觀測變量X和隱含變量Z的模型產生，模型參數爲\(\theta\)，咱們要最大化下面這個似然:

\(l(\theta)=\displaystyle\sum_{i}^{m}logp(x_i;\theta)=\displaystyle\sum_{i}^{m}log\sum_{z_i}p(x_i,z_i;\theta)\)。

直接求解這個優化問題很是困難。EM算法是經過迭代的方式求解，分爲Expectation步和Maximization步。它的主要思想先找到目標函數的下邊界，而後逐步提升這個下邊界，進而獲得一個最優解，可是這個最優解不必定是全局最優的。

下面看一下這個下界是怎麼推導出來的---

\(\displaystyle\sum_{i}^{m}logp(x_i;\theta)\)

\(=\displaystyle\sum_{i}^{m}log\sum_{z_i}p(x_i,z_i;\theta)\)

對於i,假設\(Q_i\)是Z上的某個機率分佈

\(=\displaystyle\sum_{i}^{m}log\sum_{z_i}Q_i(z_i)\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}\)

\(>=\displaystyle\sum_{i}^{m}\sum_{z_i}Q_i(z_i)log\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}\) ---- (eq1)

這一步用到了Jensen不等式，由於log函數是concave的(二階導數小於0)，因此有log(E(X))>=E(log(x))。而\(Q_i\)是機率分佈，因此能夠把\(\sum_{z_i}Q_i(z_i)\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}\)看作是指望，而後能夠套用Jensen不等式[2]，便可獲得上述結果。

如今有了下限，可是裏面的Qi還不知道。怎麼肯定Qi呢？若是咱們已經有了\(\theta\)的一個猜想值，那麼這裏天然讓下限在\(\theta\)處值和似然函數在\(\theta\)處的值越接近越好，讓不等式eq1在\(\theta\)處取得等號。由於log函數是嚴格凹函數，因此只有在E(X)==X(恆等於)的時候等號纔會成立，好比當X是個常數的時候。基於上面的性質，令

\(\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}=c\)

基於此，能夠推出

\(\frac{\sum_zp(x_i,z;\theta)}{\sum_zQ_i(z)}=c\) (這個很容易推出來,a1/b1=c,a2/b2=c,a3/b3=c => (a2+a2+a3)/(b1+b2+b3)=c)

有，

\(Q_i(z_i)=\frac{p(x_i,z_i;\theta)}{\sum_zp(x_i,z;\theta)}\)

\(=\frac{p(x_i,z_i;\theta)}{p(x_i;\theta)}\)

\(=p(z_i|x_i;\theta)\)

所以，Qi是給定xi和\(\theta\)下zi的後驗機率。

這就是E步驟，總結一下，假設已知\(\theta\)，先求出似然函數的下限，而後求出隱含變量的分佈Qi。

在接下來的M步驟，由於E步驟已經獲得了一個Qi，這個步驟求最大化eq1的\(\theta\)值，也就是求下限的最大值點。

而後把M步驟求得的\(\theta\)輸入E步驟，循環往復，直至收斂。

Repeat until convergence{

E-step:for each i,set

\(Q_i(z_i):=p(z_i|x_i;\theta)\)

M-step:set

\(\theta:=argmax_{\theta}\displaystyle\sum_{i}^{m}\sum_{z_i}Q_i(z_i)log\frac{p(x_i,z_i|\theta)}{Q_i(z_i)}\)

}

下面圖片更直觀的描述了EM過程，圖片來自[4],E步挪動下界到\(\theta\)值與目標函數相同，M步求下界函數的最大值點作爲新的\(\theta\)

 

若是定義\(J(Q,\theta)=\displaystyle\sum_{i}^{m}\sum_{z_i}Q_i(z_i)log\frac{p(x_i,z_i;\theta)}{Q_i(z_i)}\)

那麼，EM算法能夠看作函數J的座標軸降低過程，E步最大化Q，M步最大化\(\theta\)。

EM算法是會收斂的，具體的證實參見參考[3]，可是EM算法有可能陷入局部最優的，它對初始值敏感。

 

下面嘗試用EM算法解決文章開始的三硬幣問題。

假設已經通過了j步的迭代，如今已經有了\(\theta^j=(r^j,c^j,q^j)\) (爲了不寫起來混淆，把參數裏面B的正面機率p變成了c)

E步:

這裏要求\(p(z_i|x_i;\theta^j)\)，由於是二分問題，爲了描述簡便，能夠直接求正面的機率，根據貝葉斯機率公式：

(爲了書寫簡單，下面把表示迭代次數的上角標j去掉了，內心記得，r,c,q是已知的)

\(p(z_i=1|x_i;\theta)=\frac{p(x_i|z_i=1;\theta)p(z_i=1;\theta)}{p(x_i|z_i=1;\theta)p(z_i=1;\theta)+p(z_i=0;\theta)p(x_i|z_i=0;\theta)}\)

\(=\frac{rc^{x_i}(1-c)^{(1-x_i)}}{rc^{x_i}(1-c)^{1-x_i}+(1-r)q^{x_i}(1-q)^{(1-x_i)}}\)

把\(p(z_i=1|x_i;\theta)\)記作\(\mu^{(j+1)}\)表示是第j+1次迭代獲得的值，爲了書寫清楚一下(cnblog對公式支持的有些差啊)，仍是把上角標去掉了

M步驟:

如今\(p(z_i|x_i;\theta)\)已經知道，就開始解決下面這個優化問題了

\(J(\theta)=\sum\mu_ilog\frac{p(x_i,z_i=1;\theta)}{\mu_i}+(1-\mu_i)log\frac{p(x_i,z_i=0;\theta)}{1-\mu_i}\)

\(=\sum\mu_ilog\frac{rc^{x_i}(1-c)^{(1-x_i)}}{\mu_i}+(1-\mu_i)log\frac{(1-r)q^{x_i}(1-q)^{(1-x_i)}}{1-\mu_i}\)

令\(\frac{\partial J(\theta)}{r}=0\)

很容易能夠獲得\(r=\frac{1}{m}\sum\mu_i\)

令\(\frac{\partial J(\theta)}{c}=0\)

一樣容易獲得\(c=\frac{\sum\mu_ix_i}{\sum\mu_i}\)

令\(\frac{\partial J(\theta)}{q}=0\)

一樣容易獲得\(c=\frac{\sum(1-\mu_i)x_i}{\sum(1-\mu_i)}\)  參考[1][5]

 

 

參考：

[1]李航《統計學習方法》

[2]Jensen不等式：http://www.cnblogs.com/naniJser/p/5642288.html

[3]Andrew Ng機器學習課程的講義:http://cs229.stanford.edu/notes/cs229-notes8.pdf

[4]介紹EM算法的blog:http://blog.csdn.net/zouxy09/article/details/8537620

[5]http://chenrudan.github.io/blog/2015/12/02/emexample.html