EM演算法與GMM（高斯混合聚類）

EM Algorithm

EM(Expectation maximization)演算法，也即期望最大化演算法，作為「隱變數」（屬性變數不可知）估計的利器在自然語言處理（如HMM中的Baum-Welch演算法）、高斯混合聚類、心理學、定量遺傳學等含有隱變數的概率模型參數極大似然估計中有著十分廣泛的應用。EM演算法於1977年由Arthur Dempster, Nan Laird和Donald Rubin總結提出，其主要通過E步（exceptation），M步（maximization）反覆迭代直至似然函數收斂至局部最優解。由於其方法簡潔、操作有效，EM演算法曾入選「數據挖掘十大演算法」，可謂是機器學習經典演算法之一。

Introduction

EM演算法推導一

對於概率模型，當模型中的變數均為觀測變數時，我們可以直接使用給定數據通過最大似然估計（頻率學派）或貝葉斯估計（貝葉斯學派）這兩種方法求解。然而當我們的模型中存在隱變數時，我們將無法使用最大似然估計直接求解，這時即導出EM演算法。

假設一個概率模型中同時存在隱變數和可觀測變數，我們學習的目標是極大化觀測變數關於模型參數的對數似然，即：

$L( heta)=logP(Y| heta)=logsum_ZP(Y,Z| heta)=log(sum_zP(Y|Z, heta)P(Z| heta)) ag{1}$

式（1）中我們假設直接優化是很困難的，但是優化完整數據的似然函數相對容易，同時利用概率乘法公式將展開。然而由於未觀測變數的存在，上式仍求解困難，因此我們通過迭代逐步最大對數似然，這裡假設第次迭代後的估計值為。根據要求，我們希望新估計的參數使增加，即，且逐步使達到最大，因此考慮兩者之差：

$L( heta)-L( heta^i)=log(sum_zP(Y|Z, heta)P(Z| heta))-log(P(Y| heta^i))\ =log(sum_ZP(Z|Y, heta^i)frac{P(Y|Z, heta)P(Z| heta)}{P(Z|Y, heta^i)})-log(P(Y| heta^i)) ag{2}$

這裡我們根據Jensen（琴生）不等式：，其中有：

$log(sum_ZP(Z|Y, heta^i)frac{P(Y|Z, heta)P(Z| heta)}{P(Z|Y, heta^i)})-log(P(Y| heta^i))geqsum_ZP(Z|Y, heta^i)log(frac{P(Y|Z, heta)P(Z| heta)}{P(Z|Y, heta^i)})-log(P(Y| heta^i) ag{3}$

同時由於，式（3）可進一步寫為：

$sum_ZP(Z|Y, heta^i)log(frac{P(Y|Z, heta)P(Z| heta)}{P(Z|Y, heta^i)})-log(P(Y| heta^i)=sum_ZP(Z|Y, heta^i)log(frac{P(Y|Z, heta)P(Z| heta)}{P(Z|Y, heta^i)P(Y| heta^i)}) ag{4}$

因此有：

$L( heta)geq L( heta^i)+sum_ZP(Z|Y, heta^i)log(frac{P(Y|Z, heta)P(Z| heta)}{P(Z|Y, heta^i)P(Y| heta^i)})=B( heta, heta^i) ag{5}$

因此，即為的下界。故當增大時也將同時增加，為使取得最大，則我們必須在次迭代時選擇的 $heta^{i+1}$ 為使第次迭代取得最大的，即：

$heta^{i+1}=arg;max_{ heta}B( heta, heta^i)=arg;max_{ heta}(L( heta^i)+sum_ZP(Z|Y, heta^i)log(frac{P(Y|Z, heta)P(Z| heta)}{P(Z|Y, heta^i)P(Y| heta^i)}))\ =arg;max_{ heta}(sum_ZP(Z|Y, heta^i)log(P(Y|Z, heta)P(Z| heta)))\ =arg;max_{ heta}(sum_ZP(Z|Y, heta^i)logP(Y,Z| heta)=arg;max_{ heta}Q( heta, heta^i) ag{6}$

在上式的求解中我們略去了對求解極大化而言的常數項和。

因此在EM演算法的每一迭代中，我們均需求解使得取得最大值的，使得下一不迭代的 $heta^{i+1}= heta$ ，這樣如此反覆提高最大似然的下界，直至逼近的最優解（最大值）。

EM演算法推導二

這裡我們採用變分的方法，假設隱變數服從任一分布為，則。故對於同樣有：

$L( heta)=logP(Y| heta)=sum_Z q(Z)logP(Y| heta)=sum_Zq(Z)logfrac{P(Z|Y, heta)P(Y| heta)}{P(Z|Y, heta)}=sum_Zq(Z)logfrac{P(Y,Z| heta)}{P(Z|Y, heta)}\ =sum_Zq(Z)log(frac{P(Y,Z| heta)}{q(Z)}frac{q(Z)}{P(Z|Y, heta)})=sum_Zq(Z)(logfrac{P(Y,Z| heta)}{q(Z)}-logfrac{q(Z)}{P(Z|Y, heta)})\ =underbrace{sum_Zq(Z)logfrac{P(Y,Z| heta)}{q(Z)}}_{(1)}-underbrace{sum_Zq(Z)logfrac{P(Z|Y, heta)}{q(Z)}}_{(2)} ag{7}$

記（1）為 $L(q, heta)=sum_Zq(Z)logfrac{P(Y,Z| heta)}{q(Z)}$ ，（2）為 $KL(q||p)=-sum_Zq(Z)logfrac{P(Z|Y, heta)}{q(Z)}$ 。其中即為KL散度（相對熵），主要反映變數分布的相似性，可以看出KL散度=交叉熵-信息熵，故交叉熵在某種意義上與KL散度等價。有：

由於，因此即為對數似然函數的下界。同理在每一次迭代中我們均需要最大化下界，則在第次迭代中即有：

$q(Z)=P(Z|Y, heta^i) heta^{i+1}=arg;max_{ heta}L(q, heta)=sum_Zq(Z)logfrac{P(Y,Z| heta)}{q(Z)} =arg;max_{ heta}sum_ZP(Z|Y, heta^i)logfrac{P(Y,Z| heta)}{P(Z|Y, heta^i)}\ =arg;max_{ heta}(sum_ZP(Z|Y, heta^i)(logP(Y,Z| heta)-logP(Z|Y, heta^i))=arg;max_{ heta}Q( heta, heta^i)+const ag{9}$