摘自《統計學習方法》 李航著 清華大學出版社

連結：

EM介紹請點選這裡

EM算法在高斯混合模型中的應用

EM算法的一個重要應用是高斯混合模型的參數估計。高斯混合模型應用廣泛，在許多情況下，EM算法是學習高斯混合模型的有效方法。想詳細了解單高斯模型和混合高斯模型請點選這裡

高斯混合模型

定義

高斯混合模型是指具有如下形式的機率分布模型：

P(y|Θ)=∑Kk=1αϕ(y|Θk) 公式（1）

ϕ(y|Θk)=12π√σkexp−(y−μk)22σ2k 公式（2）

其中， α 是系數， α≥0 ， ∑Kk=1αk=1 ； ϕ(y|Θk) 是高斯分布密度， Θ=(μk,σk) 。公式（2）我們稱為第K個分模型。

一般混合模型可以由任意機率分布密度代替公式（2）中的高斯分布密度，我們隻介紹最常用的高斯混合模型。

高斯混合模型參數估計的EM算法

假設觀測資料 y1,y2,...,yN 由高斯混合模型生成，

P(y|Θ)=∑Kk=1ϕ(y|Θk) 公式（3）

其中， Θ=(α1,α2,...,αk;θ1,θ2,...,θk) 。我們用EM算法估計高斯混合模型的參數 θ 。

明确隐變量，寫出完全資料的對數似然函數

可以設想觀測資料 yi ， j=1,2,...,N ，是這樣産生的：首先依機率 αk 選擇第 k 個高斯分布分模型ϕ(y|θk)；然後依第 k 個分模型的機率分布ϕ(y|θk)生成觀測資料 yj 。這是觀測資料 yj ， j=1,2,...,N ，是已隻的；反映觀測資料 yj 來自第 k 個分模型的資料是未知的，k=1,2,...,K，以隐變量 γjk 表示，其定義如下：

γk={1第j個觀察來自第k個分模型0 否則

j=1,2,...,N;k=1,2,....,K 公式（4）

其中 γjk 是0-1随機變量。

有了觀測資料 yj 及未觀測資料 γjk ，那麼完全資料是

(yj,γj1,γj2,...,γjK)，j=1,2,...,N

于是可以寫出完全資料的似然函數：

P(y,γ|θ)=∏j=1NP(yj,γj1,γj2,...,γjk|θ)

=∏Kk=1∏Nj=1[αkϕ(yj|θk)]γjk 

=∏Kk=1αnk∏Nj=1[ϕ(yj|θk)]γjk  

=∏αnkk∏Nj=1[12π√σkexp(−(yj−μk)22σ2k)]γjk  

式中， nk=∑Nj=1γjk ， ∑Kk=1nk=N

公式注釋：

上述第二個等号為什麼成立？我們思考一下。在上節EM算法介紹中，我們在硬币例題介紹了

P(y|θ)=∑zP(y,z|θ)=∑zP(z|θ)P(y|z,θ)

現在，我們來對應一下。 P(z|θ) 對應高斯混合模型的 P(γjk|θ) ，也就是 αk 。在第 k 個模型的情況下，觀測值yi屬于第 k 個模型的機率是ϕ(yj|θk)，也就是 ϕ(y|Θk)=12π√σkexp−(y−μk)22σ2k ，然而不隻是由k這種模型，還有1~K個模型，是以是每個模型算出機率後相乘。注意 γjk 的值不是1就是0，是以當 yi 不屬于第 k 類的時候，任何數的0次方為1。當當yi屬于第 k 類的時候，γjk的值為1，任何數的1次方還是其本身。

那麼，完全資料的對數似然函數為:

logP(y,γ|θ)=∑Kk=1nk{logαk+∑Nj=1γjk[log(12π√)−logσk−12σ2k(yj−μk)2]}

EM算法的E步：确定Q函數

Q(θ,θ(i))=E(logP(y,γ|θ)|y,θ(i))

=E{∑Kk=1{nklogαk+∑Nj=1γjk[log12π√−logσk−12σ2k(yj−μk)2]}}

=∑Kk=1{∑Nj=1(Eγjk)logαk+∑Nj=1(Eγrk)[log(1(√2π))−logσk−12σ2(yj−μk)2]}

公式（4）

從第2個等号到第3個等号這裡涉及到一些期望的知識，以圖檔的形式補充到下面。

資料來源：http://www.360doc.com/content/13/1124/03/9482_331690142.shtml

這裡需要計算 E(γjk|y,θ) ，記為 γ‘jk 。

γjk‘=E(γjk|y,θ)=P(γrk=1|y,θ)

=P(γ=1,yj|θ)∑Kk=1P(γjk=1,yj|θ) 

=P(yj|γjk=1,θ)P(γjk=1|θ)∑Kk=1P(yj|γjk=1,θ)P(γjk=1|θ)  

=αkϕ(yj|θk)∑Kk=1αkϕ(yj|θk) j=1, 2,…, N; k = 1, 2, …, K

γ‘jk 是在目前模型參數下第j個觀測資料來自第 k 個分模型的機率，稱為分模型k對觀測資料 yj 的響應度。

将 γ‘jk=Eγjk 及 nk=∑Nj=1Eγjk 代入公式（4）（确定Q函數的公式），即可獲得

Q(θ,θ(i))=∑Kk=1{nklogαk+∑Nk=1γ‘jk[log(12π√)−logσk−12σ2k(yj−μk)2]} 公式（5）

确定EM算法的M步

疊代的M步是求函數 Q(θ,θ(i)) 對 θ 的極大值，即求新一輪疊代的模型參數：

θ(i+1)=argmaxQ(θ,θ(i))

用 μk^,σ2k^ 及 αk^ ， k=1,2,...,K ，表示 θ(i+1) 的各個參數。求 μk^ ， σ2k^ 隻需将公式（5）分别對 μk^ ， σ2k^ 求偏導數并令其為0，即可得到；求 αk^ 是在 ∑Kk=1αk=1 條件下求偏導并令其為0得到的。其結果如下所示：

μk^=∑Nj=1γjkyj^∑Nj=1γjk^k=1,2,...,K 公式（6）

σ2k^=∑Nj=1γjk^(yi−uk)2∑Nj=1γjk^k=1,2,...,K    公式（7）

αk^=nkN=∑Nj=1γjk^Nk=1,2,...,K    公式（7）

重複以上計算，直到對數似然函數值不再有明顯變化為止。

算法總結

輸入：觀測資料 y1,y2,...,yN ，高斯混合模型；

輸出：高斯混合模型參數；

步驟：

（1）取參數的初始值開始疊代

（2）E步：依據目前模型參數，計算分模型k對觀測資料 yi 的響應度

γk^=αkϕ(yj|θk)∑Kk=1αkϕ(yj|θk)j=1,2,...,N; k=1,2,...,K

（3）M步：計算新一輪疊代的模型參數:

μk^=∑Nj=1γjkyj^∑Nj=1γjk^k=1,2,...,K

σ2k^=∑Nj=1γjk^(yi−uk)2∑Nj=1γjk^k=1,2,...,K   

αk^=nkN=∑Nj=1γjk^N k=1,2,...,K 

（4）重複第（2）步和第（3）步，直到收斂。