EM算法(Expectation maximization algorithm)

• • 極大似然估計
  • 極大似然估計定義
    • 求解過程
  • EM算法
    • EM算法和極大似然估計的區別
    • 雞生蛋，蛋生雞問題
    • EM算法思想
    • 三硬幣模型
    • jensen不等式
    • Q函數
    • EM算法的推導
  • K-means中EM思想
  • 參考資料

糾結了好幾天，總算搞清楚了EM算法的大概。因此寫下這篇博客做個筆記，由於這方面懂得不是很多，可能存在理解錯誤的地方，歡迎大家指正，好了閒話不多說。

極大似然估計

極大似然估計定義

在正式介紹EM算法之前，我們需要先來了解一下最大似然估計。這個我們應該都在概率論中學過，其實思想比較簡單，而且我們在生活中也經常用到，舉一個簡單的例子：
$$ 某位同學與一位獵人一起外出打獵，一隻野兔從前方竄過。只聽一聲槍響，野兔應聲到下，如果要你推測，這一發命中的子彈是誰打的？你就會想，只發一槍便打中，由於獵人命中的概率一般大於這位同學命中的概率，看來這一槍是獵人射中的。這便是最大似然的思想，看起來是不是非常簡單。下面來看看極大似然的定義。

定義：極大似然估計是建立在極大似然原理的基礎上的一個統計方法，是概率論在統計學中的應用。極大似然估計提供了一種給定觀察數據來評估模型參數的方法，即：「模型已定，參數未知」。通過若干次試驗，觀察其結果，利用試驗結果得到某個參數值能夠使樣本出現的概率爲最大，則稱爲極大似然估計。

求解過程

現在我們以一個正態分佈爲例，假設有一組樣本

D=(x1,x2,x3....xn) $$D=(x_1,x_2,x_3....x_n)$$

我們知道樣本 x1,x2,...xn $x_1,x_2,...x_n$獨立同分佈於一個正太分佈函數：

f(x)=12π−−√δe−(x−μ)22δ2 $$f(x)=\frac{1}{\sqrt{2\pi }\delta }{e}^{-\frac{(x-\mu {)}^2}{2{\delta }^2}}$$

並且 δ $\delta$是已知的，而 μ $\mu$未知，那麼最大似然需要最什麼呢？求出當均值爲多少時，產生這種採樣數據的概率最大。
我們令似然函數 l(θ)=P(D|θ) $l(\theta )=P(D|\theta )$，這裏我們解釋下什麼是 θ $\theta$，其實 θ $\theta$表示的是是當前概率最大的最大似然函數的模型，什麼意思呢？即取到當前樣本最大的時候對應的函數參數 (μ,δ2) $(\mu ,{\delta }^2)$。在這裏其實等同於 μ $\mu$。

由於樣本之間是獨立同分布所以

l(θ)=∏i=1nf(xi|θ) $$l(\theta )=\underset{i=1}{\overset{n}{\prod }}f(x_i|\theta )$$

下一步要做的便是找到一個和是的 θ^ $\hat{\theta }$使得 l(θ) $l(\theta )$最大，即

θ^=argmaxl(θ) $$\hat{\theta }=argmaxl(\theta )$$

具體做法，對 θ $\theta$求導，然後解出最大值。

∂L(θ)∂θ=0 $$\frac{\mathrm{\partial }L(\theta )}{\mathrm{\partial }\theta }=0$$

然後求解 θ $\theta$。

具體計算過程可以看看極大似然估計詳解，可以更好的理解 θ $\theta$以及整個流程。

所以計算過程總結下來就是

（1）寫出似然函數；

（2）對似然函數取對數，並整理；

（3）求導數，令導數爲0，得到似然方程；

（4）解似然方程，得到的參數即爲所求；

EM算法

既然我們前面講到極大似然估計，那麼EM算法到底和他有什麼關係呢？以生活中送快遞爲例子

EM算法和極大似然估計的區別

極大似然估計面臨的情況
$$一個快遞員給你送貨。若他到你家只有一條路（結果的實現依賴一個概率分佈），但卻不知道這條路今天修不修路（不知道該概率分佈的參數），修路的話今天快遞員就沒法送貨，若結果是快遞員送到貨了，那這條路修路了沒？答案很明顯：沒修路。

EM算法面臨的情況
$$快遞員到你家的路有N條（結果的實現依賴多個概率分佈），但快遞員只會選擇一條路，即，今天他不會選擇第二條路，若他選擇的路修路，那他就不給你送貨了，即使這會而讓你暴跳如雷。問：如果今天快遞員送到貨了，則他選擇的哪條路？那條路修路了嗎？對於這個，因爲你不知道他選擇的哪條路（他把貨送到就走了，根本不給你問他話的時間），所以你唯一能做的就是估計出這N條路被他選擇的N個概率（即：每個概率分佈的權值），然後在根據極大似然估計來得出：這條路沒修路（求出每個概率分佈的參數）。

一句話總結就是極大似然估計是知道概率分佈，不知道參數，現在需要通過求解參數使得當前觀測值的可能性最大，而EM算法是知道觀測值屬於哪一個概率分佈，也不知道參數

雞生蛋，蛋生雞問題

對於EM算法，我們既不知道樣本屬於那個概率分佈，也不知道具體的參數，要求在什麼參數下會使觀測概率最大？這樣就帶來了一個問題，要知道屬於哪個分佈，必須要知道具體的參數。然而參數是未知的。要求極大參數，知道概率分佈是前提。兩個相互依賴，但是又都是未知的，就造成了雞生蛋蛋生雞的問題。

EM算法思想

面對上面的問題，EM算法怎麼做的呢？舉例說明：現在你有一堆糖果，現在有兩個盤子，你需要將糖果均分到兩個盤子中，你又不想一個一個的數，嫌太麻煩了。那可以先隨機把糖果分成兩堆，分別放到兩個盤子當中，然後用手分別拿起盤子掂量一下分量，判斷哪一個重一些，然後將重的那一盤糖果中那一部分出來放入輕的當中，重複這個過程直到兩邊分量感覺起來差不多。

EM算法就是這樣，假設我們想估計知道A和B兩個參數，在開始狀態下二者都是未知的，但如果知道了A的信息就可以得到B的信息，反過來知道了B也就得到了A。可以考慮首先賦予A某種初值，以此得到B的估計值，然後從B的當前值出發，重新估計A的取值，這個過程一直持續到收斂爲止。

上面第一步賦初值對應的EM算法中的E步，求期望(expcetation)，後面的迭代表示的M步，求極大值。

三硬幣模型

假設有三枚硬幣，分別記作A,B,C。這些硬幣正面出現的概率分別爲 π,p,q $\pi ,p,q$。現在進行拋硬幣實現：先拋A，根據A的結果拋出B或C，正面選B，反面選C。然後擲出所選的硬幣。出現正面記爲1，反面爲0。獨立重複n次(n=10)，觀測結果如下:

1,1,0,1,0,0,1,0,1,1 $$1,1,0,1,0,0,1,0,1,1$$

假設只能觀測到B和C的結果，不能觀測A的結果，求三枚硬幣出現正面的概率 π,p,q $\pi ,p,q$
三硬幣模型可以寫作

p(y|θ)=∑zp(y,z|θ)=∑zp(y|z,θ)p(z|θ)(1) $$p(y|\theta )=\underset{z}{\sum }p(y,z|\theta )=\underset{z}{\sum }p(y|z,\theta )p(z|\theta )(1)$$

p(y|θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y(2) $$p(y|\theta )=\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}(2)$$

其中y表示觀測變量，即B、C的結果：1或0。隨機變量z是隱含變量，即A最後的結果，我們是無法觀測。 θ=(π,p,q) $\theta =(\pi ,p,q)$是模型參數。對公示(2)的解釋， p(y|θ) $p(y|\theta )$表示在 θ $\theta$下 y $y$的概率，假設 y=1 $y=1$那麼公式則爲 p(y|θ)=πpy+(1−π)qy $p(y|\theta )=\pi p^y+(1-\pi )q^y$,當 y=0 $y=0$的時候， p(y|θ)=π(1−p)y+(1−π)(1−q)y $p(y|\theta )=\pi (1-p{)}^y+(1-\pi )(1-q{)}^y$，由於 y=1or0 $y=1or0$，將兩種結果同一即得到(2)

將上述觀測數據表示爲 Y=(y1,y2...yn) $Y=(y_1,y_2...y_n)$，爲觀測數據表示爲 Z=(Z1,Z2,....Zn) $Z=(Z_1,Z_2,....Z_n)$。則觀測數據的似然函數爲：

P(Y|θ)=∑zp(Y,Z|θ)p(Z|θ) $$P(Y|\theta )=\underset{z}{\sum }p(Y,Z|\theta )p(Z|\theta )$$

即：

P(Y|θ)=∏j=1n[πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj] $$P(Y|\theta )=\underset{j=1}{\overset{n}{\prod }}[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]$$

模型參數 θ=(π,p,q) $\theta =(\pi ,p,q)$的極大似然估計爲:

θ^=argmaxlogP(Y|θ) $$\hat{\theta }=argmaxlogP(Y|\theta )$$

這個問題就沒有解析解，只有通過迭代方式求解。EM算法就是可以用於求解這個問題的迭代算法。EM算法首先選取參數的初始值，記作 θ(0)=(π(0),p(0),q(0)) ${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$，然後通過下面的步驟迭代計算參數的估計值，直到收斂爲止。 θ(i)=(π(i),p(i),q(i)) ${\theta }^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$標識的是第i次迭代後的模型參數。第 i+1 $i+1$次的迭代我們可以這樣表示：

E步：計算模型參數 π(i),p(i),q(i) ${\pi }^{(i)},p^{(i)},q^{(i)}$下觀測數據 yj $y_j$來自硬幣B的概率：

μ(i+1)j=π(i)(p(i))yj(1−p(i))1−yjπ(i)p(i))yj(1−p(i))1−yj+(1−π(i))(q(i))yj(1−q(i))(1−yj)(3) $${\mu }_j^{(i+1)}=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}}{{\pi }^{(i)}{p}^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_j}(1-q^{(i)}{)}^{(1-y_j)}}(3)$$

這一步其實就是前面公式(2)的第i次迭代，計算出來自B的概率之後接下來需要對參數重新估值。

M步：計算參數模型的新估值。

π(i+1)=1n∑j=1nμ(i+1)j(4) $${\pi }^{(i+1)}=\frac{1}{n}\underset{j=1}{\overset{n}{\sum }}{\mu }_j^{(i+1)}(4)$$

p(i+1)=1n∑nj=1μ(i+1)yj∑nj=1μ(i+1)(5) $$p^{(i+1)}=\frac{1}{n}\frac{\underset{j=1}{\overset{n}{\sum }}{\mu }^{(i+1)}{y}_j}{\underset{j=1}{\overset{n}{\sum }}{\mu }^{(i+1)}}(5)$$

q(i+1)=∑nj=1(1−μ(i+1))yj∑nj=1(1−μ(i+1))(6) $$q^{(i+1)}=\frac{\underset{j=1}{\overset{n}{\sum }}(1-{\mu }^{(i+1)})y_j}{\underset{j=1}{\overset{n}{\sum }}(1-{\mu }^{(i+1)})}(6)$$

這四個公式看起來有點複雜，其實理解起來沒有那麼難。公式(3)和前面公式(2)是一樣的，這裏就不在贅述。現在看看公式(4), π $\pi$其實就是觀測是 yj $y_j$來自B的概率，所以只需要將公式(3)求均值就行了。公式(5)表示的觀測值 yj $y_j$來自B並且爲正面的概率，即求在觀測值來自B的條件下，觀測爲正面的條件概率。關於這個具體計算，可以參考李航統計學習方法裏面。需要注意的一點是：EM算法的參數估計值與選取的初始值有關。

jensen不等式

設f是定義域爲實數的函數，如果對於所有的實數x。如果對於所有的實數x，f(x)的二次導數大於等於0，那麼f是凸函數。當x是向量時，如果其hessian矩陣H是半正定的，那麼f是凸函數。如果只大於0，不等於0，那麼稱f是嚴格凸函數。

Jensen不等式表述如下：

如果f是凸函數，X是隨機變量，那麼： E[f(X)]≥f(E[X]) $E[f(X)]\ge f(E[X])$

特別地，如果f是嚴格凸函數，當且僅當X是常量時，上式取等號。

Q函數

定義： 完全數據的對數似然函數 logP(Y,Z|θ) $logP(Y,Z|\theta )$關於在給定的觀測數據 Y $Y$和當前參數 θ(i) ${\theta }^{(i)}$下對爲觀測數據Z的條件概率分佈 P(Z|Y,θ(i)) $P(Z|Y,{\theta }^{(i)})$的期望，即

Q(θ,θ(i))=EZ[logP(Y,Z|θ)|Y,θ(i)] $$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]$$

=∑ZP(Z|Y,θ(i))logP(Y,Z|θ) $$=\underset{Z}{\sum }P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$

什麼意思呢？其實簡化之後及爲 Q(θ,θ(i)) $$=\underset{Z}{\sum }P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$

什麼意思呢？其實簡化之後及爲 Q(θ,θ(i))=∑ZP(Z)logP(Y,Z) $Q(\theta ,{\theta }^{(i)})=\underset{Z}{\sum }P(Z)logP(Y,Z)$

EM算法的推導

現在我們來看看EM算法的一般推導，一般的，用 Y $Y$表示觀測隨機變量的數據，Z表示隱隨機變量的數據， YY $Y$和 Z $Z$連在一起被稱爲完全數據，觀測數據 Y $Y$又稱爲不完全觀測數。給定觀測數據Y，其概率分佈爲 P(Y|θ) $P(Y|\theta )$， θ $\theta$爲模型參數，那麼不完全數據 Y $Y$的對數似然函數爲 Y