[貝葉斯九]之EM演算法

03-03

<個人網頁blog已經上線，一大波乾貨即將來襲：https://faiculty.com/>

一、簡單介紹

EM(Expectaion Maximization)演算法(又稱為期望最大化方法)是一種迭代演算法，Dempster等人在1977年總結提出來的。簡單來說EM演算法就是一種含有隱變數的概率模型參數的極大似然估計。EM演算法的每次迭代由兩步組成：第一是求期望，第二是求極大。EM演算法在機器學習中有極為廣泛的應用。如常被用來學習高斯混合模型(Gaussian mixture model, 簡稱GMM)的參數。

那麼什麼是含有 $color{red}{隱變數的概率模型？}$ 這裡舉一個常用的三硬幣例子，假設我們有三枚硬幣：A、B和C，他們的質地都是不均勻的，假設他們正面朝上的概率分別是：a、b和c。現在弄一個拋硬幣的規則，先拋A硬幣，如果A正面朝上，那麼就拋B硬幣，否則就拋C硬幣。最後記下最終結果，正面朝上記為1，否則記為0。現在進行10次該實驗，假如得到的結果如下： 1，0，0，1，1，1，0，1，0，0。這個時候我們其實只得到了最終的結果，並不知道是B還是C硬幣的結果，因為不知道每次A硬幣的結果。這個時候A硬幣的拋擲就可以認為是一個隱含變數。但是問題是如何根據這個結果來估計這三個參數呢？

二、理論推導

2.1 演算法思想

在解決例子問題之前，我們先進行一些所謂枯燥的數學化定義，這樣或許能幫助理解和記憶。首先，假設 $Y$ 是最終觀測到的變數集(上述硬幣中的最終正反面結果)， $Z$ 是隱變數集(A硬幣結果)， $Theta(a,b,c)$ 是我們待求的參數集。根據之前我們對於極大似然估計的解釋，假設拋開隱變數集不管，我們最終的目的就是根據最終觀測到的變數集採用極大似然估計的方法來求解出參數集。所以，我們的目標函數就是最大化似然（這裡取似然函數的對數）。

$L(Theta|Y,Z) = log P(Y,Z|Theta)$

如果沒有 $Z$ 變數，如上所說，直接可以用極大似然估計的方法來估計參數。但是這裡多了一個隱變數 $Z$ ，所以EM演算法的精髓思想出來了：

初始化一個 $Theta^0$ 參數
E步：根據 $Theta^t$ 我們可以計算出 $Z$ 的期望值，記為 $Z^t$
M步：根據 $Z^t$ 我們可以利用極大似然法估計出參數 $Theta$ ，記為 $Theta^{t+1}$
重複上述EM步直到收斂

簡單闡述就是：其實這裡有兩類變數，一類是隱變數，一類是待求的參數變數。那麼普通的思路該怎麼求這個參數變數呢？由上述闡述可以知道，如果我們事先知道了隱變數就能利用極大似然來估計參數，如果我們知道了參數，那麼我們可以計算出隱變數集的期望。這裡就形成了一個制約，只要我們給出隱變數的初始值就能通過迭代達到兩類變數之間的平衡，也就是收斂。類似於我們在生活中的稱重，如果要將一類物品分為兩部分(比如糖果)，在沒有稱的情況下，往往我們在左右手進行掂量(這就有點像兩類變數)，如果左手上重了就分點到右手上，否則，從右手上扒拉點分到左手，直到感覺兩隻手上重量差不多。

所以這裡就落下了兩個最主要的問題：

為什麼是求 $Z$ 的期望？
一定會收斂？收斂到哪裡？

2.2 演算法推導

對於第一個問題，為什麼是求 $Z$ 的期望？其實上述只是簡單的假設，根據我們的目標需求，需要極大似然估計參數，所以需要知道隱變數，然後採取了這種假設、迭代的思路。實際上EM演算法並不是直接求解 $Z$ 的期望，而是將E、M兩步聯合。

假設迭代進行到第 $t$ 步，得到了參數 $Theta^{t}$ ，這個時候 $Z$ 的概率分布為： $P(Z|Y,Theta^t)$ 。我們不直接計算這個概率分布的期望值，而是計算對數似然函數 $L$ 關於 $Z$ 的期望值。同樣分為如下兩步。

E步(Expectation):計算對數似然函數 $L$ 關於 $Z$ 的期望值（一個條件概率的期望值，已知觀測結果 $Y$ 和參數 $Theta$ 的條件下，對數似然函數 $L$ 的期望值）。這個函數稱為 $Q$ 函數。

$egin{align} Q(Theta|Theta^t) &= E_z [log P(Y,Z|Theta)|Y,Theta^t] \ &= sum_z{P(Z|Y,Theta^t) log P(Y,Z|Theta)} ag 1 end{align}$

M步(Maximization):最大化期望似然函數。 $Theta^{t+1} = underset{Theta}{argmax} Q(Theta|Theta^t)$

重複上述兩個步驟，直到收斂到局部最優解。詳細的有關收斂問題，可以參考李航老師的《統計學習方法》。

這裡就通過簡單的方法介紹一下EM演算法的核心思路，然後主要是通過以下幾個例子來感受一下EM演算法。

三、三枚硬幣

這裡我們首先借用一下李航老師在《統計學習方法》中的三枚硬幣模型的例子進行闡述。

Step1: 先定義一些數學變數

$Y = {y_1, y_2, y_3, ........ , y_n}$ 代表最終拋擲硬幣的結果，每次的拋擲結果 $y in {0,1}$
$Z = {z_1, z_2, z_3, ....., z_n}$ 表示隱含變數，其實就是A硬幣的拋擲結果，所以 $z in {0, 1}$
$Theta$ 是模型參數，也就是 $Theta = pi, p, q$

Step2: 定義整個實驗的模型(似然函數)

在進行數學計算之前，模型定義是首要的。上述三硬幣的拋擲實驗可以用一個概率模型來定義。假設我們只拋擲一次，所以用小寫字母表示， $x$ 表示結果和 $z$ 隱變數。

$P(y|Theta) = P(y,z|Theta)$

因為 $z$ 只能取兩個值，由條件概率分布，我們展開來看。

$egin{align} P(y|Theta) &= P(y,z|Theta)\ &=underset{z}{sum}p(z|Theta)p(y|z,Theta)\ end{align}$

我們拆開來看，其實兩項都是二分布（簡單的拋擲硬幣過程, 假設A為正面朝上，然後進行一次拋硬幣，A反面朝上同樣）。所以我們可以繼續寫。

$egin{align} P(y|Theta) &= P(y,z|Theta)\ &=underset{z}{sum}p(z|Theta)p(y|z,Theta)\ &= pi p^y (1-p)^{1-y} + (1-pi) q^y (1-q)^{1-y} end{align}$

我們可以簡單的來核對一下這個概率模型寫得對不對。我們畫出一次拋擲硬幣的整個過程，並計算出相應的概率。然後帶入到上面的公式中就可以知道模型構建是否正確了。

$egin{array}{c|llll} case & ext{A(z)} & ext{B} & ext{C} & ext{概率}\ hline 1 & 1 & 1 & - & pi*p\ 2 & 1 & 0 & - & pi*(1-p) \ 3 & 0 & - & 1 & (1-pi)*q\ 4 & 0 & - & 0 & (1-pi)*(1-q)\ end{array}$

重複多次的原理也是如此，只不過因為進行的多次獨立實驗，所以計算概率直接用連乘累積。多次獨立實驗概率模型如下。由此我們首先得到了似然函數 $P(Y,Z|Theta)$ 。

$egin{align} P(Y|Theta) &= P(Y,Z|Theta)\ &=prod_{i=1}^nunderset{z_i}{sum}p(z_i|Theta)p(y_i|z_i,Theta)\ &= prod_{i=1}^npi p^{y_i} (1-p)^{1-y_i} + (1-pi) q^{y_i} (1-q)^{1-y_i} ag 2 end{align}$

這裡我們的目標就是求得參數 $Theta$ ，但是這個似然函數（聯乘）求導將非常複雜，所以在極大似然估計中一般都轉換為對數似然函數。但是依然非常複雜，所以用EM演算法迭代的思想來求解。

Step3：隱變數分布 $P(Z|Y,Theta^t)$

為了求得 $Q$ 函數，下一步我們的目標是得到隱變數的概率分布函數 $P(Z|Y,Theta^t)$ ，大家看到這個函數的時候千萬不要害怕，其實要從實際實驗的角度出發來理解這個公式。這個公式是隱變數 $Z$ 的分布函數，隱變數在實驗中其實就是一個0、1值， $z in {0,1}$ 。假設我們考慮只進行了一次拋硬幣的實驗，那麼當 $z$ 為1，就代表最後觀測結果來自B硬幣，否則來自C硬幣。所以我們可以拆開來寫這個函數。下面假設觀測結果來自B硬幣，即：A硬幣為1， $z=1$ 。後面緊隨一個拋B硬幣的二分布。

$egin{align} mu^{t+1} = P(z=1|y,Theta^t) = frac {pi p^y(1-p)^{1-y}}{pi p^y (1-p)^{1-y} + (1-pi) q^y (1-q)^{1-y}} end{align}$

我們將這個定義為 $mu^{t+1}$ ，按照之前EM演算法的迭代方法可知， $mu^{t+1}$ 是根據參數 $Theta^t$ 計算得到。

同理，如果 $z$ 為0也是一樣的。那麼 $P(Z|Y,Theta^t)$ 這個函數就呼之欲出了。這裡我們只需要考慮一次實驗， $z_i$ .

Step4: E步，得到 $Q$ 函數由上述理論推導部分可以得到。然後將上述的似然函數 $(2)$ 和隱變數的分布函數 $(3)$ 帶入 $Q$ 函數 $(1)$ 中就行。

$egin{align} Q(Theta|Theta^t) &= E_z [log P(Y,Z|Theta)|Y,Theta^t] \ &= sum_z{P(Z|Y,Theta^t) log P(Y,Z|Theta)} \ &= sum_{j=1}^{N} sum_z {P(z_j|y_j,Theta^t) log P(y_j,z_j|Theta^t)} \ &= sum_{j=1}^{N} mu_j^{t+1} log (pi p^y (1-p)^{1-y}) + (1-mu_j^{t+1}) log ((1-pi) q^y (1-q)^{1-y}) ag 4 end{align}$

Step5: M步,對 $Q$ 函數求導，極大似然估計

這裡根據Step4我們得到了隱變數的期望（也就是得到了隱變數），由此可以直接對上述的(4)式進行求導，得到相應的參數。

$egin{align} frac{partial Q}{partial pi} &= (frac{mu_1^{t+1}}{pi} - frac{1-mu_1^{t+1}}{1-pi}) + ...... +(frac{mu_N^{t+1}}{pi} - frac{1-mu_N^{t+1}}{1-pi}) \ &= frac{mu_1^{t+1}-pi}{pi(1-pi)} + ...... + frac{mu_N^{t+1} - pi}{pi(1-pi)}\ &=0 end{align}$

由此得到 $pi^{t+1} = frac{1}{N} sum_{j=1}^{N} mu_j^{t+1}$

同理可得：

$p^{t+1} = frac{sum_{j=1}^{N} mu_j^{t+1} y_j}{sum_{j=1}^{N} mu_j^{t+1}}$

$q^{t+1} = frac{sum_{j=1}^{N} (1-mu_j^{t+1}) y_j}{sum_{j=1}^{N} (1-mu_j^{t+1})}$

Step6: 進行數值計算并迭代

如何判斷停止呢？一般是給出一個較小的整數 $varepsilon_1$ ， $varepsilon_2$ ，若滿足如下條件則停止迭代。

$egin{align} ||Theta^{t+1} - Theta^{t}|| < varepsilon_1 \ or qquad ||Q(Theta^{t+1}|Theta^t) - Q(Theta^t | Theta^t) || < varepsilon_2 end{align}$

四、一個簡單例子

下面我們將編程實現一個最簡單的正態分布參數估計的例子，感受一下EM演算法。

case1:

如下圖所示是兩組一維正態分布數據。這兩組是 $color{red}{不帶隱含變數的正態分布}$ ，我們可以很明顯的看出這兩組正態分布數據，紅色的點是一類，藍色的是另外一類。也可以大概估計出這個例子中兩組正態分布的均值，比如紅色類別大概是3左右。這就是極大似然估計所處理的一個場景： $color{blue}{不帶隱變數的參數估計方法.}$

case2:如下圖所示也告訴是兩組一維正態分布的數據。但是是帶有隱變數的，所以我們完全看不出兩個類別，也不太好利用極大似然估計的方法來找到兩組參數。這個時候就是EM演算法的角斗場， $color{blue}{利用反覆迭代的思路。}$ 為了方便理解，我們完全可以理解為用一個濾波器一樣的東西在給出的數據上滑動，看哪個濾波器（一組相應的正態分布的參數）能產生最小的誤差。

下面我們將按照上述三枚硬幣的例子寫一小段代碼來感受EM演算法迭代求解的過程。

Step0：隨機數據生成首先，顯然先生成這樣兩組正態分布的數據。

import numpy as npfrom scipy import statsnp.random.seed(110) # for reproducible random results# set parametersred_mean = 3red_std = 0.8blue_mean = 7blue_std = 2# draw 20 samples from normal distributions with red/blue parametersred = np.random.normal(red_mean, red_std, size=20)blue = np.random.normal(blue_mean, blue_std, size=20)# sort all the data for later useboth_colours = np.sort(np.concatenate((red, blue))) # for later use...

我們可以輸出結果看看，方便之後與估計的參數進行對比。

>>> np.mean(red)2.802>>> np.std(red)0.871>>> np.mean(blue)6.932>>> np.std(blue)2.195

Step2: 初始化參數

這裡我們的 $Theta$ 參數就是兩組正態分布的參數： $均值mu，方差：std$

# estimates for the meanred_mean_guess = 1.1blue_mean_guess = 9# estimates for the standard deviationred_std_guess = 2blue_std_guess = 1.7

初始化的參數繪製的圖像如下所示，看見這張圖是不是有點似曾相識，在貝葉斯決策中我們就畫過這樣的圖，判斷是否屬於某個類別的時候，分界面是一個點，比如下圖中，紅色和藍色正態分布圖的交叉點就是分界點，小於這個點就是屬於紅色，否則是藍色。但是EM演算法更加強，不僅能找出這個分界面，而且能估計出參數。

Step2: 計算似然函數

likelihood_of_red = stats.norm(red_mean_guess, red_std_guess).pdf(both_colours)likelihood_of_blue = stats.norm(blue_mean_guess, blue_std_guess).pdf(both_colours)likelihood_total = likelihood_of_red + likelihood_of_bluered_weight = likelihood_of_red / likelihood_totalblue_weight = likelihood_of_blue / likelihood_total

Step3：估計參數

def estimate_mean(data, weight): return np.sum(data * weight) / np.sum(weight)def estimate_std(data, weight, mean): variance = np.sum(weight * (data - mean)**2) / np.sum(weight) return np.sqrt(variance)# new estimates for standard deviationblue_std_guess = estimate_std(both_colours, blue_weight, blue_mean_guess)red_std_guess = estimate_std(both_colours, red_weight, red_mean_guess)# new estimates for meanred_mean_guess = estimate_mean(both_colours, red_weight)blue_mean_guess = estimate_mean(both_colours, blue_weight)

**Step4：迭代20次後結果 **

下圖是迭代的過程繪製在圖像上，可以看出擬合程度。

最終的結果圖如下所示。分界點大概在4.2左右，參考最開始給出的那個紅藍分解點可以看出，估計還是比較準確的。

最下面是估計結果參數對比實際的參數的表格。

| EM guess | Actual | Delta----------+----------+--------+-------Red mean | 2.910 | 2.802 | 0.108Red std | 0.854 | 0.871 | -0.017Blue mean | 6.838 | 6.932 | -0.094Blue std | 2.227 | 2.195 | 0.032

參考文獻

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

[2] 周志華，《機器學習》

[3]機器學習演算法系列之一】EM演算法實例分析

[4]從最大似然到EM演算法淺解

[5]CS229 Lecture notes by Andrew Ng

[6]What is an intuitive explanation of the Expectation Maximization technique?