作為N大機(jī)器學(xué)習(xí)方法的一員，EM算法在各種書籍、博客、網(wǎng)上視頻上被描述或者介紹，每次看完總感覺很多地方含糊不清，不能讓一個初學(xué)者（有一定統(tǒng)計(jì)概率基礎(chǔ)）接受。最近再B站上，看到徐亦達(dá)老師的課程，EM算法這塊講解易于理解和接受，再結(jié)合PRML一書的關(guān)于混合模型和EM章節(jié)內(nèi)容，對整個EM算法從具體的原理上面有了更深入的理解。
在下文中，更多的是通過公式推導(dǎo)和一些文字說明來梳理EM算法，盡量做到大家一看就明白。

極大似然估計(jì)

極大似然估計(jì)是概率統(tǒng)計(jì)中，估計(jì)模型參數(shù)的一種方法，如果對似然概念以及極大似然估計(jì)的概念不理解，可參考維基百科https://zh.wikipedia.org/wiki/%E6%9C%80%E5%A4%A7%E4%BC%BC%E7%84%B6%E4%BC%B0%E8%AE%A1

這里，我們給出一個一維高斯函數(shù)的例子。如圖1所示，一維空間里面離散分布的樣本點(diǎn)其分布特點(diǎn)可以看成是一種高斯分布。那這些樣本的高斯分布函數(shù)參數(shù)怎么求解呢？可以通過極大似然估計(jì)獲得。

圖1 一維高斯示例

假設(shè)我們獲取圖1中數(shù)據(jù)樣本集合為 $\{x_1,x_2,..x_i,...,x_n\}$ ，其分布函數(shù)假設(shè)為一維高斯分布，即：
$p(x) =\frac{1}{\sqrt{2\pi}\sigma^2}exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$
那所有數(shù)據(jù)的聯(lián)合概率,即似然函數(shù)為：
$P(\boldsymbol{x})=\prod_{i=1}^n p(x_i)=\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma}exp\left(-\frac{(x_i-\mu)^2}{2\sigma^2}\right)$
對等式兩邊求log，即log-likelihood：
$Ln(P(\boldsymbol{x}))=-nLn(\sqrt{2\pi}\sigma)-\sum_{i=1}^n(x_i-\mu)^2/(2\sigma^2)$
分別對參數(shù)求導(dǎo)：
$\frac{\partial{Ln(P(\boldsymbol{x}))}}{\partial\mu}=\frac{\sum_{i=1}^n(x_i-\mu)}{\sigma^2}=0$
$\frac{\partial{Ln(P(\boldsymbol{x}))}}{\partial\sigma^2}=-\frac{n}{2\sigma^2}+\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^4}=0$
可以得到：
$\mu=\frac{1}{n}\sum_{i=1}^n{x_i}$
$\sigma^2=\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2$
通過上述方法，就可以得到基于樣本數(shù)據(jù)和假設(shè)分布函數(shù)模型情況下，得到樣本數(shù)據(jù)的分布函數(shù)。

高斯混合模型

從圖2所示中，如果樣本數(shù)據(jù)分布式藍(lán)色點(diǎn)和橙色點(diǎn)的集合，如果依然用高斯分布去進(jìn)行最大似然估計(jì)，得到的分布模型自然是不合適的。很顯然，樣本數(shù)據(jù)分布存在兩個密集區(qū)（藍(lán)色點(diǎn)和橙色點(diǎn)），如果能通過一種方法，確認(rèn)樣本數(shù)據(jù)里面的藍(lán)色點(diǎn)和橙色點(diǎn)，然后分別對藍(lán)色點(diǎn)集合進(jìn)行一個高斯估計(jì)，對橙色點(diǎn)集進(jìn)行另外一個高斯估計(jì)，兩個高斯混合到一起，是不是比單個高斯能更好的匹配樣本數(shù)據(jù)的分布情況。這種情況下的分布函數(shù)就是高斯混合模型。

圖2 一維混合高斯分布示意

這里我們先直接給出高斯混合模型的分布函數(shù)(多維)：
$p(\boldsymbol{x})=\sum_{k=1}^K\pi_k N (\boldsymbol{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)$
在圖2例子中，提到如果給每一個數(shù)據(jù)樣本能夠先進(jìn)行分類，即確定屬于哪一個集中的簇中，就能比較容易的進(jìn)行混合模型的求解。這說明了什么呢？我們可以理解，原始樣本數(shù)據(jù)是不完全的（incomplete），引入一個K維的二值隨機(jī)變量 $\boldsymbol z$ ，這個變量采用"1-of-K"的表示方法，即K維中，只有一維是1，其余所有元素等于0。那么每一個樣本數(shù)據(jù)，即有數(shù)據(jù)特征 $\boldsymbol x$ ，再匹配一個分配變量 $\boldsymbol z$ ，就可以將圖2過程完整描述，即我們認(rèn)為 $\boldsymbol x$ 和 $\boldsymbol z$ 聯(lián)合在一起才是完全的（complete）。

數(shù)學(xué)表示上，我們利用邊緣概率分布 $p(\boldsymbol z)$ 和條件概率分布 $p(\boldsymbol x|\boldsymbol z)$ 定義聯(lián)合概率分布 $p(\boldsymbol x,\boldsymbol z)$ .
$\boldsymbol z$ 的邊緣概率分布根據(jù)混合系數(shù) $\pi_k$ 進(jìn)行賦值，即 $p(z_k = 1)=\pi_k,0\leq\pi_k\leq 1,\sum_{k=1}^K\pi_k=1$ ,使得邊緣概率分布是合法，那么：
$p(\boldsymbol x) = \prod_{k=1}^K\pi_k^{z_k}$
類似的，在給定 $\boldsymbol z$ 的一個特定的值，即針對某一簇的情況， $\boldsymbol x$ 的條件概率分布是一個高斯分布，即 $p(\boldsymbol x|z_k = 1)=N(\boldsymbol x|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)$
$p(\boldsymbol x|\boldsymbol z) = \prod_{k=1}^KN(\boldsymbol x|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)^{z_k}$
那么根據(jù)貝葉斯公式得到高斯混合模型：
$p(\boldsymbol x) = \sum_{\boldsymbol z}p(\boldsymbol z)p(\boldsymbol x|\boldsymbol z)=\sum_{k=1}^K\pi_k\mathcal N(\boldsymbol x|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)$
由于我們只能觀察樣本數(shù)據(jù) $\boldsymbol x_1,...,\boldsymbol x_n$ ,無法直接獲取每個數(shù)據(jù)的分配變量 $\boldsymbol z_1,...,\boldsymbol z_n$ ,可理解 $\boldsymbol z_1,...,\boldsymbol z_n$ 為潛在變量（latent）
依據(jù)極大似然函數(shù)的求解方法，我們可以對高斯混合模型寫出數(shù)據(jù)的對數(shù)似然函數(shù)：
$Ln(p(\boldsymbol X))=\sum_{i=1}^nLn\left\{ \sum_{k=1}^K\pi_k N(\boldsymbol x|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k)\right\}$

由于對數(shù)函數(shù)內(nèi)部出現(xiàn)求和，這種情況是沒有辦法通過求導(dǎo)等于0的方式獲取估計(jì)參數(shù)的解析解的。這種情況下，就需要用到EM算法，通過迭代方式獲取估計(jì)參數(shù)。下面我們先從一般性問題上進(jìn)行EM算法的理論描述，然后再利用EM算法推導(dǎo)高斯混合模型的計(jì)算方法。

EM算法理論證明

EM算法叫做期望最大化方法，首先我們給出EM算法一般性結(jié)論或者說步驟，其具體分為兩步，即E-step和M-step。

E-step：是求期望的步驟，求誰的期望呢？對聯(lián)合概率分布的log數(shù)據(jù)利用潛在變量的條件概率函數(shù)求期望，這個期望定義為 $Q(\theta,\theta^{(g)})$ ，我們寫連續(xù)形式，離散形式類似。
$Q(\theta,\theta^{(g)})=\int_{\boldsymbol z}log(p_{\theta}(\boldsymbol x,\boldsymbol z))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
注：有些書籍或者博客中會將參數(shù)寫在 $p$ 里面，我個人看著比較不習(xí)慣，總會把變量和參數(shù)搞混了，這里寫在下標(biāo)上。但如果有了解貝葉斯估計(jì)的，那個時候參數(shù)也被考慮成變量，其有自己的分布。表達(dá)只是一種形式，背后的理論理解掌握就OK啦。
M-step:是求最大的步驟，求誰的最大呢？自然是期望 $Q$ 的最大，從而得到新的估計(jì)參數(shù)，即
$\theta^{(g+1)}=\mathop{\arg\min}_{\theta}Q(\theta,\theta^{(g)})$

EM算法的步驟，通過高斯混合模型可以直觀理解記憶。 $p(\boldsymbol z|\boldsymbol x)$ 是什么意思呢，其含義是在給定數(shù)據(jù)樣本的情況下，潛在變量的概率情況。也就是說在高斯混合模型中，給定樣本數(shù)據(jù)，我們推測樣本屬于哪一個高斯簇的概率情況。在確定分配情況后，每一個簇都用高斯進(jìn)行估計(jì)，衡量估計(jì)好還是不好，用期望形式表示。這樣可以幫助理解和記憶Q的定義。那EM算法怎么證明其收斂性呢？

我們要保證：
$p_{\theta^{(g+1)}}(\boldsymbol x)\geq p_ {\theta^{(g)}}(\boldsymbol x)$
這樣每次迭代可以使得似然函數(shù)隨著參數(shù)變大，一直到似然函數(shù)不再增大或滿足其他終止條件止。

那怎么保證呢？首先，利用貝葉斯公式有：
$p_ {\theta}(\boldsymbol x)=\frac{p_\theta(\boldsymbol x,\boldsymbol z)}{p_\theta(\boldsymbol z|\boldsymbol x)}$

兩邊同時取log
$log(p_ {\theta}(\boldsymbol x))=log(p_\theta(\boldsymbol x,\boldsymbol z))-log(p_\theta(\boldsymbol z|\boldsymbol x))$
然后兩邊同時用 $p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)$ 求期望，可以得：
$\int_{\boldsymbol z}log(p_ {\theta}(\boldsymbol x))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z=\int_{\boldsymbol z}log(p_\theta(\boldsymbol x,\boldsymbol z))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z-\int_{\boldsymbol z}log(p_\theta(\boldsymbol z|\boldsymbol x))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
等式左邊 $p_ {\theta}(\boldsymbol x)$ 和 $\boldsymbol z$ 沒有關(guān)系， $p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)$ 是概率密度函數(shù)，積分是1，所以等式左邊依然是 $log(p_ {\theta}(\boldsymbol x))$ .
等式右邊第一項(xiàng)就是E-step里面的Q函數(shù)，第二項(xiàng)我們記做H函數(shù)，則：
$log(p_ {\theta}(\boldsymbol x))=Q(\theta,\theta^{(g)})-H(\theta,\theta^{(g)})$
要保證 $p_{\theta^{(g+1)}}(\boldsymbol x)\geq p_ {\theta^{(g)}}(\boldsymbol x)$ ，首先 $Q(\theta^{(g+1)},\theta^{(g)})\geq Q(\theta^{(g)},\theta^{(g)})$ ，那么是不是保證 $H(\theta^{(g+1)},\theta^{(g)})\leq H(\theta^{(g)},\theta^{(g)})$ 滿足，就一定有 $p_{\theta^{(g+1)}}(\boldsymbol x)\geq p_ {\theta^{(g)}}(\boldsymbol x)$ ？答案是肯定的。（說明一下，這里面的H和Q函數(shù)都是關(guān)于 $\theta$ 的函數(shù)，每一次迭代 $\theta^{(g)}$ 是已知的，他不是變量）

那下面只要證明 $H(\theta^{(g+1)},\theta^{(g)})\leq H(\theta^{(g)},\theta^{(g)})$ 就可以了。

$H(\theta^{(g)},\theta^{(g)})-H(\theta,\theta^{(g)})$
$=\int_{\boldsymbol z}log(p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z-\int_{\boldsymbol z}log(p_\theta(\boldsymbol z|\boldsymbol x))p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
$=-\int_{\boldsymbol z}log\left(\frac{p_\theta(\boldsymbol z|\boldsymbol x)}{p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)}\right)p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
$\geq-log\left(\int_{\boldsymbol z}\frac{p_\theta(\boldsymbol z|\boldsymbol x)}{p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)}p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z\right )$
$=-log\left(\int_{\boldsymbol z}p_\theta(\boldsymbol z|\boldsymbol x)d\boldsymbol z\right)=-log(1)=0$

因此可以得到 $H(\theta^{(g+1)},\theta^{(g)})\leq H(\theta^{(g)},\theta^{(g)})$ ，則整個EM算法的收斂性證畢。

注：這里用到了Jessian不等式，如果函數(shù)是convex的，則有函數(shù)的期望大于期望的函數(shù)，即 $E[f(x)]\geq f(E[x])$ .

圖3 Jessian不等式示例

EM算法與ELOB和KL散度

上述EM算法的證明，有一些技巧性，而這些技巧性有沒有一種是在已知結(jié)論基礎(chǔ)上硬湊的感覺呢，尤其是用 $p_{\theta^{(g)}}(\boldsymbol z| \boldsymbol x)$ 去求期望，就是為了去構(gòu)造Q函數(shù)。那有沒有更好的解釋或者更為直觀的方式來得到相同的結(jié)論呢？答案是有的。

首先，仍然用到：
$log(p_ {\theta}(\boldsymbol x))=log(p_\theta(\boldsymbol x,\boldsymbol z))-log(p_\theta(\boldsymbol z|\boldsymbol x))$
我們稍微變個型，假設(shè)一個概率密度函數(shù) $q(\boldsymbol z)$ :
$log(p_ {\theta}(\boldsymbol x))=log\left(\frac{p_\theta(\boldsymbol x,\boldsymbol z)}{q(\boldsymbol z)}\right)-log\left(\frac{p_\theta(\boldsymbol z|\boldsymbol x)}{q(\boldsymbol z)}\right)$
兩邊同時用 $q(\boldsymbol z)$ 求期望：
$log(p_ {\theta}(\boldsymbol x))=\int_{\boldsymbol z}log\left(\frac{p_\theta(\boldsymbol x,\boldsymbol z)}{q(\boldsymbol z)}\right)q(\boldsymbol z)d\boldsymbol z-\int_{\boldsymbol z}log\left(\frac{p_\theta(\boldsymbol z|\boldsymbol x)}{q(\boldsymbol z)}\right)q(\boldsymbol z)d\boldsymbol z$

其中等式右邊第一項(xiàng)，我們記做 $\mathcal L(q,\theta)$ ，可以稱為ELOB，EvidenceLowerBound，第二項(xiàng)是 $q(\boldsymbol z)$ 和 $p_\theta(\boldsymbol z|\boldsymbol x)$ 的KL散度。

圖4 似然函數(shù)與ELOB和KL散度

如圖4所示，當(dāng)我固定參數(shù) $\theta$ 時候，ELOB就是關(guān)于 $q$ 的泛函（只聽過沒學(xué)過，可理解為函數(shù)的函數(shù)），那ELOB的上界是什么呢？那首先要理解KL散度，KL散度是衡量兩個概率分布之間的差異性，如果兩個分布沒有差異，則KL散度等于0，如果有差異則大于0，所以KL散度的最小值就是0，那ELOB的上界就是此刻的似然函數(shù)。

在EM算法中，當(dāng)前迭代時，參數(shù) $\theta^{(g)}$ ，則對于E-step，我們需要使得ELOB最大，即KL散度為0，如圖5所示，其中 $\theta^{old}$ 為 $\theta^{(g)}$ 。此時， $q(\boldsymbol z) = p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)$

圖5 EM算法 E-step

對于M-Step，我們是需要得到新的估計(jì)參數(shù)，這個時候，固定 $q(\boldsymbol z)$ ,重新獲得ELOB的最大值，這個時候的ELOB是什么樣的呢？
$\mathcal L(q,\theta)=\int_{\boldsymbol z}log\left(\frac{p_\theta(\boldsymbol x,\boldsymbol z)}{q(\boldsymbol z)}\right)q(\boldsymbol z)d\boldsymbol z$
$=\int_{\boldsymbol z}log\left(p_\theta(\boldsymbol x,\boldsymbol z)\right)p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z-\int_{\boldsymbol z}log\left(p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)\right)p_{\theta^{(g)}}(\boldsymbol z|\boldsymbol x)d\boldsymbol z$
右邊第二項(xiàng)沒有參數(shù) $\theta$ ，所以固定 $q$ 最大化ELOB，就等于最大化第一項(xiàng)，而第一項(xiàng)是什么呢？就是 $Q$ 函數(shù)。

圖6 EM算法M-step

如圖6所示，我們最大化 $Q$ 函數(shù)，也就是最大化此次迭代的ELOB。在新的參數(shù)下， $\mathcal L(q,\theta^{(g+1)}) \geq \mathcal L(q,\theta^{(g)})$ ，此刻 $q$ 不變，所以和新的 $p_{\theta^{(g+1)}}(\boldsymbol z|\boldsymbol x)$ 在沒有達(dá)到局部（或者全局）極大值的情況下，是兩個不同的概率分布，即二者KL散度是大于0的，那此刻的似然函數(shù)等于此刻的KL散度加上此刻的ELOB，自然是有 $p_{\theta^{(g+1)}}(\boldsymbol x)\geq p_ {\theta^{(g)}}(\boldsymbol x)$ 。

從ELOB和KL散度的層面可以更好的理解EM算法的迭代過程。

PRML一書中，有圖7所示的示意圖，能比較形象的說明，EM算法的迭代過程。

圖7 EM算法迭代示意

圖7中，紅色曲線表示(不完整數(shù)據(jù))對數(shù)似然函數(shù)，它的最大值是我們想要得到的。我們首先選擇某個初始的參數(shù)值 $\theta^{old}$ ,然后在第一個E步驟中,我們計(jì)算潛在變量上的后驗(yàn)概率分布,得到了 $\mathcal L(q, \theta^{old})$ 的一個更小的下界，它的值等于在 $\theta^{old}$ 處的對數(shù)似然函數(shù)值，用藍(lán)色曲線表示。注意，下界與對數(shù)似然函數(shù)在 $\theta^{old}$ 處以切線的方式連接，因此兩條曲線的梯度相同。這個界是一個凹函數(shù)，對于指數(shù)族分布的混合分布來說，有唯一的最大值。在M步驟中,下界被最大化，得到了新的值 $\theta^{new}$ ,這個值給出了比 $\theta^{old}$ 處更大的對數(shù)似然函數(shù)值。接下來的E步驟構(gòu)建了一個新的下界,它在 $\theta^{new}$ 處與對數(shù)似然函數(shù)切線連接,用綠色曲線表示。以這樣的方式不停迭代，直到滿足條件。