原貼地址：https://blog.csdn.net/lin_limin/article/details/81048411
最近在看曉川老（shi）師(shu)的博士論文，接觸了混合高斯模型（Gaussian mixture model, GMM）和EM(Expectation Maximization)算法，不禁被論文中龐大的數(shù)學(xué)公式所嚇退。本文通過(guò)查閱相關(guān)資料，在復(fù)雜巧妙的推理公式中融入了自己的理解，詳細(xì)梳理了混合高斯模型和EM算法。
1 單高斯模型(Gaussian single model, GSM）
簡(jiǎn)單回顧一下概率論講過(guò)的高斯模型。
高斯模型是一種常用的變量分布模型，在數(shù)理統(tǒng)計(jì)領(lǐng)域有著廣泛的應(yīng)用（……好吧讀了這么多年書(shū)沒(méi)白費(fèi)，教科書(shū)般的話語(yǔ)已植入骨髓）。一維高斯分布的概率密度函數(shù)如下：
f(x)=12π√σexp(?(x?μ)22σ2)(1)f(x) = \frac{1}{{\sqrt {2\pi } \sigma }}\exp ( - \frac{{{{(x - \mu )}^2}}}{{2{\sigma ^2}}})\tag{1}f(x)=2π?σ1?exp(?2σ2(x?μ)2?)(1)
?μ\muμ和σ2{\sigma ^2}σ2分別是高斯分布的均值和方差。
譬如將男生身高視為變量X, 假設(shè)男生的身高服從高斯分布，則X～N(μ,σ2)X \sim N(\mu ,{\sigma ^2})X～N(μ,σ2)，女生亦如此。只是男女生身高分布可能具有不同的均值和方差。圖1是從谷歌圖庫(kù)中搜索到的男女生身高分布圖，來(lái)源不清，個(gè)人覺(jué)得男生的均值身高虛高……四個(gè)記號(hào)分別表示0~3σ\sigmaσ準(zhǔn)則。
圖1 男女生身高分布差異
多維變量X=(x1,x2,...xn)X = ({x_1},{x_2},...{x_n})X=(x1?,x2?,...xn?)的聯(lián)合概率密度函數(shù)為：
f(X)=1(2π)d/2∣Σ∣1/2exp[?12(X?u)TΣ?1(X?u)],X=(x1,x2...xn)(2)f(X) = \frac{1}{{{{(2\pi )}^{d/2}}{{\left| \Sigma \right|}^{1/2}}}}\exp [ - \frac{1}{2}{(X - u)^T}{\Sigma ^{ - 1}}(X - u)],X = ({x_1},{x_2}...{x_n})\tag{2}f(X)=(2π)d/2∣Σ∣1/21?exp[?21?(X?u)TΣ?1(X?u)],X=(x1?,x2?...xn?)(2)
其中：
d：變量維度。對(duì)于二維高斯分布，有d=2;
?u=?????u1u2...un?????u = \left( \begin{array}{l}{u_1}\\{u_2}\\...\\{u_n}\end{array} \right)u=?????u1?u2?...un???????：各維變量的均值；
?Σ\SigmaΣ：協(xié)方差矩陣，描述各維變量之間的相關(guān)度。對(duì)于二維高斯分布，有：
Σ=[δ11δ21δ12δ22](3)\Sigma= \left[\begin{matrix}\delta _{11} & \delta _{12}\\\delta _{21}& \delta _{22} \end{matrix}\right] \tag{3}Σ=[δ11?δ21??δ12?δ22??](3)
圖2 二維高斯數(shù)據(jù)分布
圖2是二維高斯分布產(chǎn)生的數(shù)據(jù)示例，參數(shù)設(shè)定為：u=(00),Σ=[10.80.85]u = \left( \begin{array}{l}0\\0\end{array} \right),\Sigma= \left[\begin{matrix}1 & 0.8\\0.8&5 \end{matrix}\right]u=(00?),Σ=[10.8?0.85?]。關(guān)于二維高斯分布的參數(shù)設(shè)定對(duì)為高斯曲面的影響，這篇文章二維高斯分布（Two-dimensional Gaussian distribution）的參數(shù)分析有提及，主要是為了下文理解混合高斯分布做鋪墊。服從二維高斯分布的數(shù)據(jù)主要集中在一個(gè)橢圓內(nèi)部，服從三維的數(shù)據(jù)集中在一個(gè)橢球內(nèi)部。
2 混合高斯模型（Gaussian mixture model, GMM）
先來(lái)看一組數(shù)據(jù)。
圖3 混合高斯分布所產(chǎn)生數(shù)據(jù)
? 如果我們假設(shè)這組數(shù)據(jù)是由某個(gè)高斯分布產(chǎn)生的，利用極大似然估計(jì)（后文還會(huì)提及）對(duì)這個(gè)高斯分布做參數(shù)估計(jì)，得到一個(gè)最佳的高斯分布模型如下。
圖4 用單高斯模型對(duì)樣本作分析不合理示意
有什么問(wèn)題嗎？一般來(lái)講越靠近橢圓的中心樣本出現(xiàn)的概率越大，這是由概率密度函數(shù)決定的，但是這個(gè)高斯分布的橢圓中心的樣本量卻極少。顯然樣本服從單高斯分布的假設(shè)并不合理。單高斯模型無(wú)法產(chǎn)生這樣的樣本。
實(shí)際上，這是用兩個(gè)不同的高斯分布模型產(chǎn)生的數(shù)據(jù)。
圖5 混合高斯模型對(duì)樣本作分析示意
正當(dāng)單高斯模型抓耳撓腮的時(shí)候，混合高斯模型就大搖大擺地進(jìn)場(chǎng)了。它通過(guò)求解兩個(gè)高斯模型，并通過(guò)一定的權(quán)重將兩個(gè)高斯模型融合成一個(gè)模型，即最終的混合高斯模型。這個(gè)混合高斯模型可以產(chǎn)生這樣的樣本。
更一般化的描述為：假設(shè)混合高斯模型由K個(gè)高斯模型組成（即數(shù)據(jù)包含K個(gè)類），則GMM的概率密度函數(shù)如下：
p(x)=∑Kk=1p(k)p(x∣k)=∑Kk=1πkN(x∣uk,Σk)(3)p(x)=\sum\limits_{k = 1}^K {p(k)p(x|k) = } \sum\limits_{k = 1}^K {{\pi _k}N(x|{u_k},{\Sigma _k})} \tag{3}p(x)=k=1∑K?p(k)p(x∣k)=k=1∑K?πk?N(x∣uk?,Σk?)(3)
其中，p(x∣k)=N(x∣uk,Σk)p(x|k) = N(x|{u_k},{\Sigma _k})p(x∣k)=N(x∣uk?,Σk?)是第k個(gè)高斯模型的概率密度函數(shù)，可以看成選定第k個(gè)模型后，該模型產(chǎn)生x的概率；p(k)=πkp(k) = {\pi _k}p(k)=πk?是第k個(gè)高斯模型的權(quán)重，稱作選擇第k個(gè)模型的先驗(yàn)概率，且滿足∑Kk=1πk=1\sum\limits_{k = 1}^K {{\pi _k}} = 1k=1∑K?πk?=1。
所以，混合高斯模型并不是什么新奇的東西，它的本質(zhì)就是融合幾個(gè)單高斯模型，來(lái)使得模型更加復(fù)雜，從而產(chǎn)生更復(fù)雜的樣本。理論上，如果某個(gè)混合高斯模型融合的高斯模型個(gè)數(shù)足夠多，它們之間的權(quán)重設(shè)定得足夠合理，這個(gè)混合模型可以擬合任意分布的樣本。
下面通過(guò)幾張圖片來(lái)幫助理解混合高斯模型。
首先從簡(jiǎn)單的一維混合高斯模型說(shuō)起。
圖6 一維混合高斯模型
在圖6中，y1,y2和y3分別表示三個(gè)一維高斯模型，他們的參數(shù)設(shè)定如圖所示。y4表示將三個(gè)模型的概率密度函數(shù)直接相加，注意的是這并不是一個(gè)混合高斯模型，因?yàn)椴粷M足∑Kk=1πk=1\sum\limits_{k = 1}^K {{\pi _k}} = 1k=1∑K?πk?=1的條件。而y5和y6分別是由三個(gè)相同的高斯模型融合生成的不同混合模型。由此可見(jiàn)，調(diào)整權(quán)重將極大影響混合模型的概率密度函數(shù)曲線。另一方面也可以直觀地理解混合高斯模型可以更好地?cái)M合樣本的原因：它有更復(fù)雜更多變的概率密度函數(shù)曲線。理論上，混合高斯模型的概率密度函數(shù)曲線可以是任意形狀的非線性函數(shù)。
下面再給出一個(gè)二維空間3個(gè)高斯模型混合的例子。
(a) 3個(gè)類別高斯分布截面輪廓線
(b) 混合高斯分布截面輪廓線
? 二維混合高斯分布概率密度函數(shù)圖
圖7 二維混合高斯模型
(a) 圖表示的是3個(gè)高斯模型的截面輪廓圖，3個(gè)模型的權(quán)重系數(shù)已在圖中注明，由截面輪廓圖可知3個(gè)模型之間存在很多重疊區(qū)域。其實(shí)這也正是混合高斯模型所希望的。因?yàn)槿绻鼈冎g的重疊區(qū)域較少，那么生成的混合高斯模型一般較為簡(jiǎn)單，難以生成較為復(fù)雜的樣本。
設(shè)定好了3個(gè)高斯模型和它們之間的權(quán)重系數(shù)之后，就可以確定二維混合高斯分布的概率密度函數(shù)曲面，如圖?所示。圖(b)是對(duì)于圖?概率密度曲面的截面輪廓線。從圖7也可以看出，整個(gè)混合高斯分布曲面相對(duì)比于單高斯分布曲面已經(jīng)異常復(fù)雜。實(shí)際上，通過(guò)調(diào)整混合高斯分布的系數(shù)(π,μ,Σ)(\pi ,\mu ,\Sigma )(π,μ,Σ)，可以使得圖?的概率密度曲面去擬合任意的三維曲面，從而采樣生成所需要的數(shù)據(jù)樣本。
3 極大似然估計(jì)(Maximum Likehood Estimate, MLE)（最大化對(duì)數(shù)似然函數(shù)）
● 最大化對(duì)數(shù)似然函數(shù)（log-likelihood function）的意義
首先直觀化地解釋一下最大化對(duì)數(shù)似然函數(shù)要解決的是什么問(wèn)題。
假設(shè)我們采樣得到一組樣本yt{y_t}yt?,而且我們知道變量Y服從高斯分布（本文只提及高斯分布，其他變量分布模型類似），數(shù)學(xué)形式表示為Y～N(μ,Σ)Y \sim N(\mu ,\Sigma )Y～N(μ,Σ)。采樣的樣本如圖8所示，我們的目的就是找到一個(gè)合適的高斯分布（也就是確定高斯分布的參數(shù)μ,Σ\mu ,\Sigmaμ,Σ），使得這個(gè)高斯分布能產(chǎn)生這組樣本的可能性盡可能大。
圖8 最大化似然函數(shù)的意義
那怎么找到這個(gè)合適的高斯分布呢（在圖8中的表示就是1~4哪個(gè)分布較為合適）？這時(shí)候似然函數(shù)就閃亮登場(chǎng)了。
似然函數(shù)數(shù)學(xué)化：設(shè)有樣本集Y=y1,y2...yNY = {y_1},{y_2}...{y_N}Y=y1?,y2?...yN?。p(yn∣μ,Σ)p({y_n}|\mu ,\Sigma )p(yn?∣μ,Σ)是高斯分布的概率分布函數(shù)，表示變量Y=ynY = {y_n}Y=yn?的概率。假設(shè)樣本的抽樣是獨(dú)立的，那么我們同時(shí)抽到這N個(gè)樣本的概率是抽到每個(gè)樣本概率的乘積，也就是樣本集Y的聯(lián)合概率。此聯(lián)合概率即為似然函數(shù)：
L(μ,Σ)=L(y1,y2...yN;μ,Σ)=∏Nn=1p(yn;μ,Σ)(4)L(\mu ,\Sigma ) = L({y_1},{y_2}...{y_N};\mu ,\Sigma ) = \prod\limits_{n = 1}^N {p({y_n};\mu ,\Sigma )}\tag{4}L(μ,Σ)=L(y1?,y2?...yN?;μ,Σ)=n=1∏N?p(yn?;μ,Σ)(4)
對(duì)式子(4)進(jìn)行求導(dǎo)并令導(dǎo)數(shù)為0（即最大化似然函數(shù)，一般還會(huì)先轉(zhuǎn)化為對(duì)數(shù)似然函數(shù)再最大化），所求出的參數(shù)就是最佳的高斯分布對(duì)應(yīng)的參數(shù)。
所以最大化似然函數(shù)的意義就是：通過(guò)使得樣本集的聯(lián)合概率最大來(lái)對(duì)參數(shù)進(jìn)行估計(jì)，從而選擇最佳的分布模型。
對(duì)于圖8產(chǎn)生的樣本用最大化似然函數(shù)的方法，最終可以得到序號(hào)1對(duì)應(yīng)的高斯分布模型是最佳的模型。
4.1 為什么要有EM算法（EM算法與極大似然估計(jì)分別適用于什么問(wèn)題）
● 嘗試用極大似然估計(jì)的方法來(lái)解GMM模型
解GMM模型，實(shí)際上就是確定GMM模型的參數(shù)(μ,Σ,π)(\mu ,\Sigma ,\pi )(μ,Σ,π)，使得由這組參數(shù)確定的GMM模型最有可能產(chǎn)生采樣的樣本。
先試試看用極大似然估計(jì)的方法來(lái)解GMM模型會(huì)出現(xiàn)什么樣的問(wèn)題。
如第3小節(jié)所述，要利用極大似然估計(jì)求解模型最重要的一步就是求出似然函數(shù)，即樣本集出現(xiàn)的聯(lián)合概率。而對(duì)于混合高斯模型，如何求解某個(gè)樣本yt{y_t}yt?的概率？顯然我們得先知道這個(gè)樣本來(lái)源于哪一類高斯模型，然后求這個(gè)高斯模型生成這個(gè)樣本的概率p(yt)p({y_t})p(yt?)。
但是問(wèn)題來(lái)了：我們只有樣本。不知道樣本到底來(lái)源于哪一類的高斯模型。那么如何求解樣本的生成概率p(yt)p({y_t})p(yt?)？
先引入一個(gè)隱變量γ\gammaγ。它是一個(gè)K維二值隨機(jī)變量，在它的K維取值中只有某個(gè)特定的元素γk{\gamma _k}γk?的取值為1，其它元素的取值為0。實(shí)際上，隱變量描述的就是：每一次采樣，選擇第k個(gè)高斯模型的概率，故有：
p(γk=1)=πk(5)p({\gamma _k} = 1) = {\pi _k}\tag{5}p(γk?=1)=πk?(5)
當(dāng)給定了γ\gammaγ的一個(gè)特定的值之后（也就是知道了這個(gè)樣本從哪一個(gè)高斯模型進(jìn)行采樣），可以得到樣本y的條件分布是一個(gè)高斯分布，滿足：
p(y∣γk=1)=N(y∣μk,Σk)(6)p(y|{\gamma _k} = 1) = N(y|{\mu _k},{\Sigma _k})\tag{6}p(y∣γk?=1)=N(y∣μk?,Σk?)(6)
而實(shí)際上，每個(gè)樣本到底是從這K個(gè)高斯模型中哪個(gè)模型進(jìn)行采樣的，是都有可能的。故樣本y的概率為：
p(y)=∑γp(γ)p(y∣γ)=∑Kk=1πkN(y∣μk,Σk)(7)p(y) = \sum\nolimits_\gamma {p(\gamma )} p(y|\gamma ){\rm{ = }}\sum\limits_{{\rm{k}} = 1}^K {{\pi _k}N(y|{\mu _k},{\Sigma _k})} \tag{7}p(y)=∑γ?p(γ)p(y∣γ)=k=1∑K?πk?N(y∣μk?,Σk?)(7)
樣本集Y(n個(gè)樣本點(diǎn))的聯(lián)合概率為：
L(μ,Σ,π)=L(y1,y2...yN;μ,Σ,π)=∏Nn=1p(yn;μ,Σ,π)=∏Nn=1∑Kk=1πkN(yn∣μk,Σk)(8)L(\mu ,\Sigma ,\pi ) = L({y_1},{y_2}...{y_N};\mu ,\Sigma ,\pi ) = \prod\limits_{n = 1}^N {p({y_n};\mu ,\Sigma ,\pi )} = \prod\limits_{n = 1}^N {\sum\limits_{{\rm{k}} = 1}^K {{\pi _k}N({y_n}|{\mu _k},{\Sigma _k})} } \tag{8}L(μ,Σ,π)=L(y1?,y2?...yN?;μ,Σ,π)=n=1∏N?p(yn?;μ,Σ,π)=n=1∏N?k=1∑K?πk?N(yn?∣μk?,Σk?)(8)
對(duì)數(shù)似然函數(shù)表示為：
lnL(μ,Σ,π)=∑Nn=1ln∑Kk=1πkN(yn∣μk,Σk)(9)\ln L(\mu ,\Sigma ,\pi ) = \sum\limits_{n = 1}^N {\ln \sum\limits_{{\rm{k}} = 1}^K {{\pi _k}N({y_n}|{\mu _k},{\Sigma _k})} } \tag{9}lnL(μ,Σ,π)=n=1∑N?lnk=1∑K?πk?N(yn?∣μk?,Σk?)(9)
好了，然后求導(dǎo)，令導(dǎo)數(shù)為0，得到模型參數(shù)(μ,Σ,π)(\mu ,\Sigma ,\pi )(μ,Σ,π)。
貌似問(wèn)題已經(jīng)解決了，喜大普奔。
然而仔細(xì)觀察可以發(fā)現(xiàn)，對(duì)數(shù)似然函數(shù)里面，對(duì)數(shù)里面還有求和。實(shí)際上沒(méi)有辦法通過(guò)求導(dǎo)的方法來(lái)求這個(gè)對(duì)數(shù)似然函數(shù)的最大值。
MLE（極大似然估計(jì)）略顯沮喪。這時(shí)候EM算法走過(guò)來(lái)，安慰著說(shuō)：兄弟別哭，老哥幫你。
下面先闡述一下極大似然估計(jì)與EM算法分別適用于解決什么樣的問(wèn)題。
圖9 極大似然估計(jì)適用問(wèn)題
圖10 EM算法適用問(wèn)題
如果我們已經(jīng)清楚了某個(gè)變量服從的高斯分布，而且通過(guò)采樣得到了這個(gè)變量的樣本數(shù)據(jù)，想求高斯分布的參數(shù)，這時(shí)候極大似然估計(jì)可以勝任這個(gè)任務(wù)；而如果我們要求解的是一個(gè)混合模型，只知道混合模型中各個(gè)類的分布模型（譬如都是高斯分布）和對(duì)應(yīng)的采樣數(shù)據(jù)，而不知道這些采樣數(shù)據(jù)分別來(lái)源于哪一類（隱變量），那這時(shí)候就可以借鑒EM算法。EM算法可以用于解決數(shù)據(jù)缺失的參數(shù)估計(jì)問(wèn)題（隱變量的存在實(shí)際上就是數(shù)據(jù)缺失問(wèn)題，缺失了各個(gè)樣本來(lái)源于哪一類的記錄）。
下面將介紹EM算法的兩個(gè)步驟：E-step(expectation-step，期望步)和M-step(Maximization-step,最大化步);
我們現(xiàn)有樣本集Y=(y1,y2...yT)Y = ({y_1},{y_2}...{y_T})Y=(y1?,y2?...yT?)，通過(guò)隱變量γt,k{\gamma _{t,k}}γt,k?（表示yt{y_t}yt?這個(gè)樣本來(lái)源于第k個(gè)模型）的引入，可以將數(shù)據(jù)展開(kāi)成完全數(shù)據(jù)：
(yt,γt,1,γt,2...γt,K),t=1,2...T({y_t},{\gamma _t}_{,1},{\gamma _{t,2}}...{\gamma _{t,K}}),t = 1,2...T(yt?,γt?,1?,γt,2?...γt,K?),t=1,2...T
所謂的完全數(shù)據(jù)，就是不缺失的數(shù)據(jù)。只有樣本集Y=(y1,y2...yT)Y = ({y_1},{y_2}...{y_T})Y=(y1?,y2?...yT?)的數(shù)據(jù)是不完整的，存在信息缺失的。若yt{y_t}yt?由第1類采樣而來(lái)，則有γt,1=1,γt,2=0...γt,K=0{\gamma _t}_{,1} = 1,{\gamma _{t,2}} = 0...{\gamma _{t,K}} = 0γt?,1?=1,γt,2?=0...γt,K?=0，表示為(yt,1,0,...0)({y_t},1,0,...0)(yt?,1,0,...0)。
所以要求能采到這組數(shù)據(jù)的可能性，需要分兩步走：①第t個(gè)樣本由哪一類產(chǎn)生？②如果第t個(gè)樣本由第k類產(chǎn)生，那么第k類產(chǎn)生第t個(gè)樣本的概率為多少？
綜合考慮上面兩步，有了完全數(shù)據(jù)的似然函數(shù)：
p(y,γ∣μ,Σ,π)=∏Tt=1p(yt,γt,1,γt,2...γt,K∣μ,Σ,π) =∏Tt=1∏Kk=1(πkN(yt;μk,Σk))γt,k =∏Kk=1π∑Tt=1γt,kk∏Tt=1(N(yt;μk,Σk))γt,k(10)\begin{array}{l}p(y,\gamma |\mu ,\Sigma ,\pi ) = \prod\limits_{t = 1}^T {p({y_t},{\gamma _t}_{,1},{\gamma _{t,2}}...{\gamma _{t,K}}|\mu ,\Sigma ,\pi )} \\{\rm{ }} = \prod\limits_{t = 1}^T {\prod\limits_{k = 1}^K {{{({\pi _k}N({y_t};{\mu _k},{\Sigma _k}))}^{{\gamma _{t,k}}}}} } \\{\rm{ }} = \prod\limits_{k = 1}^K {\pi _k^{\sum\nolimits_{t = 1}^T {^{{\gamma _{t,k}}}} }} \prod\limits_{t = 1}^T {{{(N({y_t};{\mu _k},{\Sigma _k}))}^{{\gamma _{t,k}}}}} \end{array} \tag{10}p(y,γ∣μ,Σ,π)=t=1∏T?p(yt?,γt?,1?,γt,2?...γt,K?∣μ,Σ,π) =t=1∏T?k=1∏K?(πk?N(yt?;μk?,Σk?))γt,k? =k=1∏K?πk∑t=1T?γt,k??t=1∏T?(N(yt?;μk?,Σk?))γt,k??(10)
第1個(gè)等號(hào)到第2個(gè)等號(hào)的理解：若yt{y_t}yt?由第1類采樣而來(lái)，則有γt,1=1,γt,2=0...γt,K=0{\gamma _t}_{,1} = 1,{\gamma _{t,2}} = 0...{\gamma _{t,K}} = 0γt?,1?=1,γt,2?=0...γt,K?=0，
p(yt,γt,1,γt,2...γt,K∣μ,Σ,π)=∏Kk=1(πkN(yt;μk,Σk))γt,k =(π1N(yt;μ1,Σ1))γt,1(π2N(yt;μ2,Σ2))γt,2...(πKN(yt;μK,ΣK))γt,K =(π1N(yt;μ1,Σ1))1(π2N(yt;μ2,Σ2))0...(πKN(yt;μK,ΣK))0 =π1N(yt;μ1,Σ1)(11)\begin{array}{l}p({y_t},{\gamma _t}_{,1},{\gamma _{t,2}}...{\gamma _{t,K}}|\mu ,\Sigma ,\pi ){\rm{ = }}\prod\limits_{k = 1}^K {{{({\pi _k}N({y_t};{\mu _k},{\Sigma _k}))}^{{\gamma _{t,k}}}}} \\{\rm{ = }}{({\pi _1}N({y_t};{\mu _1},{\Sigma _1}))^{{\gamma _{t,1}}}}{({\pi _2}N({y_t};{\mu _2},{\Sigma _2}))^{{\gamma _{t,2}}}}...{({\pi _K}N({y_t};{\mu _K},{\Sigma _K}))^{{\gamma _{t,K}}}}\\{\rm{ }} = {({\pi _1}N({y_t};{\mu _1},{\Sigma _1}))^1}{({\pi _2}N({y_t};{\mu _2},{\Sigma _2}))^0}...{({\pi _K}N({y_t};{\mu _K},{\Sigma _K}))^0}\\{\rm{ }} = {\pi _1}N({y_t};{\mu _1},{\Sigma _1})\end{array} \tag{11}p(yt?,γt?,1?,γt,2?...γt,K?∣μ,Σ,π)=k=1∏K?(πk?N(yt?;μk?,Σk?))γt,k? =(π1?N(yt?;μ1?,Σ1?))γt,1?(π2?N(yt?;μ2?,Σ2?))γt,2?...(πK?N(yt?;μK?,ΣK?))γt,K? =(π1?N(yt?;μ1?,Σ1?))1(π2?N(yt?;μ2?,Σ2?))0...(πK?N(yt?;μK?,ΣK?))0 =π1?N(yt?;μ1?,Σ1?)?(11)
注意式子(11)與式子(7)的差別：如果求p(yt)p({y_t})p(yt?)則需要考慮yt{y_t}yt?有可能來(lái)源于k個(gè)類；而如果求的是p(yt,γt,1,γt,2...γt,K)p({y_t},{\gamma _t}_{,1},{\gamma _{t,2}}...{\gamma _{t,K}})p(yt?,γt?,1?,γt,2?...γt,K?)則已經(jīng)限定了yt{y_t}yt?只會(huì)來(lái)源于某個(gè)類。
第2個(gè)等式到第3個(gè)等式的理解：先交換累乘符號(hào)。由于πk{\pi _k}πk?與t無(wú)關(guān)，故而可以從內(nèi)部的累乘符號(hào)中提取出來(lái)。
實(shí)際上完全數(shù)據(jù)的似然函數(shù)描述的就是采集到這些樣本的可能性。
完全數(shù)據(jù)的對(duì)數(shù)似然函數(shù)為：
lnp(y,γ∣μ,Σ,π)=∑Kk=1(∑Tt=1γt,k)lnπk+∑Tt=1γt,k(?ln(2π)?12ln∣Σk∣?12(yt?μt)T(Σk)?1(yt?μt))(12)\ln p(y,\gamma |\mu ,\Sigma ,\pi ) = \sum\limits_{k = 1}^K {(\sum\limits_{t = 1}^T {{\gamma _{t,k}}} )\ln {\pi _k}} + \sum\limits_{t = 1}^T {{\gamma _{t,k}}} ( - \ln (2\pi ) - \frac{1}{2}\ln \left| {{\Sigma _k}} \right| - \frac{1}{2}{({y_t} - {\mu _t})^T}{({\Sigma _k})^{ - 1}}({y_t} - {\mu _t})) \tag{12}lnp(y,γ∣μ,Σ,π)=k=1∑K?(t=1∑T?γt,k?)lnπk?+t=1∑T?γt,k?(?ln(2π)?21?ln∣Σk?∣?21?(yt??μt?)T(Σk?)?1(yt??μt?))(12)
這一步應(yīng)該沒(méi)啥問(wèn)題吧。。。注意的是，此處考慮的是二維高斯分布的情況，對(duì)應(yīng)于式子(2)中的d=2。
我們的目標(biāo)就是找出一組參數(shù)(μ?,Σ?,π?)(\mu *,\Sigma *,\pi *)(μ?,Σ?,π?)使得lnp(y,γ∣μ,Σ,π)\ln p(y,\gamma |\mu ,\Sigma ,\pi )lnp(y,γ∣μ,Σ,π)最大。
那么問(wèn)題來(lái)了：lnp(y,γ∣μ,Σ,π)\ln p(y,\gamma |\mu ,\Sigma ,\pi )lnp(y,γ∣μ,Σ,π)中含有隱變量γ\gammaγ，γ\gammaγ的存在使得我們沒(méi)法最大化lnp(y,γ∣μ,Σ,π)\ln p(y,\gamma |\mu ,\Sigma ,\pi )lnp(y,γ∣μ,Σ,π)?。如果我們知道了γ\gammaγ，那么最大化lnp(y,γ∣μ,Σ,π)\ln p(y,\gamma |\mu ,\Sigma ,\pi )lnp(y,γ∣μ,Σ,π)就顯得水到渠成。
但是坑爹的就是：我們只有采樣數(shù)據(jù)yt{y_t}yt?，γ\gammaγ未知。
那么怎么辦呢？對(duì)γ\gammaγ來(lái)一個(gè)估計(jì)。
猜想我們給了一組起始參數(shù)(μ0,Σ0,π0)({\mu ^0},{\Sigma ^0},{\pi ^0})(μ0,Σ0,π0)或者優(yōu)化過(guò)的第i次迭代的參數(shù)(μi,Σi,πi)({\mu ^i},{\Sigma ^i},{\pi ^i})(μi,Σi,πi)，也就是說(shuō)每一個(gè)高斯分布的參數(shù)我們都有了，γ\gammaγ做的事不就是決定每個(gè)樣本由哪一個(gè)高斯分布產(chǎn)生的嘛，有了每個(gè)高斯分布的參數(shù)那我們就可以猜想每個(gè)樣本最有可能來(lái)源于哪個(gè)高斯分布沒(méi)錯(cuò)吧！Done!
為此我們不最大化lnp(y,γ∣μ,Σ,π)\ln p(y,\gamma |\mu ,\Sigma ,\pi )lnp(y,γ∣μ,Σ,π)（也無(wú)法最大化它），而是最大化Q函數(shù)。Q函數(shù)如下：
Q(μ,Σ,π,μi,Σi,πi)=Eγ[lnp(y,γ∣μ,Σ,π)∣Y,μi,Σi,πi]=Eγ[∑Kk=1(∑Tt=1γt,k∣yt,μi,Σi,πi)lnπk+∑Tt=1(γt,k∣yt,μi,Σi,πi)(?ln(2π)?12ln∣Σk∣?12(yt?μt)T(Σk)?1(yt?μt))]=∑Kk=1(∑Tt=1E(γt,k∣yt,μi,Σi,πi)lnπk+∑Tt=1E(γt,k∣yt,μi,Σi,πi)(?ln(2π)?12ln∣Σk∣?12(yt?μt)T(Σk)?1(yt?μt)))\begin{array}{l}Q(\mu ,\Sigma ,\pi ,{\mu ^i},{\Sigma ^i},{\pi ^i}) = {E_\gamma }[\ln p(y,\gamma |\mu ,\Sigma ,\pi )|Y,{\mu ^i},{\Sigma ^i},{\pi ^i}]\\= {E_\gamma }[\sum\limits_{k = 1}^K {(\sum\limits_{t = 1}^T {{\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i}} )\ln {\pi _k}} + \sum\limits_{t = 1}^T {({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})} ( - \ln (2\pi ) - \frac{1}{2}\ln \left| {{\Sigma _k}} \right| - \frac{1}{2}{({y_t} - {\mu _t})^T}{({\Sigma _k})^{ - 1}}({y_t} - {\mu _t}))]\\= \sum\limits_{k = 1}^K {(\sum\limits_{t = 1}^T {E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})} \ln {\pi _k}} + \sum\limits_{t = 1}^T {E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})} ( - \ln (2\pi ) - \frac{1}{2}\ln \left| {{\Sigma _k}} \right| - \frac{1}{2}{({y_t} - {\mu _t})^T}{({\Sigma _k})^{ - 1}}({y_t} - {\mu _t})))\end{array}Q(μ,Σ,π,μi,Σi,πi)=Eγ?[lnp(y,γ∣μ,Σ,π)∣Y,μi,Σi,πi]=Eγ?[k=1∑K?(t=1∑T?γt,k?∣yt?,μi,Σi,πi)lnπk?+t=1∑T?(γt,k?∣yt?,μi,Σi,πi)(?ln(2π)?21?ln∣Σk?∣?21?(yt??μt?)T(Σk?)?1(yt??μt?))]=k=1∑K?(t=1∑T?E(γt,k?∣yt?,μi,Σi,πi)lnπk?+t=1∑T?E(γt,k?∣yt?,μi,Σi,πi)(?ln(2π)?21?ln∣Σk?∣?21?(yt??μt?)T(Σk?)?1(yt??μt?)))?
其中，E(γt,k∣yt,μi,Σi,πi)E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})E(γt,k?∣yt?,μi,Σi,πi)就是對(duì)γ\gammaγ的估計(jì)：
E(γt,k∣yt,μi,Σi,πi)=p(γt,k=1∣yt,μi,Σi,πi) =p(γt,k=1,yt∣μi,Σi,πi)p(yt) =p(γt,k=1,yt∣μi,Σi,πi)∑Kk=1p(γt,k=1,yt∣μi,Σi,πi) =p(yt∣γt,k=1,μi,Σi,πi)p(γt,k=1∣μi,Σi,πi)∑Kk=1p(yt∣γt,k=1,μi,Σi,πi)p(γt,k=1∣μi,Σi,πi) =πikN(yt;μik,Σik)∑Kk=1πikN(yt;μik,Σik)(14)\begin{array}{l}E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i}){\rm{ = }}p({\gamma _{t,k}} = 1|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})\\ = \frac{{p({\gamma _{t,k}} = 1,{y_t}|{\mu ^i},{\Sigma ^i},{\pi ^i})}}{{p({y_t})}}\\ = \frac{{p({\gamma _{t,k}} = 1,{y_t}|{\mu ^i},{\Sigma ^i},{\pi ^i})}}{{\sum\nolimits_{k = 1}^K {p({\gamma _{t,k}} = 1,{y_t}|{\mu ^i},{\Sigma ^i},{\pi ^i})} }}\\ = \frac{{p({y_t}|{\gamma _{t,k}} = 1,{\mu ^i},{\Sigma ^i},{\pi ^i})p({\gamma _{t,k}} = 1|{\mu ^i},{\Sigma ^i},{\pi ^i})}}{{\sum\nolimits_{k = 1}^K {p({y_t}|{\gamma _{t,k}} = 1,{\mu ^i},{\Sigma ^i},{\pi ^i})p({\gamma _{t,k}} = 1|{\mu ^i},{\Sigma ^i},{\pi ^i})} }}\\ = \frac{{\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)}}{{\sum\nolimits_{k = 1}^K {\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)} }}\end{array} \tag{14}E(γt,k?∣yt?,μi,Σi,πi)=p(γt,k?=1∣yt?,μi,Σi,πi) =p(yt?)p(γt,k?=1,yt?∣μi,Σi,πi)? =∑k=1K?p(γt,k?=1,yt?∣μi,Σi,πi)p(γt,k?=1,yt?∣μi,Σi,πi)? =∑k=1K?p(yt?∣γt,k?=1,μi,Σi,πi)p(γt,k?=1∣μi,Σi,πi)p(yt?∣γt,k?=1,μi,Σi,πi)p(γt,k?=1∣μi,Σi,πi)? =∑k=1K?πki?N(yt?;μki?,Σki?)πki?N(yt?;μki?,Σki?)??(14)
這公式是不是很可怕？？別急，帶上幾點(diǎn)聲明，再去看公式就很好理解了！
① Q函數(shù)描述的其實(shí)就是在給定(μi,Σi,πi)({\mu ^i},{\Sigma ^i},{\pi ^i})(μi,Σi,πi)參數(shù)下，先對(duì)樣本Y做一個(gè)最有可能的劃分（每個(gè)樣本來(lái)源于各個(gè)類的可能性，即對(duì)γ\gammaγ的估計(jì)E(γt,k∣yt,μi,Σi,πi)E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})E(γt,k?∣yt?,μi,Σi,πi)），再描述能夠產(chǎn)生這組樣本的可能性（Q函數(shù)）；
② 有了對(duì)于γ\gammaγ的估計(jì)之后，Q函數(shù)只和樣本有關(guān)（傳統(tǒng)意義上的似然函數(shù)亦如此，完全數(shù)據(jù)的似然函數(shù)還與γ\gammaγ有關(guān)），而不再含有隱變量，從而使得最大化Q函數(shù)成為可能；
③ 最大化Q函數(shù)的過(guò)程實(shí)則就是使得能夠產(chǎn)生這組樣本的可能性最大，與最大化似然函數(shù)的思路如出一轍。
有個(gè)Q函數(shù)，就可以對(duì)Q函數(shù)進(jìn)行最大化，得到下一次迭代的模型參數(shù)了，即：
μi+1,Σi+1,πi+1=argmaxQ(μ,Σ,π,μi,Σi,πi)(15){\mu ^{i{\rm{ + }}1}},{\Sigma ^{i{\rm{ + }}1}},{\pi ^{i{\rm{ + }}1}}{\rm{ = }}\arg \max Q(\mu ,\Sigma ,\pi ,{\mu ^i},{\Sigma ^i},{\pi ^i}) \tag{15}μi+1,Σi+1,πi+1=argmaxQ(μ,Σ,π,μi,Σi,πi)(15)
對(duì)Q函數(shù)進(jìn)行求導(dǎo)，并另其導(dǎo)數(shù)為0，可得：
μi+1k=∑Tt=1πikN(yt;μik,Σik)∑Kk=1πikN(yt;μik,Σik)ytE(γt,k∣yt,μi,Σi,πi),k=1,2...K(16)\mu _k^{i + 1} = \frac{{\sum\nolimits_{t = 1}^T {\frac{{\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)}}{{\sum\nolimits_{k = 1}^K {\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)} }}} {y_t}}}{{E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})}},k = 1,2...K \tag{16}μki+1?=E(γt,k?∣yt?,μi,Σi,πi)∑t=1T?∑k=1K?πki?N(yt?;μki?,Σki?)πki?N(yt?;μki?,Σki?)?yt??,k=1,2...K(16)
Σi+1k=∑Tt=1πikN(yt;μik,Σik)∑Kk=1πikN(yt;μik,Σik)(yt?μik)2E(γt,k∣yt,μi,Σi,πi),k=1,2...K(17)\Sigma _k^{i + 1} = \frac{{\sum\nolimits_{t = 1}^T {\frac{{\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)}}{{\sum\nolimits_{k = 1}^K {\pi _k^iN({y_t};\mu _k^i,\Sigma _k^i)} }}} {{({y_t} - \mu _k^i)}^2}}}{{E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})}},k = 1,2...K\tag{17}Σki+1?=E(γt,k?∣yt?,μi,Σi,πi)∑t=1T?∑k=1K?πki?N(yt?;μki?,Σki?)πki?N(yt?;μki?,Σki?)?(yt??μki?)2?,k=1,2...K(17)
πi+1k=E(γt,k∣yt,μi,Σi,πi)T,k=1,2...K(18)\pi _k^{i + 1} = \frac{{E({\gamma _{t,k}}|{y_t},{\mu ^i},{\Sigma ^i},{\pi ^i})}}{T},k = 1,2...K\tag{18}πki+1?=TE(γt,k?∣yt?,μi,Σi,πi)?,k=1,2...K(18)
其中μi+1k,Σi+1k,πi+1k\mu _k^{i + 1},\Sigma _k^{i + 1},\pi _k^{i + 1}μki+1?,Σki+1?,πki+1?分別表示第（i+1）次迭代，第k個(gè)類的均值，協(xié)方差矩陣和所占的權(quán)重。
4.4 一個(gè)例子梳理EM算法的整個(gè)過(guò)程
EM算法的核心思想是：通過(guò)迭代的過(guò)程來(lái)找到一組最優(yōu)的參數(shù)(μ?,Σ?,π?)(\mu *,\Sigma *,\pi *)(μ?,Σ?,π?)，使得這組參數(shù)表示的模型最有可能產(chǎn)生現(xiàn)有的采樣數(shù)據(jù)。每次迭代的過(guò)程就是參數(shù)矯正的過(guò)程。
圖11 EM算法參數(shù)優(yōu)化過(guò)程
現(xiàn)假設(shè)初始化一組參數(shù)(μ0,Σ0,π0)({\mu ^0},{\Sigma ^0},{\pi ^0})(μ0,Σ0,π0)。在這組參數(shù)下，2類二維高斯分布如圖11綠色橢圓所示。然后利用現(xiàn)有的參數(shù)，E-step開(kāi)始對(duì)樣本數(shù)據(jù)進(jìn)行劃分（對(duì)γ\gammaγ進(jìn)行估計(jì)）。藍(lán)色的樣本大多都被劃分給第1類模型，橘黃色的樣本大多都被劃分給第2類模型。但是第1類模型還有優(yōu)化空間：第1類模型還不能使得藍(lán)色樣本出現(xiàn)的聯(lián)合概率達(dá)到最大。第2類模型也是如此。M-step便優(yōu)化了2類模型的參數(shù)，得到新的參數(shù)(μ1,Σ1,π1)({\mu ^1},{\Sigma ^1},{\pi ^1})(μ1,Σ1,π1)，使得優(yōu)化后2類高斯分布如圖11紅色橢圓所示。其中，第1類模型主要優(yōu)化的是模型均值（即橢圓的中心），第二類模型主要優(yōu)化的是模型協(xié)方差矩陣（即橢圓的長(zhǎng)軸、短軸和長(zhǎng)短軸的方向）。然后重復(fù)進(jìn)行E-step和M-step，直到參數(shù)(μ,Σ,π)(\mu ,\Sigma ,\pi )(μ,Σ,π)收斂。
最后談?wù)劵旌细咚鼓Ｐ偷膮?shù)π\(zhòng)piπ。
混合高斯模型的參數(shù)μ,Σ\mu ,\Sigmaμ,Σ比較好理解，用于描述各個(gè)高斯分布的形狀，對(duì)于它們的調(diào)整也比較直觀：使得本高斯分布能夠更好地接納被劃分到這類分布的樣本。而為什么要有參數(shù)π\(zhòng)piπ？它描述的是各個(gè)高斯分布所占的比重，如果不加“歧視”的話（樣本來(lái)源于各個(gè)高斯分布的可能性一致），則有πk=1/K{\pi _k} = 1/Kπk?=1/K；而如果對(duì)于某一類高斯分布（即為i）有側(cè)重的話，則相應(yīng)的πi{\pi _i}πi?較大，體現(xiàn)在圖11中就是被分配給各個(gè)類的樣本數(shù)占樣本總數(shù)的比例。如果一輪優(yōu)化后，某一類高斯分布又接納了更多樣本，則其πi{\pi _i}πi?變大，反之變?。ㄋ詧D11從綠色橢圓調(diào)整為紅色橢圓實(shí)際上兩個(gè)類所對(duì)應(yīng)的權(quán)重也被優(yōu)化了）。
而從本質(zhì)上來(lái)看參數(shù)π\(zhòng)piπ，則是為了混合高斯模型能有更好的曲面擬合能力。當(dāng)參數(shù)π\(zhòng)piπ退化為某一類高斯分布的權(quán)重遠(yuǎn)遠(yuǎn)大于其他類高斯分布的時(shí)候，混合高斯模型就退化成了單高斯模型！
圖12和圖13梳理了高斯分布和混合高斯分布參數(shù)估計(jì)的邏輯流程。
圖12 高斯分布參數(shù)估計(jì)邏輯流程
圖13 混合高斯分布參數(shù)估計(jì)邏輯流程
相對(duì)比于高斯分布的參數(shù)估計(jì)，混合高斯分布的參數(shù)估計(jì)更加復(fù)雜。主要原因在于隱變量的存在。而為什么混合高斯分布的參數(shù)估計(jì)需要多次迭代循環(huán)進(jìn)行？是因?yàn)镋M算法中對(duì)于γ\gammaγ的估計(jì)利用的是初始化或者第i步迭代的參數(shù)(μi,Σi,πi)({\mu ^i},{\Sigma ^i},{\pi ^i})(μi,Σi,πi)，這對(duì)于樣本的分類劃分是有誤差的。所以它只能通過(guò)多次迭代優(yōu)化尋找更佳的參數(shù)來(lái)抵消這一誤差。
終于把這篇文章梳理完了。世界杯要結(jié)束了，偽球迷也想見(jiàn)證一下冠軍誕生。至此，本文結(jié)束。
