EM算法詳解

• • 1 摘要
  • 2 EM算法簡介
  • 3 預備知識
  • • 3.1 極大似然估計
    • 3.2 Jenson不等式
  • 4 EM詳解
  • • 4.1 問題描述
    • 4.2 EM算法推導流程
    • 4.3 EM算法流程
  • 5. EM算法的收斂性思考
  • • 5.1 對EM算法的初始化研究
    • 5.2 EM算法能保證收斂嗎？
    • 5.3 EM算法能保證收斂到全局最大值嗎？
  • 6. EM算法的一些思考
  • 7. EM算法在非監督學習中的應用

1 摘要

EM算法是最常見的隐變量估計方法，常被用來學習高斯混合模型（Gaussian mixture model，簡稱GMM）的參數；隐式馬爾科夫算法（HMM）、LDA主題模型的變分推斷等等。本文就對EM算法的原理做一個詳細的總結。

2 EM算法簡介

EM算法是一種疊代優化政策，由于它的計算方法中每一次疊代都分兩步，其中一個為期望步（E步），另一個為極大步（M步），是以算法被稱為EM算法（Expectation-Maximization Algorithm）。EM算法受到缺失思想影響，最初是為了解決資料缺失情況下的參數估計問題。通過疊代，不斷求解下界的極大化，來逐漸求解對數似然函數極大化。其基本思想是：

1. 首先根據己經給出的觀測資料，估計出模型參數的值；（根據 x x x 估計 θ θ θ）
2. 然後再依據上一步估計出的參數值估計缺失資料的值;（根據 θ θ θ 估計标簽或者類簇指派 y y y ）
3. 再根據估計出的缺失資料加上之前己經觀測到的資料重新再對參數值進行估計；（根據 y y y + x x x 估計 θ θ θ ）
4. 然後反複疊代1-3，直至最後收斂，疊代結束。

優點： 簡單性和普适性，可看作是一種非梯度優化方法（解決梯度下降等優化方法的缺陷：求和的項數将随着隐變量的數目以指數級上升，會給梯度計算帶來麻煩）

缺點： 對初始值敏感，不同的初值可能得到不同的參數估計值；不能保證找到全局最優值。

3 預備知識

3.1 極大似然估計

（1）定義

目的：利用已知樣本結果，反推最有可能（最大機率）導緻這樣結果的參數值θ。

隻是一種機率論在統計學的應用，它是參數估計的方法之一。已知某個随機樣本 x x x 滿足某種機率分布，但是其中具體的參數不清楚；參數估計就是通過若幹次試驗，觀察其結果，利用結果推出參數的大概值。

最大似然估計是建立在這樣的思想上：已知某個參數能使這個樣本出現的機率最大，我們當然不會再去選擇其他小機率的樣本，是以幹脆就把這個參數作為估計的真實值。

（2）問題描述

假如我們需要調查學校的男生和女生的身高分布 ，我們抽取100個男生和100個女生，将他們按照性别劃分為兩組。然後，統計抽樣得到100個男生的身高資料和100個女生的身高資料。如果我們知道他們的身高服從正态分布，但是這個分布的均值 μ \mu μ 和方差 δ 2 \delta^2 δ2r455 是不知道，這兩個參數就是我們需要估計的。

問題：我們知道樣本所服從的機率分布模型和一些樣本，我們需要求解該模型的參數。如圖1所示。

對數似然方程：

l ( θ ) = l n L ( θ ) = l n ∏ i = 1 n p ( x i ; θ ) = ∑ i = 1 n l n p ( x i ; θ ) l(\theta)=lnL(\theta)=ln\prod_{i=1}^{n}p(x_{i};\theta)=\sum_{i=1}^{n}{lnp(x_{i};\theta)} l(θ)=lnL(θ)=lni=1∏n​p(xi​;θ)=i=1∑n​lnp(xi​;θ)

對于n個樣本觀察資料 x = ( x 1 , x 2 , . . . , x n ) x=(x_{1},x_{2},...,x_{n}) x=(x1​,x2​,...,xn​) ，找出樣本的模型參數θ, 極大化模型分布的對數似然函數如下：

θ ^ = a r g m a x ∑ i = 1 n l o g p ( x i ; θ ) \hat{\theta}=argmax\sum_{i=1}^{n}{logp(x_{i};\theta)} θ^=argmaxi=1∑n​logp(xi​;θ)

我們已知的條件有兩個：樣本服從的分布模型、随機抽取的樣本。我們需要求解模型的參數。根據已知條件，通過極大似然估計，求出未知參數。總的來說：極大似然估計就是用來估計模型參數的統計學方法。

3.2 Jenson不等式

（1）定義

設f是定義域為實數的函數，如果對所有的實數x，f(x)的二階導數都大于0，那麼f是凸函數，X是随機變量，那麼： E [ f ( X ) ] ≥ f ( E [ X ] ) E[f(X)] ≥ f(E[X]) E[f(X)]≥f(E[X]) 。當且僅當自變量 X X X 是常量時，該式取等号。其中，E(X)表示X的數學期望。

圖2中，實線f表示凸函數，X是随機變量，有0.5的機率是a，有0.5的機率是b。X的期望值就是a和b的中值，從圖中可以看到 E [ f ( X ) ] ≥ f ( E [ X ] ) E[f(X)] ≥ f(E[X]) E[f(X)]≥f(E[X]) 成立。

由于對數函數是凹函數，如果f(x)是凹函數, 有:

f ( E ( x ) ) ≥ E ( f ( x ) ) f(E(x))≥E(f(x)) f(E(x))≥E(f(x))

即

l o g ∑ j λ j y j ≥ ∑ j λ j l o g y j , log∑_jλ_jy_j≥∑_jλ_jlogy_j, logj∑​λj​yj​≥j∑​λj​logyj​,

需 要 滿 足 條 件 ： λ j ≥ 0 , ∑ j λ j = 1 需要滿足條件：λ_j≥0, ∑_jλ_j=1 需要滿足條件：λj​≥0,j∑​λj​=1

4 EM詳解

4.1 問題描述

我們目前有100個男生和100個女生的身高，但是我們不知道這200個資料中哪個是男生的身高，哪個是女生的身高，即抽取得到的每個樣本都不知道是從哪個分布中抽取的。這個時候，對于每個樣本，就有兩個未知量需要估計：

（1）這個身高資料是來自于男生資料集合還是來自于女生？（即先驗機率為 π i π_i πi​）

（2）男生、女生身高資料集的正态分布的參數分别是多少？（即期望 μ μ μ 和方差 δ δ δ）

4.2 EM算法推導流程

對于n個樣本觀察資料 x = ( x 1 , x 2 , . . . x n ) x=(x_{1},x_{2},...x_{n}) x=(x1​,x2​,...xn​) ，找出樣本的模型參數θ, 極大化模型分布的對數似然函數如下：

θ ^ = a r g m a x ∑ i = 1 n l o g p ( x i ; θ ) \hat{\theta}=argmax\sum_{i=1}^{n}{logp(x_{i};\theta)} θ^=argmaxi=1∑n​logp(xi​;θ)

如果我們得到的觀察資料有未觀察到的隐含資料 z = ( z 1 , z 2 , . . . , z n ) z=(z_{1},z_{2},...,z_{n}) z=(z1​,z2​,...,zn​) ，即上文中每個樣本屬于哪個分布是未知的，此時我們極大化模型分布的對數似然函數如下：

θ ^ = a r g m a x ∑ i = 1 n l o g p ( x i ; θ ) = a r g m a x ∑ i = 1 n l o g ∑ z i p ( x i , z i ; θ ) \hat{\theta}=argmax\sum_{i=1}^{n}{logp(x_{i};\theta)}=argmax\sum_{i=1}^{n}{log\sum_{z_{i}}^{}{p(x_{i},z_{i};\theta)}} θ^=argmaxi=1∑n​logp(xi​;θ)=argmaxi=1∑n​logzi​∑​p(xi​,zi​;θ)

上面這個式子是根據 x i x_{i} xi​ 的邊緣機率計算得來，沒有辦法直接求出θ。是以需要一些特殊的技巧，使用Jensen不等式對這個式子進行縮放如下：

∑ i = 1 n l o g ∑ z i p ( x i , z i ; θ ) = ∑ i = 1 n l o g ∑ z i Q i ( z i ) p ( x i , z i ; θ ) Q i ( z i ) ( 1 ) ≥ ∑ i = 1 n ∑ z i Q i ( z i ) l o g p ( x i , z i ; θ ) Q i ( z i ) ( 2 ) \sum_{i=1}^{n}{log\sum_{z_{i}}^{}{p(x_{i},z_{i};\theta)}}=\sum_{i=1}^{n}{log\sum_{z_{i}}^{}{Q_{i}(z_{i})\frac{p(x_{i},z_{i};\theta)}{Q_{i}(z_{i})}}} (1) \\ \geq \sum_{i=1}^{n}{\sum_{z_{i}}^{}{Q_{i}(z_{i})log\frac{p(x_{i},z_{i};\theta)}{Q_{i}(z_{i})}}} (2) i=1∑n​logzi​∑​p(xi​,zi​;θ)=i=1∑n​logzi​∑​Qi​(zi​)Qi​(zi​)p(xi​,zi​;θ)​(1)≥i=1∑n​zi​∑​Qi​(zi​)logQi​(zi​)p(xi​,zi​;θ)​(2)

• (1)式是引入了一個未知的新的分布 Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​) ，分子分母同時乘以它得到的。
• (2)式是由(1)式根據Jensen不等式得到的：
  • 由于 ∑ z i Q i ( z i ) l o g [ p ( x i , z i ; θ ) Q i ( z i ) ] \sum_{z_{i}}^{}{Q_{i} (z_{i})log\left[ \frac{p(x_{i},z_{i};\theta)}{Q_{i}(z_{i})} \right]} ∑zi​​Qi​(zi​)log[Qi​(zi​)p(xi​,zi​;θ)​] 為 p ( x i , z i ; θ ) Q i ( z i ) \frac{p(x_{i},z_{i};\theta)}{Q_{i}(z_{i})} Qi​(zi​)p(xi​,zi​;θ)​ 的期望；
  • 且log(x)為凹函數, f ( E ( x ) ) ≥ E ( f ( x ) ) f(E(x))≥E(f(x)) f(E(x))≥E(f(x))，根據Jensen不等式可由(1)式得到(2)式。
• 上述過程可以看作是對 l o g l ( θ ) logl(\theta) logl(θ) 求了下界（ l ( θ ) = ∑ i = 1 n l o g p ( x i ; θ ) l(\theta)=\sum_{i=1}^{n}{logp(x_{i};\theta)} l(θ)=∑i=1n​logp(xi​;θ) ）

問題1：對于 Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​) 我們如何選擇呢？

假設 θ θ θ 已經給定，那麼 l o g l ( θ ) logl(\theta) logl(θ) 的值取決于 Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​)和 p ( x i , z i ) p(x_{i},z_{i}) p(xi​,zi​) 。我們可以通過調整這兩個機率使(2)式 下界不斷上升，來逼近 l o g l ( θ ) logl(\theta) logl(θ) 的真實值。由下文可知， Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​) 即後驗機率

問題2：那麼如何算是調整好呢？

當不等式變成等式時，說明我們調整後的機率能夠等價于 l o g l ( θ ) logl(\theta) logl(θ) 了。按照這個思路，我們要找到等式成立的條件。

此時如果要滿足Jensen不等式的等号，則有：

p ( x i , z i ; θ ) Q i ( z i ) = c ， c 為 常 數 \frac{p(x_{i},z_{i};\theta)}{Q_{i}(z_{i})}=c，c為常數 Qi​(zi​)p(xi​,zi​;θ)​=c，c為常數

由于 Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​)是一個分布，是以滿足：

∑ z Q i ( z i ) = 1 \sum_{z}^{}{Q_{i}(z_{i})}=1 z∑​Qi​(zi​)=1

由此可得（分子分母同時加上 ∑ z \sum_{z} ∑z​）

∑ z p ( x i , z i ; θ ) = c \sum_{z}^{}{p(x_{i},z_{i};\theta)}=c z∑​p(xi​,zi​;θ)=c

從上面兩式，我們可以得到：

　　　　 Q i ( z i ) = p ( x i , z i ; θ ) ∑ z p ( x i , z i ; θ ) = p ( x i , z i ; θ ) p ( x i ; θ ) = p ( z i ∣ x i ; θ ) Q_{i}(z_{i})=\frac{p(x_{i},z_{i};\theta)}{\sum_{z}^{}{p(x_{i},z_{i};\theta)}}=\frac{p(x_{i},z_{i};\theta)}{p(x_{i};\theta)}=p(z_{i}|x_{i};\theta) Qi​(zi​)=∑z​p(xi​,zi​;θ)p(xi​,zi​;θ)​=p(xi​;θ)p(xi​,zi​;θ)​=p(zi​∣xi​;θ)

• 至此，我們推出了在固定其他參數 θ θ θ 後， Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​) 的計算公式就是後驗機率，解決了 Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​) 如何選擇的問題。

• 則第(2)式是我們的包含隐藏資料的對數似然的一個下界。如果我們能極大化這個下界，等價于極大化對數似然。即我們需要最大化下式：

  a r g m a x ∑ i = 1 n ∑ z i Q i ( z i ) l o g p ( x i , z i ; θ ) Q i ( z i ) argmax\sum_{i=1}^{n}{\sum_{z_{i}}^{}{Q_{i}(z_{i})log\frac{p(x_{i},z_{i};\theta)}{Q_{i}(z_{i})}}} argmaxi=1∑n​zi​∑​Qi​(zi​)logQi​(zi​)p(xi​,zi​;θ)​

上式也就是EM算法的M步，去掉分母中的常數項，則我們需要極大化的對數似然下界為：

a r g m a x ∑ i = 1 n ∑ z i Q i ( z i ) l o g p ( x i , z i ; θ ) argmax\sum_{i=1}^{n}{\sum_{z_{i}}^{}{Q_{i}(z_{i})logp(x_{i},z_{i};\theta)}} argmaxi=1∑n​zi​∑​Qi​(zi​)logp(xi​,zi​;θ)

那E步呢？

注意到上式中 Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​) 是一個分布，是以 Q i ( z i ) l o g p ( x i , z i ; θ ) {Q_{i}(z_{i})logp(x_{i},z_{i};\theta)} Qi​(zi​)logp(xi​,zi​;θ) 可以了解為 l o g p ( x i , z i ; θ ) logp(x_{i},z_{i};\theta) logp(xi​,zi​;θ) 基于條件機率分布 Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​) 的期望。這一步就是E步，該步建立了 l ( θ ) l(\theta) l(θ) 的下界。

4.3 EM算法流程

現在我們總結一下EM算法流程。

輸入：觀察到的資料 x = ( x 1 , x 2 , . . . x n ) x=(x_{1},x_{2},...x_{n}) x=(x1​,x2​,...xn​)，聯合分布 p ( x , z ; θ ) p(x,z;\theta) p(x,z;θ)，條件分布 p ( z ∣ x , θ ) p(z|x,\theta) p(z∣x,θ) ，最大疊代次數 J J J。

算法步驟：

（1）随機初始化模型參數 θ θ θ 的初值 θ 0 \theta_{0} θ0​ 。

（2） j = 1 , 2 , . . . , J j=1,2,...,J j=1,2,...,J 開始EM算法疊代：

• E步：計算聯合分布的條件機率期望：

  注意： θ θ θ 與 θ j θ_j θj​ 的位置， θ θ θ 在聯合分布中， θ j θ_j θj​ 在條件分布中。

  Q i ( z i ) = p ( z i ∣ x i , θ j ) Q_{i}(z_{i})=p(z_{i}|x_{i},\theta_{j}) Qi​(zi​)=p(zi​∣xi​,θj​)

l ( θ , θ j ) = ∑ i = 1 n ∑ z i Q i ( z i ) l o g p ( x i , z i ; θ ) l(\theta,\theta_{j})=\sum_{i=1}^{n}{\sum_{z_{i}}^{}{Q_{i}(z_{i})logp(x_{i},z_{i};\theta)}} l(θ,θj​)=i=1∑n​zi​∑​Qi​(zi​)logp(xi​,zi​;θ)

• M步：極大化 l ( θ , θ j ) l(\theta,\theta_{j}) l(θ,θj​) ,得到 θ j + 1 \theta_{j+1} θj+1​ :

  此處使用拉格朗日函數的對偶法求解得到：

  θ j + 1 = a r g m a x l ( θ , θ j ) \theta_{j+1}=argmaxl(\theta,\theta_{j}) θj+1​=argmaxl(θ,θj​)

• 如果 θ j + 1 \theta_{j+1} θj+1​ 已經收斂（兩輪 θ θ θ 或者 Q i ( z i ) Q_{i}(z_{i}) Qi​(zi​) 的絕對內插補點小于門限值），則算法結束。否則繼續進行E步和M步進行疊代。

（3）輸出：模型參數 θ θ θ。

5. EM算法的收斂性思考

5.1 對EM算法的初始化研究

上面介紹的傳統EM算法對初始值敏感，聚類結果随不同的初始值而波動較大。總的來說，EM算法收斂的優劣很大程度上取決于其初始參數。

5.2 EM算法能保證收斂嗎？

要證明EM算法收斂，則我們需要證明我們的對數似然函數的值在疊代的過程中一直在增大。即：

∑ i = 1 n l o g p ( x i , θ j + 1 ) > = ∑ i = 1 n l o g p ( x i , θ j ) \sum_{i=1}^{n}logp(x_{i},\theta_{j+1}) >= \sum_{i=1}^{n}logp(x_{i},\theta_{j}) i=1∑n​logp(xi​,θj+1​)>=i=1∑n​logp(xi​,θj​)

參考：EM算法原理總結 - 劉建平Pinard - 部落格園 ；李航《統計機器學習》

5.3 EM算法能保證收斂到全局最大值嗎？

從上面的推導可以看出，EM算法可以保證收斂到一個穩定點，但是卻不能保證收斂到全局的極大值點，是以它是局部最優的算法，當然，如果我們的優化目标 L ( θ , θ j ) L(θ,θj) L(θ,θj) 是凸的，則EM算法可以保證收斂到全局最大值，這點和梯度下降法這樣的疊代算法相同。

　是以在應用中，初值的選擇非常重要。常用的辦法是選取幾個不同的初值進行疊代，然後對得到的各個估計值加以比較，從中選擇最好的。

6. EM算法的一些思考

如果我們從算法思想的角度來思考EM算法，我們可以發現：已知的是觀察資料，未知的是隐含資料和模型參數，

• 在E步，固定模型參數的值，優化隐含資料的分布（類簇的指派 / 所屬）；
• 在M步，固定隐含資料分布，優化模型參數的值（類簇分布的具體參數）。

一個最直覺了解EM算法思路的是K-Means算法。在K-Means聚類時，每個聚類簇的質心是隐含資料。k-means算法是高斯混合聚類在混合成分方差相等，且每個樣本僅指派一個混合成分時候的特例。k-means算法與EM算法的關系是這樣的：

• k-means是兩個步驟交替進行:确定中心點，對每個樣本選擇最近中心點–> E步和M步。
• 假設K個初始化質心
• E 步中： 将每個點選擇最近的類優化目标函數，分給中心距它最近的類(硬配置設定)，可以看成是EM算法中E步（軟配置設定，計算各個類簇的質心對觀測資料的響應度）的近似。
• M步中： 更新每個類的中心點（即更新新一輪疊代的模型參數），可以認為是在「各類分布均為機關方差的高斯分布」的假設下，最大化似然值；

當然，K-Means算法是比較簡單的，高斯混合模型（GMM）也是EM算法的一個應用。

7. EM算法在非監督學習中的應用

EM算法可用于生成模型的非監督學習。生成模型由聯合機率分布 P ( X , Y ) P(X,Y) P(X,Y) 表示，可以認為非監督學習訓練資料是聯合機率分布産生的資料。X為觀測資料，Y為未觀測資料。

參考：

1. 知乎 EM算法詳解 [EM算法原理總結 -
2. 劉建平Pinard - 部落格園](http://www.cnblogs.com/pinard/p/6912636.html)
3. 李航《統計機器學習》
4. https://blog.csdn.net/qq_34896915/article/details/75040578