人人都懂EM演算法

05-20

人人都懂EM演算法

來自專欄人工智慧隨筆

估計有很多入門機器學習的同學在看到EM演算法的時候會有種種疑惑：EM演算法到底是個什麼玩意？它能做什麼？它的應用場景是什麼？網上的公式推導怎麼看不懂？

下面我從極大似然估計開始，過渡到EM演算法，講解EM演算法最核心的idea，以及EM演算法的具體步驟。鑒於網上很多博客文章都是直接翻譯吳恩達的課程筆記內容，有很多推導步驟都是跳躍性的，我會把這些中間步驟彌補上，讓大家都能看懂EM演算法的推導過程。最後以一個二硬幣模型作為EM演算法的一個實例收尾。希望閱讀本篇文章之後能對EM演算法有更深的了解和認識。

極大似然和EM(Expectation Maximization)演算法，與其說是一種演算法，不如說是一種解決問題的思想，解決一類問題的框架，和線性回歸，邏輯回歸，決策樹等一些具體的演算法不同，極大似然和EM演算法是更加抽象，是很多具體演算法的基礎。

1. 從極大似然到EM

1.1 極大似然

1.1.1 問題描述

假設我們需要調查我們學校學生的身高分布。我們先假設學校所有學生的身高服從正態分布 $N(mu,sigma^2)$ 。(注意：極大似然估計的前提一定是要假設數據總體的分布，如果不知道數據分布，是無法使用極大似然估計的)，這個分布的均值 $mu$ 和方差 $sigma^2$ 未知，如果我們估計出這兩個參數，那我們就得到了最終的結果。那麼怎樣估計這兩個參數呢？

我們可以先對學生進行抽樣。假設我們隨機抽到了 200 個人（也就是 200 個身高的樣本數據，為了方便表示，下面，「人」的意思就是對應的身高）。然後統計抽樣這 200 個人的身高。根據這 200 個人的身高估計均值 $mu$ 和方差 $sigma^2$ 。

用數學的語言來說就是：為了統計學校學生的身高分布，我們獨立地按照概率密度 $p(x|θ)$ 抽取了 200 個（身高），組成樣本集 $X={x_1,x_2,…,x_N}$ (其中 $x_i$ 表示抽到的第 $i$ 個人的身高，這裡 N 就是 200，表示樣本個數)，我們想通過樣本集 X 來估計出未知參數 $θ$ 。這裡概率密度 $p(x|θ)$ 服從高斯分布 $N(mu,sigma^2)$ ，其中的未知參數是 $θ=[mu, sigma]^T$ 。

那麼問題來了怎樣估算參數 $heta$ 呢？

1.1.2 估算參數

我們先回答幾個小問題：

問題一：抽到這 200 個人的概率是多少呢？

$L( heta) = L(x_1, x_2, cdots , x_n; heta) = prod limits _{i=1}^{n}p(x_i| heta), quad heta in Theta$

n為抽取的樣本的個數，本例中 $n=200$ ，這個概率反映了，在概率密度函數的參數是 $θ$ 時，得到 X 這組樣本的概率。因為這裡 L 是已知的，也就是說我抽取到的這 200 個人的身高可以測出來。而 $θ$ 是未知了，則上面這個公式只有 $θ$ 是未知數，所以 L 是 $θ$ 的函數。

這個函數反映的是在不同的參數 $θ$ 取值下，取得當前這個樣本集的可能性，因此稱為參數 $θ$ 相對於樣本集 X 的似然函數（likelihood function），記為 $θ$ 。

為了便於分析，對 L 取對數，將其變成連加的，稱為對數似然函數，如下式：

$H( heta) = ext{ln} L( heta) = ext{ln} prod limits _{i=1}^{n}p(x_i; heta) = sum limits _{i=1}^{n} ext{ln} p(x_i; heta)$

問題二：學校那麼多學生，為什麼就恰好抽到了這 200 個人 ( 身高) 呢？

在學校那麼學生中，我一抽就抽到這 200 個學生（身高），而不是其他人，那是不是表示在整個學校中，這 200 個人（的身高）出現的概率極大啊，也就是其對應的似然函數 $L(θ)$ 極大，即

$hat heta = ext{argmax} L( heta)$

$hat heta$ 這個叫做 $θ$ 的極大似然估計量，即為我們所求的值。

問題三：那麼怎麼極大似然函數？

求 $L( heta)$ 對所有參數的偏導數，然後讓這些偏導數為 0，假設有 $n$ 個參數，就有 $n$ 個方程組成的方程組，那麼方程組的解就是似然函數的極值點了，從而得到對應的 $θ$ 了。

1.1.3 極大似然估計總結

極大似然估計你可以把它看作是一個反推。多數情況下我們是根據已知條件來推算結果，而極大似然估計是已經知道了結果，然後尋求使該結果出現的可能性極大的條件，以此作為估計值。

比如說，

假如一個學校的學生男女比例為 9:1 (條件)，那麼你可以推出，你在這個學校里更大可能性遇到的是男生 (結果)；
假如你不知道那女比例，你走在路上，碰到100個人，發現男生就有90個 (結果)，這時候你可以推斷這個學校的男女比例更有可能為 9:1 (條件)，這就是極大似然估計。

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

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

1.1.4 求極大似然函數估計值的一般步驟：

（1）寫出似然函數；

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

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

（4）解似然方程，得到的參數。

1.1.5 極大似然函數的應用

應用一：回歸問題中的極小化平方和 （極小化代價函數）

假設線性回歸模型具有如下形式: $h(x) = sum limits _{i=1}^{d} heta_jx_j + epsilon = heta^Tx + epsilon$ ，其中 $x in R^{1 imes d}, heta in R^{1 imes d}，誤差 epsilon in R$ ，誤差 $X = (x_1, cdots, x_m)^T in R^{m imes d}, y in R^{m imes 1}$ ，如何求 $θ$ 呢？

最小二乘估計：最合理的參數估計量應該使得模型能最好地擬合樣本數據，也就是估計值和觀測值之差的平方和最小，其推導過程如下所示：

$J( heta)= sum limits _{i=1}^n (h_{ heta}(x_i)- y_i) ^2$

求解方法是通過梯度下降演算法，訓練數據不斷迭代得到最終的值。

極大似然法：最合理的參數估計量應該使得從模型中抽取 m 組樣本觀測值的概率極大，也就是似然函數極大。

假設誤差項 $epsilon in N(0, sigma^2)$ ，則 $y_i in N( heta x_i, sigma^2)$ (建議複習一下正態分布的概率密度函數和相關的性質)

$p(y_i|x_i; heta) = frac{1}{sqrt{2 pi}sigma}exp(-frac{(y_i- heta^Tx_{i})^2}{2sigma^2}) \ egin {align*}L( heta) &= prod limits _{i=1}^mp(y_i|x_i; heta) \ &= prod limits _{i=1}^mfrac{1}{sqrt{2 pi}sigma}exp(-frac{(y_i- heta^Tx_{i})^2}{2sigma^2})end{align*}$

$egin {align*}H( heta) &= log(L( heta)) \ &= ext{log} prod limits _{i=1}^mfrac{1}{sqrt{2 pi}sigma}exp(-frac{(y_i- heta^Tx_{i})^2}{2sigma^2}) \&= sum limits _{i=1}^m( ext{log} frac{1}{sqrt{2 pi}sigma}exp(-frac{(y_i- heta^Tx_{i})^2}{2sigma^2})) \ & = -frac{1}{2sigma^2} sum limits _{i=1}^m(y_i- heta^Tx_{i})^2 - m ext{ln} sigma sqrt{2pi} end{align*}$

令 $J( heta) = frac{1}{2} sum limits _{i=1}^m(y_i- heta^Tx_{i})^2$ 則 $arg max limits_{ heta} H( heta) Leftrightarrow arg min limits_{ heta} J( heta)$ ，即將極大似然函數等價於極小化平方和。

這時可以發現，此時的極大化似然函數和最初的最小二乘損失函數的估計結果是等價的。但是要注意這兩者只是恰好有著相同的表達結果，原理和出發點完全不同。

應用二：分類問題中極小化交叉熵 （極小化代價函數）

在分類問題中，交叉熵的本質就是似然函數的極大化，邏輯回歸的假設函數為：

$h(x) = hat y = frac 1 {1+e^{- heta^Tx + b}}$

根據之前學過的內容我們知道 $hat y = p(y=1|x, heta)$ ，

當 $y=1$ 時， $p_1 = p(y = 1|x, heta) = hat y$

當 $y=0$ 時， $p_0 = p(y = 0|x, heta) = 1- hat y$

合併上面兩式子，可以得到

$p(y|x， heta) = hat y^y(1- hat y)^{1- y}$ $egin {align*}L( heta) &= prod limits _{i=1}^mp(y_i|x_i; heta) \ &= prod limits _{i=1}^mhat y_i^{y_i}(1- hat y_i)^{1- y_i}end{align*}$

$egin {align*}H( heta) &= ext{log}(L( heta)) \ &= ext{log}prod limits _{i=1}^mhat y_i^{y_i}(1- hat y_i)^{1- y_i} \&= sum limits _{i=1}^m ext{log} hat y_i^{y_i}(1- hat y_i)^{1- y_i} \ & =sum limits _{i=1}^m y_i ext{log} hat y_i + (1-y_i) ext{log} (1 - hat y_i) end{align*}$

令 $J( heta) = -H( heta) = -sum limits _{i=1}^m y_i ext{log} hat y_i + (1-y_i) ext{log} (1 - hat y_i)$ 則 $arg max limits_{ heta} H( heta) Leftrightarrow arg min limits_{ heta} J( heta)$ ，即將極大似然函數等價於極小化平方和。

1.2 EM演算法

1.2.1 問題描述

上面我們先假設學校所有學生的身高服從正態分布 $N(mu,sigma^2)$ 。實際情況並不是這樣的，男生和女生分別服從兩種不同的正態分布，即男生 $in N(mu_1, sigma_1^2)$ ，女生 $in N(mu_2, sigma_2^2)$ ，(注意：EM演算法和極大似然估計的前提是一樣的，都要假設數據總體的分布，如果不知道數據分布，是無法使用EM演算法的)。那麼該怎樣評估學生的身高分布呢？

簡單啊，我們可以隨便抽 100 個男生和 100 個女生，將男生和女生分開，對他們單獨進行極大似然估計。分別求出男生和女生的分布。

假如某些男生和某些女生好上了，糾纏起來了。咱們也不想那麼殘忍，硬把他們拉扯開。這時候，你從這 200 個人（的身高）裡面隨便給我指一個人（的身高），我都無法確定這個人（的身高）是男生（的身高）還是女生（的身高）。用數學的語言就是，抽取得到的每個樣本都不知道是從哪個分布來的。那怎麼辦呢？

1.2.2 EM 演算法

這個時候，對於每一個樣本或者你抽取到的人，就有兩個問題需要估計了，一是這個人是男的還是女的，二是男生和女生對應的身高的正態分布的參數是多少。這兩個問題是相互依賴的：

當我們知道了每個人是男生還是女生，我們可以很容易利用極大似然對男女各自的身高的分布進行估計。
反過來，當我們知道了男女身高的分布參數我們才能知道每一個人更有可能是男生還是女生。例如我們已知男生的身高分布為 $N(mu_1 = 172, sigma^2_1=5^2)$ ，女生的身高分布為 $N(mu_2 = 162, sigma^2_2=5^2)$ ，一個學生的身高為180，我們可以推斷出這個學生為男生的可能性更大。

但是現在我們既不知道每個學生是男生還是女生，也不知道男生和女生的身高分布。這就成了一個先有雞還是先有蛋的問題了。雞說，沒有我，誰把你生出來的啊。蛋不服，說，沒有我，你從哪蹦出來啊。為了解決這個你依賴我，我依賴你的循環依賴問題，總得有一方要先打破僵局，不管了，我先隨便整一個值出來，看你怎麼變，然後我再根據你的變化調整我的變化，然後如此迭代著不斷互相推導，最終就會收斂到一個解。這就是EM演算法的基本思想了。

EM的意思是「Expectation Maximization」，具體方法為：

先設定男生和女生的身高分布參數(初始值)，例如男生的身高分布為 $N(mu_1 = 172, sigma^2_1=5^2)$ ，女生的身高分布為 $N(mu_2 = 162, sigma^2_2=5^2)$ ，當然了，剛開始肯定沒那麼准；
然後計算出每個人更可能屬於第一個還是第二個正態分布中的（例如，這個人的身高是180，那很明顯，他極大可能屬於男生），這個是屬於Expectation 一步；
我們已經大概地按上面的方法將這 200 個人分為男生和女生兩部分，我們就可以根據之前說的極大似然估計分別對男生和女生的身高分布參數進行估計（這不變成了極大似然估計了嗎？極大即為Maximization）這步稱為 Maximization；
然後，當我們更新這兩個分布的時候，每一個學生屬於女生還是男生的概率又變了，那麼我們就再需要調整E步；
……如此往複，直到參數基本不再發生變化或滿足結束條件為止。

1.2.3 總結

上面的學生屬於男生還是女生我們稱之為隱含參數，女生和男生的身高分布參數稱為模型參數。

EM 演算法解決這個的思路是使用啟發式的迭代方法，既然我們無法直接求出模型分布參數，那麼我們可以先猜想隱含參數（EM 演算法的 E 步），接著基於觀察數據和猜測的隱含參數一起來極大化對數似然，求解我們的模型參數（EM演算法的M步)。由於我們之前的隱含參數是猜測的，所以此時得到的模型參數一般還不是我們想要的結果。我們基於當前得到的模型參數，繼續猜測隱含參數（EM演算法的 E 步），然後繼續極大化對數似然，求解我們的模型參數（EM演算法的M步)。以此類推，不斷的迭代下去，直到模型分布參數基本無變化，演算法收斂，找到合適的模型參數。

一個最直觀了解 EM 演算法思路的是 K-Means 演算法。在 K-Means 聚類時，每個聚類簇的質心是隱含數據。我們會假設 K 個初始化質心，即 EM 演算法的 E 步；然後計算得到每個樣本最近的質心，並把樣本聚類到最近的這個質心，即 EM 演算法的 M 步。重複這個 E 步和 M 步，直到質心不再變化為止，這樣就完成了 K-Means 聚類。

2. EM演算法推導

2.1 基礎知識

2.1.1 凸函數

設是定義在實數域上的函數，如果對於任意的實數，都有：

$f ge0$

那麼是凸函數。若不是單個實數，而是由實數組成的向量，此時，如果函數的 Hesse 矩陣是半正定的，即

$H ge 0$

是凸函數。特別地，如果 $f > 0$ 或者 $H > 0$ ，稱為嚴格凸函數。

2.1.2 Jensen不等式

如下圖，如果函數 $f$ 是凸函數， $x$ 是隨機變數，有 0.5 的概率是 a，有 0.5 的概率是 b， $x$ 的期望值就是 a 和 b 的中值了那麼：

$E[f(x)] ge f(E(x))$

其中， $E[f(x)] = 0.5f(a) + 0.5 f(b)，f(E(x)) = f(0.5a + 0.5b)$ ，這裡 a 和 b 的權值為 0.5, $f(a)$ 與 a 的權值相等， $f(b)$ 與 b 的權值相等。

特別地，如果函數 $f$ 是嚴格凸函數，當且僅當： $p(x = E(x)) = 1$ (即隨機變數是常量) 時等號成立。

註：若函數 $f$ 是凹函數，Jensen不等式符號相反。

2.1.3 期望

對於離散型隨機變數 X 的概率分布為 $p_i = p{X=x_i}$ ，數學期望 $E(X)$ 為：

$E(X) = sum limits _i x_ip_i$

$p_i$ 是權值，滿足兩個條件 $p_i ge 0，sum limits _i p_i = 1$ 。

若連續型隨機變數X的概率密度函數為 $f(x)$ ，則數學期望 $E(X)$ 為：

$E(X) = int _ {-infty} ^{+infty} xf(x) dx$

設 $Y = g(X)$ ，若 $X$ 是離散型隨機變數，則：

$E(Y) = sum limits _i g(x_i)p_i$

若 $X$ 是連續型隨機變數，則：

$E(X) = int _ {-infty} ^{+infty} g(x)f(x) dx$

2.2 EM演算法的推導

對於 $m$ 個相互獨立的樣本 $x=(x^{(1)},x^{(2)},...x^{(m)})$ ，對應的隱含數據 $z=(z^{(1)},z^{(2)},...z^{(m)})$ ，此時 $(x,z)$ 即為完全數據，樣本的模型參數為 $θ$ , 則觀察數據 $x^{(i)}$ 的概率為 $P(x^{(i)}| heta)$ ，完全數據 $(x^{(i)},z^{(i)})$ 的似然函數為 $P(x^{(i)},z^{(i)}| heta)$ 。

假如沒有隱含變數 $z$ ，我們僅需要找到合適的 $heta$ 極大化對數似然函數即可：

$heta =arg max limits_{ heta}L( heta) = arg max limits_{ heta}sumlimits_{i=1}^m logP(x^{(i)}| heta)$

增加隱含變數 $z$ 之後，我們的目標變成了找到合適的 $heta$ 和 $z$ 讓對數似然函數極大：

$heta, z = arg max limits_{ heta,z}L( heta, z) = arg max limits_{ heta,z}sumlimits_{i=1}^m logsumlimits_{z^{(i)}}P(x^{(i)}， z^{(i)}| heta)$

如果對分別對未知的 $heta$ 和 $z$ 分別求偏導，由於 $logP(x^{(i)}| heta)$ 是 $P(x^{(i)}， z^{(i)}| heta)$ 邊緣概率(建議沒基礎的同學網上搜一下邊緣概率的概念)，轉化為 $logP(x^{(i)}| heta)$ 求導後形式會非常複雜（可以想像下 $log(f_1(x)+ f_2(x)+…$ 複合函數的求導) ，所以很難求解得到 $heta$ 和 $z$ 。那麼我們想一下可不可以將加號從 log 中提取出來呢？我們對這個式子進行縮放如下： $egin{align} sumlimits_{i=1}^m logsumlimits_{z^{(i)}}P(x^{(i)}， z^{(i)}| heta) & = sumlimits_{i=1}^m logsumlimits_{z^{(i)}}Q_i(z^{(i)})frac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})} ag{1} \ & geq sumlimits_{i=1}^m sumlimits_{z^{(i)}}Q_i(z^{(i)})logfrac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})} ag{2} end{align}$

上面第(1)式引入了一個未知的新的分布 $Q_i(z^{(i)})$ ，滿足：

$sum limits _z Q_i(z)=1,0 le Q_i(z)le 1$

第(2)式用到了 Jensen 不等式 (對數函數是凹函數)：

$log(E(y)) ge E(log(y))$

其中：

$E(y) = sumlimits_ilambda_iy_i, lambda_i geq 0, sumlimits_ilambda_i =1$

$y_i = frac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})}$

$lambda_i = Q_i(z^{(i)})$

也就是說 $frac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})}$ 為第 i 個樣本， $Q_i(z^{(i)})$ 為第 i 個樣本對應的權重，那麼：

$E(frac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})}) = sumlimits_{z^{(i)}}Q_i(z^{(i)})frac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})}$

上式我實際上是我們構建了 $L( heta, z)$ 的下界，我們發現實際上就是 $frac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})}$ 的加權平均，由於上面講過權值 $Q_i(z^{(i)})$ 累積和為1，因此上式是 $frac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})}$ 的期望，這就是Expectation的來歷啦。下一步要做的就是尋找一個合適的 $Q_i(z)$ 最優化這個下界(M步)。

假設 $heta$ 已經給定，那麼 $logL( heta)$ 的值就取決於 $Q_i(z)$ 和 $p(x^{(i)},z^{(i)})$ 了。我們可以通過調整這兩個概率使下界逼近 $logL( heta)$ 的真實值，當不等式變成等式時，說明我們調整後的下界能夠等價於 $logL( heta)$ 了。由 Jensen 不等式可知，等式成立的條件是隨機變數是常數，則有： $frac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})} =c$

其中 c 為常數，對於任意 $i$ ，我們得到：

${P(x^{(i)}， z^{(i)}| heta)} =c{Q_i(z^{(i)})}$

方程兩邊同時累加和：

$sumlimits_{z} {P(x^{(i)}， z^{(i)}| heta)} = csumlimits_{z} {Q_i(z^{(i)})}$

由於 $sumlimits_{z}Q_i(z^{(i)}) =1$ 。從上面兩式，我們可以得到：

$sumlimits_{z} {P(x^{(i)}， z^{(i)}| heta)} = c$

$Q_i(z^{(i)}) = frac{P(x^{(i)}， z^{(i)}| heta)}{c} = frac{P(x^{(i)}， z^{(i)}| heta)}{sumlimits_{z}P(x^{(i)}， z^{(i)}| heta)} = frac{P(x^{(i)}， z^{(i)}| heta)}{P(x^{(i)}| heta)} = P( z^{(i)}|x^{(i)}， heta)$

其中：

邊緣概率公式： $P(x^{(i)}| heta) = sumlimits_{z}P(x^{(i)}， z^{(i)}| heta)$

條件概率公式： $frac{P(x^{(i)}， z^{(i)}| heta)}{P(x^{(i)}| heta)} = P( z^{(i)}|x^{(i)}， heta)$

從上式可以發現 $Q(z)$ 是已知樣本和模型參數下的隱變數分布。

如果 $Q_i(z^{(i)}) = P( z^{(i)}|x^{(i)}， heta))$ , 則第 (2) 式是我們的包含隱藏數據的對數似然的一個下界。如果我們能極大化這個下界，則也在嘗試極大化我們的對數似然。即我們需要極大化下式： $arg max limits_{ heta} sumlimits_{i=1}^m sumlimits_{z^{(i)}}Q_i(z^{(i)})logfrac{P(x^{(i)}， z^{(i)}| heta)}{Q_i(z^{(i)})}$

至此，我們推出了在固定參數 $heta$ 後分布 $Q_i(z^{(i)})$ 的選擇問題，從而建立了 $logL( heta)$ 的下界，這是 E 步，接下來的M 步驟就是固定 $Q_i(z^{(i)})$ 後，調整 $heta$ ，去極大化 $logL( heta)$ 的下界。

去掉上式中常數的部分 $Q_i(z^{(i)})$ ，則我們需要極大化的對數似然下界為：

$arg max limits_{ heta} sumlimits_{i=1}^m sumlimits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}， z^{(i)}| heta)}$

2.3 EM演算法流程

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

輸入：觀察數據 $x=(x^{(1)},x^{(2)},...x^{(m)})$ ，聯合分布 $p(x,z | heta)$ ，條件分布 $p(z|x, heta)$ ，極大迭代次數 $J$ 。

1) 隨機初始化模型參數 $heta$ 的初值 $heta^0$

2) $ext{for j from 1 to J}$ ：

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

$Q_i(z^{(i)}) := P( z^{(i)}|x^{(i)}， heta))$

M步：極大化 $L( heta)$ ,得到 $heta$ :

$heta : = arg max limits_{ heta}sumlimits_{i=1}^msumlimits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}， z^{(i)}| heta)}$

重複E、M步驟直到 $heta$ 收斂

輸出：模型參數 $heta$

2.4 EM演算法另一種理解

坐標上升法（Coordinate ascent）(類似於梯度下降法，梯度下降法的目的是最小化代價函數，坐標上升法的目的是最大化似然函數；梯度下降每一個循環僅僅更新模型參數就可以了，EM演算法每一個循環既需要更新隱含參數和也需要更新模型參數，梯度下降和坐標上升的詳細分析參見攀登傳統機器學習的珠峰-SVM (下))：

圖中的直線式迭代優化的路徑，可以看到每一步都會向最優值前進一步，而且前進路線是平行於坐標軸的，因為每一步只優化一個變數。

這猶如在x-y坐標系中找一個曲線的極值，然而曲線函數不能直接求導，因此什麼梯度下降方法就不適用了。但固定一個變數後，另外一個可以通過求導得到，因此可以使用坐標上升法，一次固定一個變數，對另外的求極值，最後逐步逼近極值。對應到EM上，E步：固定 θ，優化Q；M步：固定 Q，優化 θ；交替將極值推向極大。

2.5 EM演算法的收斂性思考

EM演算法的流程並不複雜，但是還有兩個問題需要我們思考：

1） EM演算法能保證收斂嗎？

2） EM演算法如果收斂，那麼能保證收斂到全局極大值嗎？　

首先我們來看第一個問題, EM 演算法的收斂性。要證明 EM 演算法收斂，則我們需要證明我們的對數似然函數的值在迭代的過程中一直在增大。即：

$sumlimits_{i=1}^m logP(x^{(i)}| heta^{j+1}) geq sumlimits_{i=1}^m logP(x^{(i)}| heta^{j})$

由於：

$L( heta, heta^{j}) = sumlimits_{i=1}^msumlimits_{z^{(i)}}P( z^{(i)}|x^{(i)}， heta^{j}))log{P(x^{(i)}， z^{(i)}| heta)}$

令：

$H( heta, heta^{j}) = sumlimits_{i=1}^msumlimits_{z^{(i)}}P( z^{(i)}|x^{(i)}， heta^{j}))log{P( z^{(i)}|x^{(i)}， heta)}$

上兩式相減得到：

$sumlimits_{i=1}^m logP(x^{(i)}| heta) = L( heta, heta^{j}) - H( heta, heta^{j})$

在上式中分別取 $θ$ 為 $θ^j$ 和 $θ^{j+1}$ ，並相減得到：

$sumlimits_{i=1}^m logP(x^{(i)}| heta^{j+1}) - sumlimits_{i=1}^m logP(x^{(i)}| heta^{j}) = [L( heta^{j+1}, heta^{j}) - L( heta^{j}, heta^{j}) ] -[H( heta^{j+1}, heta^{j}) - H( heta^{j}, heta^{j}) ]$

要證明EM演算法的收斂性，我們只需要證明上式的右邊是非負的即可。

由於 $θ^{j+1}$ 使得 $L(θ,θ^j)$ 極大，因此有：

$L( heta^{j+1}, heta^{j}) - L( heta^{j}, heta^{j}) geq 0$

而對於第二部分，我們有：

$egin{align} H( heta^{j+1}, heta^{j}) - H( heta^{j}, heta^{j}) & = sumlimits_{i=1}^msumlimits_{z^{(i)}}P( z^{(i)}|x^{(i)}， heta^{j})logfrac{P( z^{(i)}|x^{(i)}， heta^{j+1})}{P( z^{(i)}|x^{(i)}， heta^j)} ag{3} \ & leq sumlimits_{i=1}^mlog(sumlimits_{z^{(i)}}P( z^{(i)}|x^{(i)}， heta^{j})frac{P( z^{(i)}|x^{(i)}， heta^{j+1})}{P( z^{(i)}|x^{(i)}， heta^j)}) ag{4} \ & = sumlimits_{i=1}^mlog(sumlimits_{z^{(i)}}P( z^{(i)}|x^{(i)}， heta^{j+1})) = 0 ag{5} end{align}$

其中第（4）式用到了Jensen不等式，只不過和第二節的使用相反而已，第（5）式用到了概率分布累積為1的性質。

至此，我們得到了： $sumlimits_{i=1}^m logP(x^{(i)}| heta^{j+1}) - sumlimits_{i=1}^m logP(x^{(i)}| heta^{j}) geq 0$ ，證明了EM演算法的收斂性。

從上面的推導可以看出，EM 演算法可以保證收斂到一個穩定點，但是卻不能保證收斂到全局的極大值點，因此它是局部最優的演算法，當然，如果我們的優化目標 $L(θ,θ^j)$ 是凸的，則EM演算法可以保證收斂到全局極大值，這點和梯度下降法這樣的迭代演算法相同。至此我們也回答了上面提到的第二個問題。

2.6. EM演算法應用

如果我們從演算法思想的角度來思考EM演算法，我們可以發現我們的演算法里已知的是觀察數據，未知的是隱含數據和模型參數，在E步，我們所做的事情是固定模型參數的值，優化隱含數據的分布，而在M步，我們所做的事情是固定隱含數據分布，優化模型參數的值。EM的應用包括：

支持向量機的SMO演算法
混合高斯模型
K-means
隱馬爾可夫模型

3. EM演算法案例-兩硬幣模型

假設有兩枚硬幣A、B，以相同的概率隨機選擇一個硬幣，進行如下的擲硬幣實驗：共做 5 次實驗，每次實驗獨立的擲十次，結果如圖中 a 所示，例如某次實驗產生了H、T、T、T、H、H、T、H、T、H (H代表正面朝上)。a 是在知道每次選擇的是A還是B的情況下進行，b是在不知道選擇的是A還是B的情況下進行，問如何估計兩個硬幣正面出現的概率？

CASE a

已知每個實驗選擇的是硬幣A 還是硬幣 B，重點是如何計算輸出的概率分布，這其實也是極大似然求導所得。

$egin{align*} underset{ heta }{argmax}logP(Y| heta) &= log(( heta_B^5(1- heta_B)^5) ( heta_A^9(1- heta_A))( heta_A^8(1- heta_A)^2) ( heta_B^4(1- heta_B)^6) ( heta_A^7(1- heta_A)^3) ) \ &= log[( heta_A^{24}(1- heta_A)^6) ( heta_B^9(1- heta_B)^{11}) ] end{align*}$

上面這個式子求導之後發現，5 次實驗中A正面向上的次數再除以總次數作為即為 $hatθ_A$ ，5次實驗中B正面向上的次數再除以總次數作為即為，即:

$hat{ heta}_A = frac{24 }{24+6} = 0.80$

$hat{ heta}_B = frac{9}{ 9 + 11} = 0.45$

CASE b

由於並不知道選擇的是硬幣 A 還是硬幣 B，因此採用EM演算法。

E步：初始化 $hat θ_A^{(0)} = 0.60$ 和 $hat θ_B^{(0)} = 0.50$ ，計算每個實驗中選擇的硬幣是 A 和 B 的概率，例如第一個實驗中選擇 A 的概率為：

$P(z=A|y_1, heta) = frac {P(z=A, y_1| heta)}{P(z=A,y_1| heta) + P(z=B,y_1| heta)} = frac{(0.6)^5*(0.4)^5}{(0.6)^5*(0.4)^5+(0.5)^{10}} = 0.45$

$P(z=B|y_1, heta) = 1- P(z=A|y_1, heta) = 0.55$

計算出每個實驗為硬幣 A 和硬幣 B 的概率，然後進行加權求和。

M步：求出似然函數下界 $Q( heta, heta^i)$ ， $y_j$ 代表第 $j$ 次實驗正面朝上的個數， $mu_j$ 代表第 $j$ 次實驗選擇硬幣 A 的概率， $1-mu_j$ 代表第 $j$ 次實驗選擇硬幣B的概率。

$egin{align*} Q( heta, heta^i) &= sum_{j=1}^5sum_{z} P(z|y_j, heta^i)logP(y_j, z| heta)\&=sum_{j=1}^5 mu_jlog( heta_A^{y_j}(1- heta_A)^{10-y_j}) + (1-mu_j)log( heta_B^{y_j}(1- heta_B)^{10-y_j}) end{align*}$

針對L函數求導來對參數求導，例如對 $θ_A$ 求導：

$egin{align*} frac{partial Q}{partial heta_A} &= mu_1(frac{y_1}{ heta_A}-frac{10-y_1}{1- heta_A}) + cdot cdot cdot + mu_5(frac{y_5}{ heta_A}-frac{10-y_5}{1- heta_A}) = mu_1(frac{y_1 - 10 heta_A} { heta_A(1- heta_A)}) + cdot cdot cdot +mu_5(frac{y_5 - 10 heta_A} { heta_A(1- heta_A)}) \ &= frac{sum_{j=1}^5 mu_jy_j - sum_{j=1}^510mu_j heta_A} { heta_A(1- heta_A)} end{align*}$

求導等於 0 之後就可得到圖中的第一次迭代之後的參數值:

$hat{ heta}_A^{(1)} = 0.71$

$hat{ heta}_B^{(1)} = 0.58$

當然，基於Case a 我們也可以用一種更簡單的方法求得：

$hat{ heta}_A^{(1)} = frac{21.3}{21.3+8.6} = 0.71$

$hat{ heta}_B^{(1)} = frac{11.7}{ 11.7 + 8.4} = 0.58$

第二輪迭代：基於第一輪EM計算好的 $hat{ heta}_A^{(1)}, hat{ heta}_B^{(1)}$ , 進行第二輪 EM，計算每個實驗中選擇的硬幣是 A 和 B 的概率（E步），然後在計算M步，如此繼續迭代......迭代十步之後 $hat{ heta}_A^{(10)} = 0.8, hat{ heta}_B^{(10)} = 0.52$

引用文獻：

1.《從最大似然到EM演算法淺解》

2. Andrew Ng 《Mixtures of Gaussians and the EM algorithm》

3. 《What is the expectation maximization algorithm?》

歡迎關注我的博客：https://2018august.github.io/

人工智慧隨筆?

2018august.github.io

歡迎關注我的微信公眾號：人工智慧隨筆