最近看了很多關于EM算法推導的文章，包括《統計學習方法（第二版）》中關于EM算法的内容，總感覺說得不夠清楚。雖然公式都寫得挺詳細，但是沒有說清楚為什麼要這樣算，怎麼想到這樣變形的，本文總結一下我了解的EM算法推導過程，盡量把EM的核心思想說清楚。
- 參考：
- 《統計學習方法（第二版）》—— 第九章
- 淺析：從最大似然估計（MLE）、最大後驗估計（MAP）到期望最大算法（EM）
- 為什麼不把EM算法說清楚？
- EM算法詳細推導和講解
- 如何感性地了解EM算法？
- NLP —— 圖模型（零）：EM算法簡述及簡單示例（三硬币模型）
- EM算法的九重境界
文章目錄
- 0. 補充
- 1. 引入
-
- 1.1 三硬币模型
-
- 1.1.1 模組化
- 1.1.2 分析
- 1.1.2 使用EM方法解
- 1.2 問題形式化
- 2. EM算法
-
- 2.1 EM算法流程
- 2.2 EM算法的推導
-
- 2.2.1 EM核心思想：疊代優化
- 2.2.2 EM推導
-
- 2.2.2.1 兩個難點
- 2.2.2.2 解難點2：Jensen不等式
- 2.2.2.3 解難點1：分步疊代
- 2.3 EM算法的收斂性
-
- 2.3.1 定理1
- 2.3.2 定理2
- 2.3.3 小結
- 3. 用EM算法解三硬币問題
-
- 3.1 E步
-
- 3.1.1 公式法直接計算 Q 函數
- 3.1.2 從期望角度計算 Q 函數
- 3.2 M步
- 3.3 對比
- 3.4 代碼
0. 補充
- 極大似然估計、最大後驗估計：一文看懂 “極大似然估計” 與 “最大後驗估計”
- 關于Jensen不等式：Jensen不等式
1. 引入
- 根據觀測性，可以把模型變量分為
觀測變量
隐變量/潛在變量
- EM算法是一種疊代算法，用于對含有隐變量的機率模型參數進行極大似然估計或最大後驗機率估計。EM的适用情況和 MLE 或 MAP 一樣，都是在模型确定的情況下估計模型參數，差別在于EM算法允許模型中存在隐變量
1.1 三硬币模型
- 先看一個經典的三硬币模型示例，感受變量的觀測性和EM算法流程。
1.1.1 模組化
-
假設有 A,B,C 三枚硬币，正面出現的機率分别為 π , p , q \pi,p,q π,p,q。進行如下實驗：先擲硬币A，正面選B，反面選C；然後擲選出的硬币。我們不能觀測到A的結果（隐變量），但是可以觀測到選出硬币B或C的擲币結果（觀測變量）。假如10次實驗觀測結果為
正 ， 正 ， 反 ， 正 ， 反 ， 反 ， 正 ， 反 ， 正 ， 正 正，正，反，正，反，反，正，反，正，正 正，正，反，正，反，反，正，反，正，正
現在請估計三枚硬币擲出正面的機率，即三硬币模型的參數
-
用1表示正面，0表示反面，設随機變量 y , z y,z y,z 分别為觀測變量（B/C的擲币結果）和隐變量（A的擲币結果），觀測值分别為 Y , Z ∈ { 0 , 1 } Y,Z \in \{0,1\} Y,Z∈{0,1} ，模型參數為 θ = { π , p , q } \theta = \{\pi,p,q\} θ={π,p,q} 。注意到擲硬币是一種兩點分布，則三硬币問題可以模組化表示為
P ( y ∣ θ ) = ∑ z P ( y , z ∣ θ ) = ∑ z P ( z ∣ θ ) P ( y ∣ z , θ ) = π p y ( 1 − p ) 1 − y + ( 1 − π ) q y ( 1 − q ) 1 − y \begin{aligned} P(y|\theta) &= \sum_zP(y,z|\theta) = \sum_zP(z|\theta)P(y|z,\theta) \\ & = \pi p^y(1-p)^{1-y} + (1-\pi)q^y(1-q)^{1-y} \end{aligned} P(y∣θ)=z∑P(y,z∣θ)=z∑P(z∣θ)P(y∣z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y
- 将
觀測資料
未觀測資料
-
極大似然估計：觀測資料的似然函數、對數似然函數、參數 θ \theta θ 估計值分别為
P ( Y ∣ θ ) = ∑ Z P ( Z ∣ θ ) P ( Y ∣ Z , θ ) = ∏ j − 1 n [ π p Y j ( 1 − p ) 1 − Y j + ( 1 − π ) q Y j ( 1 − q ) 1 − Y j ] l o g P ( Y ∣ θ ) = ∑ j = 1 n l o g [ π p Y j ( 1 − p ) 1 − Y j + ( 1 − π ) q Y j ( 1 − q ) 1 − Y j ] \begin{aligned} P(Y|\theta) &= \sum_ZP(Z|\theta)P(Y|Z,\theta) \\ &= \prod_{j-1}^n[\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}] \\ logP(Y|\theta) &= \sum_{j=1}^nlog[\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}]\\ \end{aligned} P(Y∣θ)logP(Y∣θ)=Z∑P(Z∣θ)P(Y∣Z,θ)=j−1∏n[πpYj(1−p)1−Yj+(1−π)qYj(1−q)1−Yj]=j=1∑nlog[πpYj(1−p)1−Yj+(1−π)qYj(1−q)1−Yj]
θ ^ = arg max θ l o g P ( Y ∣ θ ) \hat{\theta} = \argmax_\theta logP(Y|\theta) θ^=θargmaxlogP(Y∣θ)
-
最大後驗估計：參數 θ \theta θ 的後驗機率、估計值分别為
P ( θ ∣ Y ) = P ( θ , Y ) P ( Y ) = P ( θ ) P ( Y ∣ θ ) P ( Y ) = P ( θ ) P ( Z ∣ θ ) P ( Y ∣ Z , θ ) P ( Y ) = P ( θ ) P ( Z ∣ θ ) P ( Y ∣ Z , θ ) ∫ P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( θ ) d θ ∝ P ( θ ) P ( Z ∣ θ ) P ( Y ∣ Z , θ ) \begin{aligned} P(\theta|Y) &= \frac{P(\theta,Y)}{P(Y)} \\ &= \frac{P(\theta)P(Y|\theta)}{P(Y)} \\ &= \frac{P(\theta)P(Z|\theta)P(Y|Z,\theta)}{P(Y)} \\ &= \frac{P(\theta)P(Z|\theta)P(Y|Z,\theta)}{\int P(Y|Z,\theta)P(Z|\theta)P(\theta)d\theta}\\ &\propto P(\theta)P(Z|\theta)P(Y|Z,\theta) \end{aligned} P(θ∣Y)=P(Y)P(θ,Y)=P(Y)P(θ)P(Y∣θ)=P(Y)P(θ)P(Z∣θ)P(Y∣Z,θ)=∫P(Y∣Z,θ)P(Z∣θ)P(θ)dθP(θ)P(Z∣θ)P(Y∣Z,θ)∝P(θ)P(Z∣θ)P(Y∣Z,θ)
θ ^ = arg max θ P ( θ ) P ( Z ∣ θ ) P ( Y ∣ Z , θ ) \hat{\theta} = \argmax_\theta P(\theta)P(Z|\theta)P(Y|Z,\theta) θ^=θargmaxP(θ)P(Z∣θ)P(Y∣Z,θ)
這裡需要 θ \theta θ 的先驗機率 p ( θ ) p(\theta) p(θ)，可以設成三個 μ = 0.5 \mu = 0.5 μ=0.5 的正态分布連乘，比較麻煩，是以這個例子後面就隻寫最大似然估計的過程了。
-
1.1.2 分析
-
先嘗試直接極大似然估計求解。令偏導為0列方程組
{ ∂ l o g P ( Y ∣ θ ) ∂ π = ∑ j = 1 n p y j ( 1 − p ) 1 − y j − q y j ( 1 − q ) 1 − y j π p y j ( 1 − p ) 1 − y j + ( 1 − π ) q y j ( 1 − q ) 1 − y j = 令 0 ∂ l o g P ( Y ∣ θ ) ∂ p = ∑ j = 1 n π p y j − 1 ( 1 − p ) − y j ( y j − p ) π p y j ( 1 − p ) 1 − y j + ( 1 − π ) q y j ( 1 − q ) 1 − y j = 令 0 ∂ l o g P ( Y ∣ θ ) ∂ q = ∑ j = 1 n ( 1 − π ) q y j − 1 ( 1 − q ) − y j ( y j − q ) π p y j ( 1 − p ) 1 − y j + ( 1 − π ) q y j ( 1 − q ) 1 − y j = 令 0 \left\{ \begin{aligned} \frac{\partial logP(Y|\theta)}{\partial \pi} & = \sum_{j=1}^n \frac{p^{y_j}(1-p)^{1-y_j}-q^{y_j}(1-q)^{1-y_j}}{\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi)q^{y_j}(1-q)^{1-y_j}} \stackrel{令}{=}0\\ \frac{\partial logP(Y|\theta)}{\partial p} & = \sum_{j=1}^n \frac{\pi p^{y_j-1}(1-p)^{-y_j}(y_j-p)}{\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi)q^{y_j}(1-q)^{1-y_j}} \stackrel{令}{=}0\\ \frac{\partial logP(Y|\theta)}{\partial q} & = \sum_{j=1}^n \frac{(1-\pi) q^{y_j-1}(1-q)^{-y_j}(y_j-q)}{\pi p^{y_j}(1-p)^{1-y_j} + (1-\pi)q^{y_j}(1-q)^{1-y_j}} \stackrel{令}{=}0 \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧∂π∂logP(Y∣θ)∂p∂logP(Y∣θ)∂q∂logP(Y∣θ)=j=1∑nπpyj(1−p)1−yj+(1−π)qyj(1−q)1−yjpyj(1−p)1−yj−qyj(1−q)1−yj=令0=j=1∑nπpyj(1−p)1−yj+(1−π)qyj(1−q)1−yjπpyj−1(1−p)−yj(yj−p)=令0=j=1∑nπpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj(1−π)qyj−1(1−q)−yj(yj−q)=令0
帶入觀測資料 Y Y Y，有
{ − 10 π q 2 − 10 π p 2 + 10 q 2 + 20 π p q − 10 p q + 6 p − 6 q = 0 ( 1 ) 10 π 2 q − 10 π 2 p − 10 π q + 6 π = 0 ( 2 ) − 10 π 2 q + 10 π 2 p + 20 π q − 10 π p − 10 q − 6 π + 6 = 0 ( 3 ) \left\{ \begin{aligned} &-10\pi q^2-10\pi p^2 + 10q^2+20\pi pq-10pq+6p-6q &= 0 \space\space\space\space (1)\\ &10\pi^2 q-10\pi^2 p-10\pi q+6\pi &= 0 \space\space\space\space(2) \\ &-10\pi^2 q+10\pi^2 p + 20\pi q -10\pi p - 10q -6\pi+6 &= 0 \space\space\space\space(3) \end{aligned} \right. ⎩⎪⎨⎪⎧−10πq2−10πp2+10q2+20πpq−10pq+6p−6q10π2q−10π2p−10πq+6π−10π2q+10π2p+20πq−10πp−10q−6π+6=0 (1)=0 (2)=0 (3)
發現式子(2)(3)有關系： ( 3 ) = π 1 − π ( 2 ) (3) =\frac{\pi}{1-\pi} (2) (3)=1−ππ(2)，這兩個式子等價，是以有無窮解。最終結果為 p = q = 0.6 p=q= 0.6 p=q=0.6， π \pi π 取任意值，這種結果是符合直覺的。由于未觀測資料 Z Z Z 不可見，直接最大化觀測資料 Y Y Y 的似然會導緻 Z Z Z 及相關參數被當作不存在，使得模型實質上退化到了一個硬币（不存在隐變量）的情況
1.1.2 使用EM方法解
- 使用EM算法可疊代求出更好的唯一解，過程如下
- 先選取參數初值， θ ( 0 ) = { π ( 0 ) , p ( 0 ) , q ( 0 ) } \theta^{(0)} = \{\pi^{(0)},p^{(0)},q^{(0)}\} θ(0)={π(0),p(0),q(0)}
- 反複執行以下兩步，疊代計算參數 θ \theta θ 估計值直到收斂為止。設第 i i i 次疊代參數估計值為 θ ( i ) = { π ( i ) , p ( i ) , q ( i ) } \theta^{(i)} = \{\pi^{(i)},p^{(i)},q^{(i)}\} θ(i)={π(i),p(i),q(i)}，則第 i + 1 i+1 i+1 次疊代如下
-
E步：計算模型在參數 θ ( i ) \theta^{(i)} θ(i) 下觀測資料 y j y_j yj 來自擲硬币 B 的機率 μ j ( i + 1 ) \mu_j^{(i+1)} μj(i+1)
μ j ( i + 1 ) = P ( B ∣ y j ) = P ( B , y j ) P ( y j ) = P ( B ) P ( y j ∣ B ) P ( B ) P ( y j ∣ B ) + P ( C ) P ( y j ∣ C ) = π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j + ( 1 − π ( i ) ) ( q ( i ) ) y j ( 1 − q ( i ) ) 1 − y j \begin{aligned} \mu_j^{(i+1)} = P(B|y_j) &= \frac{P(B,y_j)}{P(y_j)} = \frac{P(B)P(y_j|B)}{P(B)P(y_j|B)+P(C)P(y_j|C)} \\ &= \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}} \end{aligned} μj(i+1)=P(B∣yj)=P(yj)P(B,yj)=P(B)P(yj∣B)+P(C)P(yj∣C)P(B)P(yj∣B)=π(i)(p(i))yj(1−p(i))1−yj+(1−π(i))(q(i))yj(1−q(i))1−yjπ(i)(p(i))yj(1−p(i))1−yj
-
M步：計算模型參數的新估計值
π ( i + 1 ) = P ( B ∣ Y ) = 1 n ∑ i = 1 n μ j ( i + 1 ) p ( i + 1 ) = P ( B ∣ Y = 1 ) P ( B ∣ Y ) = ∑ i = 1 n μ j ( i + 1 ) y j ∑ i = 1 n μ j ( i + 1 ) q ( i + 1 ) = P ( C ∣ Y = 1 ) P ( C ∣ Y ) = ∑ i = 1 n ( 1 − μ j ( i + 1 ) ) y j ∑ i = 1 n ( 1 − μ j ( i + 1 ) ) \begin{aligned} \pi^{(i+1)} &= P(B|Y) = \frac{1}{n} \sum_{i=1}^n \mu_j^{(i+1)} \\ p^{(i+1)} &= \frac{P(B|Y=1)}{P(B|Y)} = \frac{\sum_{i=1}^n \mu_j^{(i+1)}y_j}{\sum_{i=1}^n \mu_j^{(i+1)}} \\ q^{(i+1)} &= \frac{P(C|Y=1)}{P(C|Y)} = \frac{\sum_{i=1}^n (1-\mu_j^{(i+1)})y_j}{\sum_{i=1}^n (1-\mu_j^{(i+1)})} \end{aligned} π(i+1)p(i+1)q(i+1)=P(B∣Y)=n1i=1∑nμj(i+1)=P(B∣Y)P(B∣Y=1)=∑i=1nμj(i+1)∑i=1nμj(i+1)yj=P(C∣Y)P(C∣Y=1)=∑i=1n(1−μj(i+1))∑i=1n(1−μj(i+1))yj
-
- 帶入觀測值計算
- 初值選 θ ( 0 ) = { 0.5 , 0.5 , 0.5 } \theta^{(0)} = \{0.5,0.5,0.5\} θ(0)={0.5,0.5,0.5} 時，EM方法解出的估計值為 θ ^ = { 0.5 , 0.6 , 0.6 } \hat{\theta} = \{0.5,0.6,0.6\} θ^={0.5,0.6,0.6}
- 初值選 θ ( 0 ) = { 0.4 , 0.6 , 0.7 } \theta^{(0)} = \{0.4,0.6,0.7\} θ(0)={0.4,0.6,0.7} 時，EM方法解出的估計值為 θ ^ = { 0.4064 , 0.5368 , 0.6432 } \hat{\theta} = \{0.4064,0.5368,0.6432\} θ^={0.4064,0.5368,0.6432}
1.2 問題形式化
- 一般地，用 Y Y Y 表示觀測變量的資料，用 Z Z Z 表示隐變量的資料，用 θ \theta θ 表示模型參數
-
完全資料
-
不完全資料
-
- 當完全資料已知時，可以直接使用MLE估計或MAP估計 θ \theta θ；隻知道不完全資料時，EM算法通過疊代最大化似然函數或後驗機率。每次疊代分成兩步
- E步：求Q函數（期望）
- M步：求Q函數（期望）極大
2. EM算法
2.1 EM算法流程
- 輸入：觀測變量資料 Y Y Y，隐變量資料 Z Z Z，聯合分布 P ( Y , Z ∣ θ ) P(Y,Z|\theta) P(Y,Z∣θ)，條件分布 P ( Z ∣ Y , θ ) P(Z|Y,\theta) P(Z∣Y,θ)
- 輸出：模型參數 θ \theta θ
- 流程：
- 選擇合适的初始參數 θ ( 0 ) \theta^{(0)} θ(0)，開始疊代。初值可以任意選擇，但需注意EM算法是對初值敏感的
-
E步：記 θ ( i ) \theta^{(i)} θ(i) 為第 i i i 輪疊代參數估計值，在第 i + 1 i+1 i+1 次疊代的 E 步，計算
Q ( θ , θ ( i ) ) = E Z [ l o g P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = P ( Z ∣ Y , θ ( i ) ) l o g P ( Y , Z ∣ θ ) \begin{aligned} Q(\theta,\theta^{(i)}) &= \mathbb{E}_Z[logP(Y,Z|\theta) | Y,\theta^{(i)}] \\ &= P(Z|Y,\theta^{(i)}) logP(Y,Z|\theta) \end{aligned} Q(θ,θ(i))=EZ[logP(Y,Z∣θ)∣Y,θ(i)]=P(Z∣Y,θ(i))logP(Y,Z∣θ)
這裡 P ( Z ∣ Y , θ ( i ) ) P(Z|Y,\theta^{(i)}) P(Z∣Y,θ(i)) 是給定觀測資料 Y Y Y 和目前估參數估計值 θ ( i ) \theta^{(i)} θ(i) 下隐變量資料 Z Z Z 的條件機率分布
-
M步：求使得 Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q(θ,θ(i)) 最大化的 θ \theta θ，确定第 i + 1 i+1 i+1 輪的參數估計值 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)
θ ( i + 1 ) = arg max θ Q ( θ , θ ( i ) ) \theta^{(i+1)} = \argmax_\theta Q(\theta,\theta^{(i)}) θ(i+1)=θargmaxQ(θ,θ(i))
-
重複第 2/3 步，直到收斂。收斂條件一般是：對于給定小量 ε 1 , ε 2 \varepsilon_1,\varepsilon_2 ε1,ε2，滿足下式則收斂，停止疊代
∣ θ ( i + 1 ) − θ ( i ) ∣ ∣ < ε 1 或 ∣ ∣ Q ( θ ( i + 1 ) , θ ( i ) ) − Q ( θ ( i ) , θ ( i ) ) ∣ ∣ < ε 1 |\theta^{(i+1)} - \theta^{(i)}|| < \varepsilon_1 \space\space或 \space\space ||Q(\theta^{(i+1)},\theta^{(i)}) - Q(\theta^{(i)},\theta^{(i)}) || < \varepsilon_1 \\ ∣θ(i+1)−θ(i)∣∣<ε1 或 ∣∣Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∣∣<ε1
-
Q函數
：定義完全資料的對數似然函數 l o g P ( Y , Z ∣ θ ) logP(Y,Z|\theta) logP(Y,Z∣θ) 關于給定觀測資料 Y Y Y 和目前估計參數 θ ( i ) \theta^{(i)} θ(i) 下對未觀測資料 Z Z Z 的條件機率分布 P ( Z ∣ Y , θ ( i ) ) P(Z|Y,\theta^{(i)}) P(Z∣Y,θ(i)) 的期望為Q函數，即
Q ( θ , θ ( i ) ) = E Z [ l o g P ( Y , Z ∣ θ ) ∣ Y , θ ( i ) ] = E Z ∣ Y , θ ( i ) l o g P ( Y , Z ∣ θ ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) l o g P ( Y , Z ∣ θ ) = ∑ j = 1 N ∑ Z P ( Z ∣ Y j , θ ( i ) ) l o g P ( Y j , Z ∣ θ ) \begin{aligned} Q(\theta,\theta^{(i)}) &= \mathbb{E}_Z[logP(Y,Z|\theta) | Y,\theta^{(i)}] \\ &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \\ &=\sum_Z P(Z|Y,\theta^{(i)})logP(Y,Z|\theta) \\ &= \sum_{j=1}^N \sum_Z P(Z|Y_j,\theta^{(i)})logP(Y_j,Z|\theta) \end{aligned} Q(θ,θ(i))=EZ[logP(Y,Z∣θ)∣Y,θ(i)]=EZ∣Y,θ(i)logP(Y,Z∣θ)=Z∑P(Z∣Y,θ(i))logP(Y,Z∣θ)=j=1∑NZ∑P(Z∣Yj,θ(i))logP(Yj,Z∣θ)
EM 算法的每次疊代，實際是在求Q函數這個期望及其極大，這也是EM算法全程 “期望最大化”（Expection-Maximization）的由來
2.2 EM算法的推導
- EM算法目标：極大化參數 θ \theta θ 下觀測資料 Y Y Y 的對數似然函數（或最大後驗機率）
- 方法：利用Jensen不等式推出目标對數似然函數下界，疊代優化下界以優化目标，不保證全局最優
2.2.1 EM核心思想：疊代優化
-
在推導之前，有必要探讨一下為什麼直接最大化似然函數沒有解析解。先看極大化目标：參數 θ \theta θ 下觀測資料 Y Y Y 的對數似然函數。即
L ( θ ) = l o g P ( Y ∣ θ ) = l o g ∑ Z P ( Y , Z ∣ θ ) = l o g ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) \begin{aligned} L(\theta) &= logP(Y|\theta) = log\sum_Z P(Y,Z|\theta) \\ &= log(\sum_ZP(Y|Z,\theta)P(Z|\theta)) \end{aligned} L(θ)=logP(Y∣θ)=logZ∑P(Y,Z∣θ)=log(Z∑P(Y∣Z,θ)P(Z∣θ))
回顧極大似然估計 MLE 的一般形式，假設可以觀測到完全資料 X X X，對數似然函數為
L ( θ ) = l o g ∏ X P ( X ∣ θ ) = ∑ X l o g P ( X ∣ θ ) L(\theta) = log \prod_X P(X|\theta) = \sum_X log P(X|\theta) L(θ)=logX∏P(X∣θ)=X∑logP(X∣θ)
要得到 arg max θ L ( θ ) \argmax_\theta L(\theta) θargmaxL(θ) 的解析解，需要保證 L ( θ ) L(\theta) L(θ) 中隻有 θ \theta θ 一個變量。在隐變量 z z z 不能觀測時，優化目标中同時存在估計量 P ( Z ∣ θ ) P(Z|\theta) P(Z∣θ) 和 θ \theta θ 兩個變量，直接算自然會得到無窮多解。
- 事實上，未觀測資料 Z Z Z 和參數 θ \theta θ 是互相依賴的，估計一個需要借助另一個，更新一個的估計值會導緻另一個的估計值也變化，這就導緻了一個先有雞還是先有蛋的問題。還是以三硬币模型為例，如果我們知道 Z Z Z，就能知道哪些觀測結果來自硬币B，哪些來自硬币C，進而能分别使用極大似然估計法得到 p , q p,q p,q。在隻有不完全資料 Y Y Y 時，為了區分觀測結果的來源以估計 p , q p,q p,q，我們需要先估計每一個觀測資料 Y i Y_i Yi 對應的 Z i Z_i Zi 取值（或取值的分布），這個估計又是依賴于 p , q p,q p,q 的，如圖所示 對于這種含有兩個未知變量的優化問題，可以使用疊代的方法。選取一個合适的初始值開始疊代，每輪疊代過程中，先固定一個變量優化另一個，再固定另一個變量估計優化這個，直到二者逐漸收斂為止。EM 算法在此思想基礎上更進一步，在估計未觀測資料時，不止估計其值，而是估計其分布，進而把優化目标轉換完全資料似然函數為對未觀測資料分布的期望形式。此處推薦參考一篇文章：如何感性地了解EM算法？
EM算法（1）—— 徹底了解EM算法推導0. 補充1. 引入2. EM算法3. 用EM算法解三硬币問題 - 強化學習中 model-based 方法 policy-iteration 使用了這種分步疊代思想，不妨把它的示意圖貼過來 policy-iteration 的疊代分為兩個階段，evaluation步驟對目前政策下的狀态價值進行評估，improvement步驟利用政策價值貪心地優化政策。可以把每個流程看成二維空間中的一條線，代表對于其目标的解決方案。每個流程都把價值 V V V 或 政策 π \pi π 推向其中一條線，由于這兩條線不是正交的，是以兩個目标間會産生互相作用。evaluation步會導緻政策逐漸偏離最優，improvement步會導緻價值偏離新政策，直接沖向某個目标會導緻某種程度上偏離另一個目标，然而整體流程的結果最終會更接近最優的總目标。EM算法中的E步可以類比為evaluation步驟，M步可以類比為improvement步驟。
EM算法（1）—— 徹底了解EM算法推導0. 補充1. 引入2. EM算法3. 用EM算法解三硬币問題
2.2.2 EM推導
2.2.2.1 兩個難點
-
極大化目标：參數 θ \theta θ 下觀測資料 Y Y Y 的對數似然函數。即
L ( θ ) = l o g P ( Y ∣ θ ) = l o g ∑ Z P ( Y , Z ∣ θ ) = l o g ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) \begin{aligned} L(\theta) &= logP(Y|\theta) = log\sum_Z P(Y,Z|\theta) \\ &= log(\sum_ZP(Y|Z,\theta)P(Z|\theta)) \\ \end{aligned} L(θ)=logP(Y∣θ)=logZ∑P(Y,Z∣θ)=log(Z∑P(Y∣Z,θ)P(Z∣θ))
- 注意到 θ ^ = arg max θ L ( θ ) \hat{\theta} = \argmax_\theta L(\theta) θ^=θargmaxL(θ) 這一極大化的主要困難有兩個
- 需要同時完成對 P ( Z ∣ θ ) P(Z|\theta) P(Z∣θ) 的估計和 θ \theta θ 的優化，但兩個目标會互相影響
- 求和符号 ∑ \sum ∑ 在 l o g log log 裡，沒法進一步變形處理
2.2.2.2 解難點2：Jensen不等式
-
首先處理第二個難點，我們希望能把 ∑ \sum ∑ 符号拿到 l o g log log 外面，這時就要用到 jensen 不等式，即
∑ i = 1 M λ i f ( x i ) ≥ f ( ∑ i = 1 M λ i x i ) \sum_{i=1}^M \lambda_if(x_i) \geq f(\sum_{i=1}^M\lambda_ix_i) i=1∑Mλif(xi)≥f(i=1∑Mλixi)
此式要求 f f f 是凸函數， λ i ≥ 0 \lambda_i\geq0 λi≥0 且 ∑ i = 1 M λ i = 1 \sum_{i=1}^M \lambda_i = 1 ∑i=1Mλi=1，當随機變量 X X X 是常數時等号成立。由于 l o g log log 是凹函數，是以需要應用 jensen 不等式的凹函數形式，即把 ≥ \geq ≥ 換成 ≤ \leq ≤。
-
為了應用 Jensen 不等式，還要拼湊一個 λ \lambda λ。注意到 l o g log log 裡的和是關于隐變量 Z Z Z 的和，而 λ \lambda λ 需要滿足 λ i ≥ 0 \lambda_i\geq0 λi≥0 且 ∑ i = 1 M λ i = 1 \sum_{i=1}^M \lambda_i = 1 ∑i=1Mλi=1，不妨設 λ \lambda λ 為隐變量 z z z 的分布函數 K ( Z ∣ θ ) K(Z|\theta) K(Z∣θ)，于是有
L ( θ ) = l o g P ( Y ∣ θ ) = l o g ( ∑ Z P ( Y , Z ∣ θ ) ) = l o g ( ∑ Z K ( Z ∣ θ ) P ( Y , Z ∣ θ ) K ( Z ∣ θ ) ) ≥ ∑ Z K ( Z ∣ θ ) l o g P ( Y , Z ∣ θ ) K ( Z ∣ θ ) \begin{aligned} L(\theta) &= logP(Y|\theta) \\ &=log(\sum_Z P(Y,Z|\theta)) \\ &= log(\sum_Z K(Z|\theta) \frac{P(Y,Z|\theta)}{K(Z|\theta)}) \\ &\geq \sum_Z K(Z|\theta)log\frac{P(Y,Z|\theta)}{K(Z|\theta)} \end{aligned} L(θ)=logP(Y∣θ)=log(Z∑P(Y,Z∣θ))=log(Z∑K(Z∣θ)K(Z∣θ)P(Y,Z∣θ))≥Z∑K(Z∣θ)logK(Z∣θ)P(Y,Z∣θ)
這樣我們就找到了似然函數的下界形式，并成功把 ∑ \sum ∑ 拿到了 l o g log log 外邊，接下來我們希望通過最大化下界來間接地最大化 L ( θ ) L(\theta) L(θ)。為了保證優化效果，我們要求 ∑ Z K ( Z ∣ θ ) l o g ( P ( Y , Z ∣ θ ) K ( Z ∣ θ ) ) \sum_Z K(Z|\theta)log(\frac{P(Y,Z|\theta)}{K(Z|\theta)}) ∑ZK(Z∣θ)log(K(Z∣θ)P(Y,Z∣θ)) 是一個緊的下界，也就是至少有一點使不等式中等号成立。因為如果這個下界特别小，離似然函數太遠，優化它也沒有什麼用
-
現在來讨論 K ( Z ∣ θ ) K(Z|\theta) K(Z∣θ) 的形式，為了得到緊下界，回頭看 Jensen 不等式中使得等号成立的條件，需要滿足
P ( Y , Z ∣ θ ) K ( Z ∣ θ ) = C \frac{P(Y,Z|\theta)}{K(Z|\theta)} = C K(Z∣θ)P(Y,Z∣θ)=C
因為 K ( Z ∣ θ ) K(Z|\theta) K(Z∣θ) 是 Z Z Z 的分布函數，有
∵ ∑ Z K ( Z ∣ θ ) = ∑ Z P ( Y , Z ∣ θ ) C = 1 ∴ C = ∑ Z P ( Y , Z ∣ θ ) = P ( Y ∣ θ ) ∴ K ( Z ∣ θ ) = P ( Y , Z ∣ θ ) C = P ( Y , Z ∣ θ ) P ( Y ∣ θ ) = P ( Z ∣ Y , θ ) \begin{aligned} &\because \sum_ZK(Z|\theta) = \sum_Z \frac{P(Y,Z|\theta)}{C} = 1 \\ &\therefore C = \sum_Z P(Y,Z|\theta) = P(Y|\theta) \\ &\therefore K(Z|\theta) = \frac{P(Y,Z|\theta)}{C} = \frac{P(Y,Z|\theta)}{P(Y|\theta)} = P(Z|Y,\theta) \end{aligned} ∵Z∑K(Z∣θ)=Z∑CP(Y,Z∣θ)=1∴C=Z∑P(Y,Z∣θ)=P(Y∣θ)∴K(Z∣θ)=CP(Y,Z∣θ)=P(Y∣θ)P(Y,Z∣θ)=P(Z∣Y,θ)
-
這樣我們就得到了緊下界的表示形式，由于始終滿足等号成立條件 P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ) = P ( Y ∣ θ ) ≡ C \frac{P(Y,Z|\theta)}{P(Z|Y,\theta)} = P(Y|\theta) \equiv C P(Z∣Y,θ)P(Y,Z∣θ)=P(Y∣θ)≡C，是以不等号轉換為恒等号，可以證明如下
∑ Z P ( Z ∣ Y , θ ) l o g P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ) = ∑ Z P ( Z ∣ Y , θ ) l o g P ( Y ∣ θ ) = l o g P ( Y ∣ θ ) = L ( θ ) \begin{aligned} \sum_Z P(Z|Y,\theta)log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta)} &= \sum_Z P(Z|Y,\theta)logP(Y|\theta) \\ &= logP(Y|\theta) \\ &= L(\theta) \end{aligned} Z∑P(Z∣Y,θ)logP(Z∣Y,θ)P(Y,Z∣θ)=Z∑P(Z∣Y,θ)logP(Y∣θ)=logP(Y∣θ)=L(θ)
2.2.2.3 解難點1：分步疊代
- 難點1的本質就是沒法同時估計 P ( Z ∣ Y , θ ) P(Z|Y,\theta) P(Z∣Y,θ) 和優化 θ \theta θ ，我們上面的一堆分析并沒有解決這個問題。牢記 2.2.1 節的分析，解決這種問題可以使用分步疊代法，每輪疊代先控制一個變量優化另一個變量，再控制另一個變量優化這個變量。具體到EM算法中，第 i i i 輪疊代分兩步
- 控制 θ \theta θ 估計 Z Z Z：固定模型參數 θ ( i ) \theta^{(i)} θ(i)，估計 Z Z Z 條件機率分布 P ( Z ∣ Y , θ ( i ) ) P(Z|Y,\theta^{(i)}) P(Z∣Y,θ(i))
-
控制 Z Z Z 優化 θ \theta θ：利用估計的 P ( Z ∣ Y , θ ( i ) ) P(Z|Y,\theta^{(i)}) P(Z∣Y,θ(i)) 重新構造緊下界函數，最大化緊下界函數得到新的參數值 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)。由于
P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) 不 恒 等 于 C \frac{P(Y,Z|\theta)}{P(Z|Y,\theta^{(i)})} 不恒等于 C P(Z∣Y,θ(i))P(Y,Z∣θ)不恒等于C
此時等号恒成立條件不再滿足，還原到不等号
L ( θ ) = l o g ( ∑ Z P ( Y , Z ∣ θ ) ) ≥ ∑ Z P ( Z ∣ Y , θ ( i ) ) l o g P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) L(\theta) =log(\sum_ZP(Y,Z|\theta)) \geq \sum_Z P(Z|Y,\theta^{(i)}) log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^{(i)})} L(θ)=log(Z∑P(Y,Z∣θ))≥Z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y,Z∣θ)
不妨設緊下界為 l ( θ , θ ( i ) ) l(\theta,\theta^{(i)}) l(θ,θ(i))，即
l ( θ , θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) l o g P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) L ( θ ) ≥ l ( θ , θ ( i ) ) l(\theta,\theta^{(i)}) = \sum_Z P(Z|Y,\theta^{(i)}) log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^{(i)})} \\ L(\theta) \geq l(\theta,\theta^{(i)}) l(θ,θ(i))=Z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y,Z∣θ)L(θ)≥l(θ,θ(i))
注意到當 θ = θ ( i ) \theta = \theta^{(i)} θ=θ(i) 時， L ( θ ( i ) ) = l ( θ ( i ) , θ ( i ) ) L(\theta^{(i)}) =l(\theta^{(i)},\theta^{(i)}) L(θ(i))=l(θ(i),θ(i)) 等号成立
如圖所示，求 θ ( i + 1 ) = arg max θ l ( θ , θ ( i ) ) \theta^{(i+1)} = \argmax_\theta l(\theta,\theta^{(i)}) θ(i+1)=θargmaxl(θ,θ(i)) ，就可以通過增大緊下界間接地增大 L ( θ ) L(\theta) L(θ) EM算法（1）—— 徹底了解EM算法推導0. 補充1. 引入2. EM算法3. 用EM算法解三硬币問題
-
忽略掉對 l ( θ , θ ( i ) ) l(\theta,\theta^{(i)}) l(θ,θ(i)) 中對極大化而言無關的項，得到 Q函數如下
Q ( θ , θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) l o g P ( Y , Z ∣ θ ) = E Z ∣ Y , θ ( i ) l o g P ( Y , Z ∣ θ ) \begin{aligned} Q(\theta,\theta^{(i)}) &= \sum_Z P(Z|Y,\theta^{(i)}) logP(Y,Z|\theta) \\ &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \end{aligned} Q(θ,θ(i))=Z∑P(Z∣Y,θ(i))logP(Y,Z∣θ)=EZ∣Y,θ(i)logP(Y,Z∣θ)
-
大功告成！至此我們已經對EM算法流程有了清晰的了解，并完整地推導了EM的核心函數 Q Q Q。不妨比較一下 Q Q Q 函數和原始的優化目标 L ( θ ) L(\theta) L(θ)
Q ( θ , θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) l o g P ( Y , Z ∣ θ ) = E Z ∣ Y , θ ( i ) l o g ( P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) = E Z ∣ Y , θ ( i ) l o g P ( Y , Z ∣ θ ) L ( θ ) = l o g ( ∑ Z P ( Y , Z ∣ θ ) ) = l o g ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) = l o g ( E Z ∣ θ P ( Y ∣ Z , θ ) ) \begin{aligned} Q(\theta,\theta^{(i)}) &= \sum_Z P(Z|Y,\theta^{(i)}) logP(Y,Z|\theta) \\ &= \mathbb{E}_{Z|Y,\theta^{(i)}}log(P(Y|Z,\theta)P(Z|\theta)) \\ &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \\ L(\theta) &= log(\sum_ZP(Y,Z|\theta)) \\ &= log(\sum_ZP(Y|Z,\theta)P(Z|\theta))\\ &= log(\mathbb{E}_{Z|\theta}P(Y|Z,\theta)) \end{aligned} Q(θ,θ(i))L(θ)=Z∑P(Z∣Y,θ(i))logP(Y,Z∣θ)=EZ∣Y,θ(i)log(P(Y∣Z,θ)P(Z∣θ))=EZ∣Y,θ(i)logP(Y,Z∣θ)=log(Z∑P(Y,Z∣θ))=log(Z∑P(Y∣Z,θ)P(Z∣θ))=log(EZ∣θP(Y∣Z,θ))
主要的變化有兩個
- 使用 Jensen 不等式把 ∑ \sum ∑ 拿到 l o g log log 外邊
- 引入後驗機率 P ( Z ∣ Y , θ ( i ) ) P(Z|Y,\theta^{(i)}) P(Z∣Y,θ(i))，使疊代求解成為可能。 這裡可能會有疑惑， Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q(θ,θ(i)) 裡有個 P ( Y , Z ∣ θ ) = P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P(Y,Z|\theta) = P(Y|Z,\theta)P(Z|\theta) P(Y,Z∣θ)=P(Y∣Z,θ)P(Z∣θ)，這不還是要用 θ \theta θ 估計未觀測資料分布 P ( Z ∣ θ ) P(Z|\theta) P(Z∣θ) 嗎。我個人了解是：式子中存在 P ( Z ∣ θ ) P(Z|\theta) P(Z∣θ) 不是關鍵，關鍵要看優化目标是關于什麼的期望，或者說優化目标這個随機變量服從的分布是什麼。
2.3 EM算法的收斂性
- EM算法提供了一種近似計算含有隐變量機率模型的極大似然估計的方法。其最大優點是簡單性和普适性。
- EM算法得到的估計序列是否收斂，如果收斂，是否收斂到全局最大值或局部最大值？下面給出關于EM算法收斂性的兩個定理。證明見《統計學習方法（第二版）》第九章
2.3.1 定理1
-
設 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ) 為觀測資料的似然函數， θ ( i ) ( i = 1 , 2 , . . . ) \theta^{(i)}(i=1,2,...) θ(i)(i=1,2,...) 為EM算法得到的參數估計序列， P ( Y ∣ θ ( i ) ) ( i = 1 , 2 , . . . ) P(Y|\theta^{(i)})(i=1,2,...) P(Y∣θ(i))(i=1,2,...) 為對應的似然函數序列，則 P ( Y ∣ θ ( i ) ) P(Y|\theta^{(i)}) P(Y∣θ(i)) 是單調遞增的，即
P ( Y ∣ θ ( i + 1 ) ) ≥ P ( Y ∣ θ ( i ) ) P(Y|\theta^{(i+1)}) \geq P(Y|\theta^{(i)}) P(Y∣θ(i+1))≥P(Y∣θ(i))
2.3.2 定理2
- 設 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ) 為觀測資料的似然函數， θ ( i ) ( i = 1 , 2 , . . . ) \theta^{(i)}(i=1,2,...) θ(i)(i=1,2,...) 為EM算法得到的參數估計序列， P ( Y ∣ θ ( i ) ) ( i = 1 , 2 , . . . ) P(Y|\theta^{(i)})(i=1,2,...) P(Y∣θ(i))(i=1,2,...) 為對應的似然函數序列， L ( θ ( i ) ) L(\theta^{(i)}) L(θ(i)) 為對應的對數似然函數序列，則
- 如果 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ) 有上界，則 L ( θ ( i ) ) = l o g P ( Y ∣ θ ( i ) ) L(\theta^{(i)}) = logP(Y|\theta^{(i)}) L(θ(i))=logP(Y∣θ(i)) 收斂到某一值 L ∗ L^* L∗
- 在函數 Q ( θ , θ ′ ) Q(\theta,\theta') Q(θ,θ′) 與 L ( θ ) L(\theta) L(θ) 滿足一定條件下（大多數情況下都滿足），由 EM 算法得到的參數估計序列 θ ( i ) \theta^{(i)} θ(i) 的收斂值 θ ∗ \theta^* θ∗ 是 L ( θ ) L(\theta) L(θ) 的穩定點。
- 參考文獻：Wu C . On the Convergence Properties of the EM Algorithm[J]. Annals of Statistics, 1983, 11(1):95-103.
2.3.3 小結
- EM算法的收斂性包含兩層意思，前者不蘊含後者
- 對數似然函數 L ( θ ( i ) ) L(\theta^{(i)}) L(θ(i)) 的收斂性
- 參數估計序列 θ ( i ) \theta^{(i)} θ(i) 的收斂性
- 定理僅保證參數估計序列 θ ( i ) \theta^{(i)} θ(i) 收斂到對數似然函數序列 L ( θ ( i ) ) L(\theta^{(i)}) L(θ(i)) 的穩定點（就是駐點，導數為0的點），不能保證收斂到極大值點。如果 L ( θ ( i ) ) L(\theta^{(i)}) L(θ(i)) 是凹的，則EM算法可以保證收斂到全局極大值，這點和梯度下降法這樣的疊代算法相同。
- 在應用中，初值的選擇變得非常重要，常用的辦法是選取幾個不同的初值進行疊代，對各個初始值加以比較，從中選擇最好的。
3. 用EM算法解三硬币問題
- 了解EM算法原理之後，回看 1.1節 最後的三硬币問題求解過程，分析M步和E步的計算
- 把三硬币問題的 “觀測資料” 和 “未觀測資料” 分别表示為來自總體 y y y 和 z z z 的容量為 n n n 的簡單随機樣本 y = { y 1 , y 2 , . . . , y n } T , z = { z 1 , z 2 , . . . , z n } T y = \{y_1,y_2,...,y_n\}^T,z = \{z_1,z_2,...,z_n\}^T y={y1,y2,...,yn}T,z={z1,z2,...,zn}T ，觀測值和未觀測值分别為 Y = { Y 1 , Y 2 , . . . , Y n } T , Z = { Z 1 , Z 2 , . . . , Z n } T Y = \{Y_1,Y_2,...,Y_n\}^T,Z = \{Z_1,Z_2,...,Z_n\}^T Y={Y1,Y2,...,Yn}T,Z={Z1,Z2,...,Zn}T，考慮求模型參數 θ \theta θ
- 注意： y i , Y i y_i,Y_i yi,Yi 和 z i , Z i z_i,Z_i zi,Zi 都是随機變量及其觀測值的關系，式子中大部分時候等價。下面這部分規範一下符号，約定
- 在帶入觀測值前（計算Q函數前），機率相關計算一律使用 y i , z i y_i,z_i yi,zi；帶入觀測值後（計算Q函數時）一律使用 Y i , Z i Y_i,Z_i Yi,Zi
- 求和符号 ∑ j = 1 N \sum_{j=1}^N ∑j=1N 是對每個随機變量 y j y_j yj的取值或觀測值 Y j Y_j Yj求和； ∑ z \sum_z ∑z 是對隐變量的可能取值（0或1）求和
3.1 E步
-
先求 P ( z ∣ y j , θ ( i ) ) P(z|y_j,\theta^{(i)}) P(z∣yj,θ(i))
P ( z ∣ y j , θ ( i ) ) = { μ j ( i + 1 ) = π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j π ( i ) ( p ( i ) ) y j ( 1 − p ( i ) ) 1 − y j + ( 1 − π ( i ) ) ( q ( i ) ) y j ( 1 − q ( i ) ) 1 − y j , i f z = 1 1 − μ j ( i + 1 ) , i f z = 0 P(z|y_j,\theta^{(i)}) = \left\{ \begin{aligned} &\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}} & ,if \space z = 1\\ &1-\mu_j^{(i+1)} &,if \space z=0 \end{aligned} \right. P(z∣yj,θ(i))=⎩⎪⎨⎪⎧μj(i+1)=π(i)(p(i))yj(1−p(i))1−yj+(1−π(i))(q(i))yj(1−q(i))1−yjπ(i)(p(i))yj(1−p(i))1−yj1−μj(i+1),if z=1,if z=0
-
再求 P ( z , y j ∣ θ ) = P ( y j ∣ z , θ ) P ( z ∣ θ ) P(z,y_j|\theta) = P(y_j|z,\theta)P(z|\theta) P(z,yj∣θ)=P(yj∣z,θ)P(z∣θ)
P ( z , y j ∣ θ ) = { π p y j ( 1 − p ) 1 − y j , i f z = 1 ( 1 − π ) q y j ( 1 − q ) 1 − y j , i f z = 0 P(z,y_j|\theta) = \left\{ \begin{aligned} &\pi p^{y_j} (1-p)^{1-y_j} & ,if \space z = 1 \\ &(1-\pi) q^{y_j} (1-q)^{1-y_j} & ,if \space z = 0 \end{aligned} \right. P(z,yj∣θ)={πpyj(1−p)1−yj(1−π)qyj(1−q)1−yj,if z=1,if z=0
3.1.1 公式法直接計算 Q 函數
-
直接帶入公式計算 Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q(θ,θ(i))
Q ( θ , θ ( i ) ) = E Z ∣ Y , θ ( i ) l o g P ( Y , Z ∣ θ ) = ∑ j = 1 N E Z ∣ Y j , θ ( i ) l o g P ( Y j , Z ∣ θ ) = ∑ j = 1 N { ∑ Z P ( Z ∣ Y j , θ ( i ) ) l o g P ( Y j , Z ∣ θ ) } = ∑ j = 1 N { P ( Z = 1 ∣ Y j , θ ( i ) ) l o g ( Y j , Z = 1 ∣ θ ) + P ( Z = 0 ∣ Y j , θ ( i ) ) l o g ( Y j , Z = 0 ∣ θ ) } = ∑ j = 1 N { μ j ( i + 1 ) l o g [ π p Y j ( 1 − p ) 1 − Y j ] + ( 1 − μ j ( i + 1 ) ) l o g [ ( 1 − π ) q Y j ( 1 − q ) 1 − Y j ] } \begin{aligned} Q(\theta,\theta^{(i)}) &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \\ &= \sum_{j=1}^N \mathbb{E}_{Z|Y_j,\theta^{(i)}}logP(Y_j,Z|\theta) \\ &= \sum_{j=1}^N\{\sum_ZP(Z|Y_j,\theta^{(i)})logP(Y_j,Z|\theta)\} \\ &= \sum_{j=1}^N\{P(Z=1|Y_j,\theta^{(i)})log(Y_j,Z=1|\theta) + P(Z=0|Y_j,\theta^{(i)})log(Y_j,Z=0|\theta)\} \\ &= \sum_{j=1}^N\{\mu_j^{(i+1)}log[\pi p^{Y_j} (1-p)^{1-Y_j}] + (1-\mu_j^{(i+1)})log[(1-\pi) q^{Y_j} (1-q)^{1-Y_j}]\} \end{aligned} Q(θ,θ(i))=EZ∣Y,θ(i)logP(Y,Z∣θ)=j=1∑NEZ∣Yj,θ(i)logP(Yj,Z∣θ)=j=1∑N{Z∑P(Z∣Yj,θ(i))logP(Yj,Z∣θ)}=j=1∑N{P(Z=1∣Yj,θ(i))log(Yj,Z=1∣θ)+P(Z=0∣Yj,θ(i))log(Yj,Z=0∣θ)}=j=1∑N{μj(i+1)log[πpYj(1−p)1−Yj]+(1−μj(i+1))log[(1−π)qYj(1−q)1−Yj]}
3.1.2 從期望角度計算 Q 函數
-
從期望角度計算
∵ P ( y j , z ∣ θ ) = { π p y j ( 1 − p ) 1 − y j , i f z = 1 ( 1 − π ) q y j ( 1 − q ) 1 − y j , i f z = 0 ∴ l o g P ( y j , z ∣ θ ) = { l o g [ π p y j ( 1 − p ) 1 − y j ] , i f z = 1 l o g [ ( 1 − π ) q y j ( 1 − q ) 1 − y j ] , i f z = 0 ∴ E z ∣ y i , θ ( i ) l o g P ( y j , z ∣ θ ) = { E z = 1 ∣ y i , θ ( i ) l o g P ( y j , z ∣ θ ) , i f z = 1 E z = 0 ∣ y i , θ ( i ) l o g P ( y j , z ∣ θ ) , i f z = 0 = { P ( z = 1 ∣ y i , θ ( i ) ) l o g [ π p y j ( 1 − p ) 1 − y j ] , i f z = 1 P ( z = 0 ∣ y i , θ ( i ) ) l o g [ ( 1 − π ) q y j ( 1 − q ) 1 − y j ] , i f z = 0 = { μ j ( i + 1 ) l o g [ π p y j ( 1 − p ) 1 − y j ] , i f z = 1 ( 1 − μ j ( i + 1 ) ) l o g [ ( 1 − π ) q y j ( 1 − q ) 1 − y j ] , i f z = 0 ∴ Q ( θ , θ ( i ) ) = E Z ∣ Y , θ ( i ) l o g P ( Y , Z ∣ θ ) = ∑ j = 1 N { E Z = 0 ∣ Y j , θ ( i ) l o g P ( Y j , Z ∣ θ ) + E Z = 1 ∣ Y j , θ ( i ) l o g P ( Y j , Z ∣ θ ) } = ∑ j = 1 N { μ j ( i + 1 ) l o g [ π p Y j ( 1 − p ) 1 − Y j ] + ( 1 − μ j ( i + 1 ) ) l o g [ ( 1 − π ) q Y j ( 1 − q ) 1 − Y j ] } \begin{aligned} \because P(y_j,z|\theta) &= \left\{ \begin{aligned} &\pi p^{y_j} (1-p)^{1-y_j} & ,if \space z= 1 \\ &(1-\pi) q^{y_j} (1-q)^{1-y_j} & ,if \space z= 0 \end{aligned} \right. \\ \therefore logP(y_j,z|\theta) &= \left\{ \begin{aligned} &log[\pi p^{y_j} (1-p)^{1-y_j}] & ,if \space z= 1 \\ &log[(1-\pi) q^{y_j} (1-q)^{1-y_j}] & ,if \space z= 0 \end{aligned} \right. \\ \therefore \mathbb{E}_{z|y_i,\theta^{(i)}}logP(y_j,z|\theta) &= \left\{ \begin{aligned} & \mathbb{E}_{z=1|y_i,\theta^{(i)}}logP(y_j,z|\theta) & ,if \space z= 1 \\ & \mathbb{E}_{z=0|y_i,\theta^{(i)}}logP(y_j,z|\theta) & ,if \space z= 0 \end{aligned} \right. \\ &= \left\{ \begin{aligned} & P(z=1|y_i,\theta^{(i)}) log[\pi p^{y_j} (1-p)^{1-y_j}] & ,if \space z= 1 \\ & P(z=0|y_i,\theta^{(i)}) log[(1-\pi) q^{y_j} (1-q)^{1-y_j}] & ,if \space z= 0 \end{aligned} \right. \\ &= \left\{ \begin{aligned} & \mu_j^{(i+1)} log[\pi p^{y_j} (1-p)^{1-y_j}] & ,if \space z= 1 \\ & (1-\mu_j^{(i+1)}) log[(1-\pi) q^{y_j} (1-q)^{1-y_j}] & ,if \space z= 0 \end{aligned} \right. \\ \therefore Q(\theta,\theta^{(i)}) &= \mathbb{E}_{Z|Y,\theta^{(i)}}logP(Y,Z|\theta) \\ &= \sum_{j=1}^N\{\mathbb{E}_{Z=0|Y_j,\theta^{(i)}}logP(Y_j,Z|\theta)+ \mathbb{E}_{Z=1|Y_j,\theta^{(i)}}logP(Y_j,Z|\theta)\}\\ &= \sum_{j=1}^N\{\mu_j^{(i+1)}log[\pi p^{Y_j} (1-p)^{1-Y_j}] + (1-\mu_j^{(i+1)})log[(1-\pi) q^{Y_j} (1-q)^{1-Y_j}]\} \end{aligned} ∵P(yj,z∣θ)∴logP(yj,z∣θ)∴Ez∣yi,θ(i)logP(yj,z∣θ)∴Q(θ,θ(i))={πpyj(1−p)1−yj(1−π)qyj(1−q)1−yj,if z=1,if z=0={log[πpyj(1−p)1−yj]log[(1−π)qyj(1−q)1−yj],if z=1,if z=0={Ez=1∣yi,θ(i)logP(yj,z∣θ)Ez=0∣yi,θ(i)logP(yj,z∣θ),if z=1,if z=0={P(z=1∣yi,θ(i))log[πpyj(1−p)1−yj]P(z=0∣yi,θ(i))log[(1−π)qyj(1−q)1−yj],if z=1,if z=0=⎩⎨⎧μj(i+1)log[πpyj(1−p)1−yj](1−μj(i+1))log[(1−π)qyj(1−q)1−yj],if z=1,if z=0=EZ∣Y,θ(i)logP(Y,Z∣θ)=j=1∑N{EZ=0∣Yj,θ(i)logP(Yj,Z∣θ)+EZ=1∣Yj,θ(i)logP(Yj,Z∣θ)}=j=1∑N{μj(i+1)log[πpYj(1−p)1−Yj]+(1−μj(i+1))log[(1−π)qYj(1−q)1−Yj]}
3.2 M步
-
對 Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q(θ,θ(i)) 求偏導
∂ Q ( θ , θ ( i ) ) ∂ π = ∑ j = 1 N { μ j ( i + 1 ) p Y j ( 1 − p ) 1 − Y j π p Y j ( 1 − p ) 1 − Y j + ( 1 − μ j ( i + 1 ) ) − q Y j ( 1 − q ) 1 − Y j ( 1 − π ) q Y j ( 1 − q ) 1 − Y j } = ∑ j = 1 N { μ j ( i + 1 ) − π π ( 1 − π ) } = ( ∑ j = 1 N μ j ( i + 1 ) ) − N π π ( 1 − π ) ∂ Q ( θ , θ ( i ) ) ∂ p = ∑ j = 1 N { μ j ( i + 1 ) π p Y j − 1 ( 1 − p ) − Y j ( Y j − p ) π p Y j ( 1 − p ) 1 − Y j } = ∑ j = 1 N { μ j ( i + 1 ) ( Y j − p ) p ( 1 − p ) } = ( ∑ j = 1 N μ j ( i + 1 ) Y j ) − ( p ∑ j = 1 N μ j ( i + 1 ) ) p ( 1 − p ) \begin{aligned} \frac{\partial Q(\theta,\theta^{(i)})}{\partial \pi} &= \sum_{j=1}^N\{\mu_j^{(i+1)}\frac{p^{Y_j}(1-p)^{1-Y_j} }{\pi p^{Y_j}(1-p)^{1-Y_j}} + (1-\mu_j^{(i+1)}) \frac{-q^{Y_j}(1-q)^{1-Y_j}}{(1-\pi)q^{Y_j}(1-q)^{1-Y_j}}\} \\ &= \sum_{j=1}^N\{\frac{\mu_j^{(i+1)}-\pi}{\pi(1-\pi)}\} \\ &=\frac{(\sum_{j=1}^N\mu_j^{(i+1)})-N\pi}{\pi(1-\pi)}\\ \frac{\partial Q(\theta,\theta^{(i)})}{\partial p} &= \sum_{j=1}^N\{\mu_j^{(i+1)}\frac{\pi p^{Y_j-1}(1-p)^{-Y_j}(Y_j-p) }{\pi p^{Y_j}(1-p)^{1-Y_j}}\} \\ &= \sum_{j=1}^N\{\frac{\mu_j^{(i+1)}(Y_j-p)}{p(1-p)}\} \\ &=\frac{(\sum_{j=1}^N\mu_j^{(i+1)}Y_j)-(p\sum_{j=1}^N\mu_j^{(i+1)})}{p(1-p)}\\ \end{aligned} ∂π∂Q(θ,θ(i))∂p∂Q(θ,θ(i))=j=1∑N{μj(i+1)πpYj(1−p)1−YjpYj(1−p)1−Yj+(1−μj(i+1))(1−π)qYj(1−q)1−Yj−qYj(1−q)1−Yj}=j=1∑N{π(1−π)μj(i+1)−π}=π(1−π)(∑j=1Nμj(i+1))−Nπ=j=1∑N{μj(i+1)πpYj(1−p)1−YjπpYj−1(1−p)−Yj(Yj−p)}=j=1∑N{p(1−p)μj(i+1)(Yj−p)}=p(1−p)(∑j=1Nμj(i+1)Yj)−(p∑j=1Nμj(i+1))
-
令偏導為 0 解 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)
∂ Q ( θ , θ ( i ) ) ∂ π = 令 0 ⇒ π ( i + 1 ) = 1 n ∑ i = 1 n μ j ( i + 1 ) ∂ Q ( θ , θ ( i ) ) ∂ p = 令 0 ⇒ { p ( i + 1 ) = ∑ i = 1 n μ j ( i + 1 ) Y j ∑ i = 1 n μ j ( i + 1 ) q ( i + 1 ) = ∑ i = 1 n ( 1 − μ j ( i + 1 ) ) Y j ∑ i = 1 n ( 1 − μ j ( i + 1 ) ) （ 對 稱 性 ） \begin{aligned} &\frac{\partial Q(\theta,\theta^{(i)})}{\partial \pi} \stackrel{令}{=}0 \Rightarrow \pi^{(i+1)} = \frac{1}{n} \sum_{i=1}^n \mu_j^{(i+1)} \\ &\frac{\partial Q(\theta,\theta^{(i)})}{\partial p}\stackrel{令}{=}0 \Rightarrow \left\{ \begin{aligned} &p^{(i+1)} = \frac{\sum_{i=1}^n \mu_j^{(i+1)}Y_j}{\sum_{i=1}^n \mu_j^{(i+1)}} \\ &q^{(i+1)} = \frac{\sum_{i=1}^n (1-\mu_j^{(i+1)})Y_j}{\sum_{i=1}^n (1-\mu_j^{(i+1)})} （對稱性） \end{aligned} \right. \end{aligned} ∂π∂Q(θ,θ(i))=令0⇒π(i+1)=n1i=1∑nμj(i+1)∂p∂Q(θ,θ(i))=令0⇒⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧p(i+1)=∑i=1nμj(i+1)∑i=1nμj(i+1)Yjq(i+1)=∑i=1n(1−μj(i+1))∑i=1n(1−μj(i+1))Yj（對稱性）
3.3 對比
-
對比 EM 和直接最大似然估計的優化目标
Q ( θ , θ ( i ) ) = ∑ j = 1 N { μ j ( i + 1 ) l o g [ π p Y j ( 1 − p ) 1 − Y j ] + ( 1 − μ j ( i + 1 ) ) l o g [ ( 1 − π ) q Y j ( 1 − q ) 1 − Y j ] } l o g P ( Y ∣ θ ) = ∑ j = 1 N l o g [ π p Y j ( 1 − p ) 1 − Y j + ( 1 − π ) q Y j ( 1 − q ) 1 − Y j ] \begin{aligned} &Q(\theta,\theta^{(i)}) = \sum_{j=1}^N\{\mu_j^{(i+1)}log[\pi p^{Y_j} (1-p)^{1-Y_j}] + (1-\mu_j^{(i+1)})log[(1-\pi) q^{Y_j} (1-q)^{1-Y_j}]\} \\ &logP(Y|\theta) = \sum_{j=1}^Nlog[\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}]\\ \end{aligned} Q(θ,θ(i))=j=1∑N{μj(i+1)log[πpYj(1−p)1−Yj]+(1−μj(i+1))log[(1−π)qYj(1−q)1−Yj]}logP(Y∣θ)=j=1∑Nlog[πpYj(1−p)1−Yj+(1−π)qYj(1−q)1−Yj]
-
對比 EM 和直接最大似然估計的偏導數方程組（因為最後要令等于0，重點對比分子即可）
{ ∂ Q ( θ , θ ( i ) ) ∂ π = ∑ j = 1 N { μ j ( i + 1 ) p Y j ( 1 − p ) 1 − Y j π p Y j ( 1 − p ) 1 − Y j + ( 1 − μ j ( i + 1 ) ) − q Y j ( 1 − q ) 1 − Y j ( 1 − π ) q Y j ( 1 − q ) 1 − Y j } ∂ Q ( θ , θ ( i ) ) ∂ p = ∑ j = 1 N { μ j ( i + 1 ) π p Y j − 1 ( 1 − p ) − Y j ( Y j − p ) π p Y j ( 1 − p ) 1 − Y j } { ∂ l o g P ( Y ∣ θ ) ∂ π = ∑ j = 1 n p Y j ( 1 − p ) 1 − Y j − q Y j ( 1 − q ) 1 − Y j π p Y j ( 1 − p ) 1 − Y j + ( 1 − π ) q Y j ( 1 − q ) 1 − Y j ∂ l o g P ( Y ∣ θ ) ∂ p = ∑ j = 1 n π p Y j − 1 ( 1 − p ) − Y j ( Y j − p ) π p Y j ( 1 − p ) 1 − Y j + ( 1 − π ) q Y j ( 1 − q ) 1 − Y j \begin{aligned} &\left\{ \begin{aligned} &\frac{\partial Q(\theta,\theta^{(i)})}{\partial \pi} = \sum_{j=1}^N\{\mu_j^{(i+1)}\frac{p^{Y_j}(1-p)^{1-Y_j} }{\pi p^{Y_j}(1-p)^{1-Y_j}} + (1-\mu_j^{(i+1)}) \frac{-q^{Y_j}(1-q)^{1-Y_j}}{(1-\pi)q^{Y_j}(1-q)^{1-Y_j}}\} \\ &\frac{\partial Q(\theta,\theta^{(i)})}{\partial p} = \sum_{j=1}^N\{\mu_j^{(i+1)}\frac{\pi p^{Y_j-1}(1-p)^{-Y_j}(Y_j-p) }{\pi p^{Y_j}(1-p)^{1-Y_j}}\} \end{aligned} \right.\\ &\left\{ \begin{aligned} \frac{\partial logP(Y|\theta)}{\partial \pi} & = \sum_{j=1}^n \frac{p^{Y_j}(1-p)^{1-Y_j}-q^{Y_j}(1-q)^{1-Y_j}}{\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}}\\ \frac{\partial logP(Y|\theta)}{\partial p} & = \sum_{j=1}^n \frac{\pi p^{Y_j-1}(1-p)^{-Y_j}(Y_j-p)}{\pi p^{Y_j}(1-p)^{1-Y_j} + (1-\pi)q^{Y_j}(1-q)^{1-Y_j}} \end{aligned} \right. \end{aligned} ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧∂π∂Q(θ,θ(i))=j=1∑N{μj(i+1)πpYj(1−p)1−YjpYj(1−p)1−Yj+(1−μj(i+1))(1−π)qYj(1−q)1−Yj−qYj(1−q)1−Yj}∂p∂Q(θ,θ(i))=j=1∑N{μj(i+1)πpYj(1−p)1−YjπpYj−1(1−p)−Yj(Yj−p)}⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧∂π∂logP(Y∣θ)∂p∂logP(Y∣θ)=j=1∑nπpYj(1−p)1−Yj+(1−π)qYj(1−q)1−YjpYj(1−p)1−Yj−qYj(1−q)1−Yj=j=1∑nπpYj(1−p)1−Yj+(1−π)qYj(1−q)1−YjπpYj−1(1−p)−Yj(Yj−p)
3.4 代碼
- 最後給出三硬币問題的代碼。網上搜到的基本都是同一份代碼的複制粘貼，其實那份代碼是适用于 jupyter notebook 那種互動式環境的，用了比較少見的
yield
- 我按照 1.1.2 節給出的流程重新寫了一遍，如下
import numpy as np import math class EM: def __init__(self, theta, y): self.theta = theta self.y = y self.mu = [0]*len(y) print('init pi={:.3f}; p={:.3f}; q={:.3f}'.format(theta[0],theta[1],theta[2])) # E_step def EStep(self): for j in range(len(self.y)): yj = self.y[j] pi,p,q = self.theta P_B_yj = pi * math.pow(p, yj) * math.pow((1-p), 1-yj) P_C_yj = (1-pi) * math.pow(q, yj) * math.pow((1-q), 1-yj) self.mu[j] = P_B_yj/(P_B_yj+P_C_yj) # m_step def MStep(self): sum_mu,sum_muy,sum_y = 0,0,0 N = len(self.y) for j in range(N): sum_mu += self.mu[j] sum_muy += self.mu[j]*self.y[j] sum_y += self.y[j] self.theta[0] = sum_mu/N # pi self.theta[1] = sum_muy/sum_mu # p self.theta[2] = (sum_y-sum_muy)/(N-sum_mu) # q def em(self): num = 1 while True: delta = self.theta.copy() self.EStep() self.MStep() print('ite{} pi={:.3f}; p={:.3f}; q={:.3f}'.format(num,self.theta[0],self.theta[1],self.theta[2])) num += 1 # 退出條件：theta變化量二範數小于門檻值 for i in range(3): delta[i] -= self.theta[i] if np.linalg.norm(delta) < 0.005: print('') break data=[1,1,0,1,0,0,1,0,1,1] test1 = EM([0.5, 0.5, 0.5],data) test1.em() test2 = EM([0.4, 0.6, 0.7],data) test2.em() ''' init pi=0.500; p=0.500; q=0.500 ite1 pi=0.500; p=0.600; q=0.600 ite2 pi=0.500; p=0.600; q=0.600 init pi=0.400; p=0.600; q=0.700 ite1 pi=0.406; p=0.537; q=0.643 ite2 pi=0.406; p=0.537; q=0.643 '''