從似然函數到EM算法(附代碼實現)

1. 什麼是EM算法

最大指望算法（Expectation-maximization algorithm，又譯爲指望最大化算法），是在機率模型中尋找參數最大似然估計或者最大後驗估計的算法，其中機率模型依賴於沒法觀測的隱性變量。

最大指望算法通過兩個步驟交替進行計算，

第一步是計算指望（E），利用對隱藏變量的現有估計值，計算其最大似然估計值；
第二步是最大化（M），最大化在E步上求得的最大似然值來計算參數的值。M步上找到的參數估計值被用於下一個E步計算中，這個過程不斷交替進行。

極大似然估計用一句話歸納就是：知道結果，反推條件θ。

1.1 似然函數

在數理統計學中，似然函數是一種關於統計模型中的參數的函數，表示模型參數中的似然性。「似然性」與「或然性」或「機率」意思相近，都是指某種事件發生的可能性。而極大似然就至關於最大可能的意思。

好比你一位同窗和一位獵人一塊兒外出打獵，一隻野兔從前方竄過。只聽一聲槍響，野兔應聲到下，若是要你推測，這一發命中的子彈是誰打的？你就會想，只發一槍便打中，因爲獵人命中的機率通常大於你那位同窗命中的機率，從而推斷出這一槍應該是獵人射中的。

這個例子所做的推斷就體現了最大似然法的基本思想。

多數狀況下咱們是根據已知條件來推算結果，而最大似然估計是已經知道告終果，而後尋求使該結果出現的可能性最大的條件，以此做爲估計值。

1.3 極大似然函數的求解步驟

假定咱們要從10萬我的當中抽取100我的來作身高統計，那麼抽到這100我的的機率就是(機率連乘)：

\[L(\theta)=L(x_1,...,x_n|\theta)=\prod_{i=1}^{n}p(x_i|\theta),\theta\in\ominus\]

如今要求的就是這個 \(\theta\) 值，也就是使得 \(L(\theta)\) 的機率最大化，那麼這時的參數 \(\theta\) 就是所求。

爲了便於分析，咱們能夠定義對數似然函數，將其變成連加的形式：

\[H(\theta)=lnL(\theta)=ln\prod_{i=1}^{n}p(x_i|\theta)=\sum_{i=1}^{n}lnp(x_i|\theta)\]

對於求一個函數的極值，經過咱們在本科所學的微積分知識，最直接的設想是求導，而後讓導數爲0，那麼解這個方程獲得的θ就是了（固然，前提是函數L(θ)連續可微）。但，若是θ是包含多個參數的向量那怎麼處理呢？固然是求L(θ)對全部參數的偏導數，也就是梯度了，從而n個未知的參數，就有n個方程，方程組的解就是似然函數的極值點了，最終獲得這n個參數的值。

求極大似然函數估計值的通常步驟：

1. 寫出似然函數；
2. 對似然函數取對數，並整理；
3. 求導數，令導數爲0，獲得似然方程；
4. 解似然方程，獲得的參數即爲所求；

1.4 EM算法

兩枚硬幣A和B，假定隨機拋擲後正面朝上機率分別爲PA，PB。爲了估計這兩個硬幣朝上的機率，我們輪流拋硬幣A和B，每一輪都連續拋5次，總共5輪：

硬幣	結果	統計
A	正正反正反	3正-2反
B	反反正正反	2正-3反
A	正反反反反	1正-4反
B	正反反正正	3正-2反
A	反正正反反	2正-3反

硬幣A被拋了15次，在第一輪、第三輪、第五輪分別出現了3次正、1次正、2次正，因此很容易估計出PA，相似的，PB也很容易計算出來(真實值)，以下：

PA = （3+1+2）/ 15 = 0.4
PB= （2+3）/10 = 0.5

問題來了，若是咱們不知道拋的硬幣是A仍是B呢（即硬幣種類是隱變量），而後再輪流拋五輪，獲得以下結果：

硬幣	結果	統計
Unknown	正正反正反	3正-2反
Unknown	反反正正反	2正-3反
Unknown	正反反反反	1正-4反
Unknown	正反反正正	3正-2反
Unknown	反正正反反	2正-3反

OK，問題變得有意思了。如今咱們的目標沒變，仍是估計PA和PB，須要怎麼作呢？

顯然，此時咱們多了一個硬幣種類的隱變量，設爲z，能夠把它認爲是一個5維的向量（z1,z2,z3,z4,z5)，表明每次投擲時所使用的硬幣，好比z1，就表明第一輪投擲時使用的硬幣是A仍是B。

• 可是，這個變量z不知道，就沒法去估計PA和PB，因此，咱們必須先估計出z，而後才能進一步估計PA和PB。
• 可要估計z，咱們又得知道PA和PB，這樣咱們才能用極大似然機率法則去估計z，這不是雞生蛋和蛋生雞的問題嗎，如何破？

答案就是先隨機初始化一個PA和PB，用它來估計z，而後基於z，仍是按照最大似然機率法則去估計新的PA和PB，而後依次循環，若是新估計出來的PA和PB和咱們真實值差異很大，直到PA和PB收斂到真實值爲止。

咱們不妨這樣，先隨便給PA和PB賦一個值，好比：
硬幣A正面朝上的機率PA = 0.2
硬幣B正面朝上的機率PB = 0.7

而後，咱們看看第一輪拋擲最多是哪一個硬幣。
若是是硬幣A，得出3正2反的機率爲 0.20.20.20.80.8 = 0.00512
若是是硬幣B，得出3正2反的機率爲0.70.70.70.30.3=0.03087
而後依次求出其餘4輪中的相應機率。作成表格以下：

輪數	如果硬幣A	如果硬幣B
1	0.00512，即0.2 0.2 0.2 0.8 0.8，3正-2反	0.03087，3正-2反
2	0.02048，即0.2 0.2 0.8 0.8 0.8，2正-3反	0.01323，2正-3反
3	0.08192，即0.2 0.8 0.8 0.8 0.8，1正-4反	0.00567，1正-4反
4	0.00512，即0.2 0.2 0.2 0.8 0.8，3正-2反	0.03087，3正-2反
5	0.02048，即0.2 0.2 0.8 0.8 0.8，2正-3反	0.01323，2正-3反

按照最大似然法則：
第1輪中最有可能的是硬幣B
第2輪中最有可能的是硬幣A
第3輪中最有可能的是硬幣A
第4輪中最有可能的是硬幣B
第5輪中最有可能的是硬幣A

咱們就把機率更大，即更多是A的，即第2輪、第3輪、第5輪出現正的次數二、一、2相加，除以A被拋的總次數15（A拋了三輪，每輪5次），做爲z的估計值，B的計算方法相似。而後咱們即可以按照最大似然機率法則來估計新的PA和PB。

PA = （2+1+2）/15 = 0.33
PB =（3+3）/10 = 0.6

就這樣，不斷迭代 不斷接近真實值，這就是EM算法的奇妙之處。

能夠期待，咱們繼續按照上面的思路，用估計出的PA和PB再來估計z，再用z來估計新的PA和PB，反覆迭代下去，就能夠最終獲得PA = 0.4，PB=0.5，此時不管怎樣迭代，PA和PB的值都會保持0.4和0.5不變，因而乎，咱們就找到了PA和PB的最大似然估計。

總結一下計算步驟：

1. 隨機初始化分佈參數θ

2. E步，求Q函數，對於每個i，計算根據上一次迭代的模型參數來計算出隱性變量的後驗機率（其實就是隱性變量的指望），來做爲隱藏變量的現估計值：

   \[Q_i(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)\]

3. M步，求使Q函數得到極大時的參數取值）將似然函數最大化以得到新的參數值

   \[\theta=argmax\sum_{i}\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\]

4. 而後循環重複二、3步直到收斂。

詳細的推導過程請參考文末的參考文獻。

2. 採用 EM 算法求解的模型有哪些？

用EM算法求解的模型通常有GMM或者協同過濾，k-means其實也屬於EM。EM算法必定會收斂，可是可能收斂到局部最優。因爲求和的項數將隨着隱變量的數目指數上升，會給梯度計算帶來麻煩。

3.代碼實現

高斯混合模型 EM 算法

4. 參考文獻

如何通俗理解EM算法

做者：@mantchs

GitHub：https://github.com/NLP-LOVE/ML-NLP

歡迎你們加入討論！共同完善此項目！羣號：【541954936】