EM演算法與混合高斯

極大似然估計

考慮一個高斯分布 p(mathbf{x}mid{	heta}) ,其中 	heta=(mu,Sigma) 。樣本集 X={x_1,...,x_N} 中每個樣本都是獨立的從該高斯分布中抽取得到的,滿足獨立同分布假設。

因此,取到這個樣本集的概率為:

egin{aligned} p(Xmid{	heta}) &= prod_{i=1}^Np(x_imid	heta) end{aligned}

我們要估計模型參數向量 	heta 的值,由於現在樣本集 X 已知,所以可以把 p(Xmid{	heta}) 看作是參數向量 	heta 的函數,稱之為樣本集 X 下的似然函數 L(	hetamid{X})

egin{aligned} L(	hetamid{X}) &= prod_{i=1}^Np(x_imid	heta) end{aligned}

極大似然估計要做的就是最大化該似然函數,找到能夠使得 p(Xmid{	heta}) 最大的參數向量 	heta ,換句話說,就是找到最符合當前觀測樣本集 X 的模型的參數向量 	heta

為方便運算,通常我們使用似然函數的對數形式(可以將乘積運算轉化為求和形式),稱之為「對數似然函數」。由於底數大於 1 的對數函數總是單調遞增的,所以使對數似然函數達到

最大值的參數向量也能使似然函數本身達到最大值。故,對數似然函數為:

egin{aligned} L(	hetamid{X}) &= logBig(prod_{i=1}^Np(x_imid	heta)Big) \ &= sum_{i=1}^Nlog{p(x_imid	heta)} end{aligned}

參數的估計值 hat{	heta}=argunderset{	heta}{max}L(	hetamid{X})

混合高斯模型及其求解困境

混合高斯模型是指由 k 個高斯模型疊加形成的一種概率分布模型,形式化表示為:

p(mathbf{x}mid{Theta}) = sum_{l=1}^{k}alpha_lp_l(mathbf{x}mid{	heta_l})

其中,參數  Theta=(alpha_1,...,alpha_k,	heta_1,...,	heta_k) ,並且有 Sigma_{l=1}^{k}alpha_l=1alpha_l 代表單個高斯分布在混合高斯模型中的權重。

現在我們假設觀測樣本集 X=(x_1,...x_N) 來自於該混合高斯模型,滿足獨立同分布假設。為了估計出該混合高斯模型的參數 Theta ,我們寫出這 n 個數據的對數似然函數:

egin{aligned} L(Thetamid{X}) &= logBig(p(Xmid{Theta})Big) \ &= logBig(prod_{i=1}^{N}p(x_imid{Theta})Big) \ &= sum_{i=1}^{N}logBig(p(x_imid{Theta})Big) \ &= sum_{i=1}^{N}logBig(sum_{l=1}^{k}alpha_lp_l(x_imid{	heta_l})Big) end{aligned}

我們的目標就是要通過最大化該似然函數從而估計出混合高斯模型的參數  Theta

觀察該式,由於對數函數中還包含求和式,因此如果仿照極大似然估計單個高斯分布參數的方法來求解這個問題,是無法得到解析解的。所以我們要尋求更好的方式來解決這個問題。

EM演算法基本過程

考慮一個樣本集 X 是從某種分布(參數為 Theta )中觀測得到的,我們稱之為不完全數據(incomplete data)。我們引入一個無法直接觀測得到的隨機變數集合 Z ,叫作隱變數。 XZ 連在一起稱作完全數據。我們可以得到它們的聯合概率分布為:

p(X,Zmid{Theta})=p(Zmid{X,Theta})p(Xmid{Theta})

我們定義一個新的似然函數叫作完全數據似然:

L(Thetamid{X,Z})=p(X,Zmid{Theta})

我們可以通過極大化這樣一個似然函數來估計參數 Theta 。然而Z是隱變數,其分布未知上式無法直接求解,我們通過計算完全數據的對數似然函數關於 Z 的期望,來最大化已觀測數據的邊際似然。我們定義這樣一個期望值為 Q 函數:

egin{aligned} Q(Theta,Theta^{g}) &= mathbb{E}_ZBig[log{p(X,ZmidTheta)}mid{X,Theta^g}Big]\ &= sum_Zlog{p(X,Zmid{Theta})p(Zmid{X,Theta^g})} end{aligned}

其中, Theta^g 表示當前的參數估計值, X 是觀測數據,都作為常量。 Theta 是我們要極大化的參數, Z 是來自於分布 p(Zmid{X,Theta^g}) 的隨機變數。 p(Zmid{X,Theta^g}) 是未觀測變數的邊緣分布,並且依賴於觀測數據 X 和當前模型參數 Theta^g

EM 演算法中的 E-setp 就是指上面計算期望值的過程。

M-step 則是極大化這樣一個期望,從而得到能夠最大化期望的參數 Theta

egin{aligned} Theta^{g+1} &= argunderset{	heta}{max}Q(Theta,Theta^{g}) \ &= argunderset{	heta}{max}sum_Zlog{p(X,Zmid{Theta})p(Zmid{X,Theta^g})} end{aligned}

重複執行 E-stepM-step 直至收斂。

收斂性證明

(略)

EM演算法應用於高斯混合模型

回到高斯混合模型的參數估計問題,我們將樣本集 X 看作不完全數據,考慮一個無法觀測到的隨機變數 Z={z_i}_{i=1}^{N} ,其中 z_iin{1,...,k} ,用來指示每一個數據點是由哪一個高斯分布產生的,(可以類比成一個類別標籤,這個標籤我們是無法直接觀測得到的)。比如, z_i=k 表示第 i 個樣本點是由混合高斯模型中的第 k 個分量模型生成的。

那麼完全數據 (X,Z) 的概率分布為:

egin{aligned} p(X,Zmid{Theta}) &= prod_{i=1}^{N}p(x_i,z_imid{Theta}) \ &= prod_{i=1}^{N}p(z_imid{Theta})p(x_imid{z_i,Theta}) \ &= prod_{i=1}^{N}alpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}}) end{aligned}

在混合高斯模型問題中, p(z_imidTheta) 是指第 z_i 個模型的先驗概率,也就是其權重 alpha_{z_i} 。給定樣本點來源的高斯分布類別後, p(x_imid{z_i,Theta}) 可以寫成對應的高斯分布下的概率密度形式,即 p_{z_i}(x_imid{	heta_{z_i}})

根據貝葉斯公式,又可以得到隱變數的條件概率分布:

egin{aligned} p(Zmid{X,Theta}) &= prod_{i=1}^Np(z_imid{x_i,Theta}) \ &= prod_{i=1}^Nfrac{p(x_i,z_i,Theta)}{p(x_i,Theta)} \ &= prod_{i=1}^Nfrac{alpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})}{sum_{z_i=1}^kalpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})} end{aligned}

因此,我們可以寫出相應的 Q 函數:

egin{aligned} Q(Theta,Theta^{g}) &= sum_Zlog{p(X,Zmid{Theta})p(Zmid{X,Theta^g})} \ &= sum_Zlog{prod_{i=1}^{N}alpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})}prod_{i=1}^Np(z_imid{x_i,Theta^g}) \ &= sum_Zsum_{i=1}^{N}logBig(alpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})Big)prod_{i=1}^Np(z_imid{x_i,Theta^g}) \ &= sum_{z_1=1}^ksum_{z_2=1}^k...sum_{z_N=1}^ksum_{i=1}^{N}logBig(alpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})Big)prod_{i=1}^Np(z_imid{x_i,Theta^g}) \ &= sum_{z_1=1}^ksum_{z_2=1}^k...sum_{z_N=1}^klogBig(alpha_{z_1}p_{z_1}(x_1mid{	heta_{z_1}})Big)prod_{i=1}^Np(z_imid{x_i,Theta^g}) \ &quad+ sum_{z_1=1}^ksum_{z_2=1}^k...sum_{z_N=1}^ksum_{i=2}^{N}logBig(alpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})Big)prod_{i=1}^Np(z_imid{x_i,Theta^g}) \ &= mathcal{A} + mathcal{B} end{aligned}

其中,

egin{aligned} mathcal{A} &= sum_{z_1=1}^ksum_{z_2=1}^k...sum_{z_N=1}^klogBig(alpha_{z_1}p_{z_1}(x_1mid{	heta_{z_1}})Big)prod_{i=1}^Np(z_imid{x_i,Theta^g}) \ &= sum_{z_1=1}^klogBig(alpha_{z_1}p_{z_1}(x_1mid{	heta_{z_1}})Big)p(z_1mid{x_1,Theta^g})underbrace{sum_{z_2=1}^k...sum_{z_N=1}^kprod_{i=2}^Np(z_imid{x_i,Theta^g})}_{result=1} \ &= sum_{z_1=1}^klogBig(alpha_{z_1}p_{z_1}(x_1mid{	heta_{z_1}})Big)p(z_1mid{x_1,Theta^g}) \ end{aligned}

mathcal{B} 式可以按照同樣的操作技巧進行分解,故不贅述。

並且,我們用 l 代替 z_i 來簡化我們的表達式。

因此, Q 函數可以化簡為:

egin{aligned} Q(Theta,Theta^{g}) &= sum_{i=1}^Nsum_{z_i=1}^klogBig(alpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})Big)p(z_imid{x_i,Theta^g}) \ &= sum_{i=1}^Nsum_{l=1}^klogBig(alpha_{l}p_{l}(x_imid{	heta_{l}})Big)p(lmid{x_i,Theta^g}) \ &= sum_{l=1}^ksum_{i=1}^NlogBig(alpha_{l}p_{l}(x_imid{	heta_{l}})Big)p(lmid{x_i,Theta^g}) \ &= sum_{l=1}^ksum_{i=1}^Nlog(alpha_l)p(lmid{x_i,Theta^g})+sum_{l=1}^ksum_{i=1}^NlogBig(p_l(x_imid{	heta_{l}})Big)p(lmid{x_i,Theta^g}) \ end{aligned}

這樣,我們可以對包含參數 alpha_l 和參數 	heta_l 的項分別進行最大化從而得到各自的估計值。

估計參數 alpha_l

這個問題可以表示為下面的約束最優化問題:

egin{aligned} underset{alpha_l}{max}&quad sum_{l=1}^ksum_{i=1}^Nlog(alpha_l)p(lmid{x_i,Theta^g}) \ s.t.& quadsum_{l=1}^k{alpha_l}=1 end{aligned}

引入拉格朗日乘子 lambda ,構建拉格朗日函數:

egin{aligned} mathcal{L}(alpha_1,...,alpha_2,lambda) &= sum_{l=1}^ksum_{i=1}^Nlog(alpha_l)p(lmid{x_i,Theta^g})-lambdaBig(sum_{l=1}^k{alpha_l}-1Big) \ &= sum_{l=1}^klog(alpha_l)sum_{i=1}^Np(lmid{x_i,Theta^g})-lambdaBig(sum_{l=1}^k{alpha_l}-1Big) end{aligned}

對參數 alpha_l 求偏導並令其為 0

frac{partialmathcal{L}}{partialalpha_l}=frac{1}{alpha_l}sum_{i=1}^{N}p(lmid{x_i,Theta^g})-lambda=0

得到:

alpha_l=frac{1}{lambda}sum_{i=1}^{N}p(lmid{x_i,Theta^g})

代回約束條件,有:

1-frac{1}{lambda}sum_{l=1}^{k}sum_{i=1}^{N}p(lmid{x_i,Theta^g})=0

在之前的推導中,我們得到過這樣一個公式,即隱變數 z 的條件分布:

egin{aligned} p(z_imid{x_i,Theta}) &= frac{p(x_i,z_i,Theta)}{p(x_i,Theta)} \ &= frac{alpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})}{sum_{z_i=1}^kalpha_{z_i}p_{z_i}(x_imid{	heta_{z_i}})} end{aligned}

同樣用 l 代替 z_i 來簡化我們的表達式,得到:

egin{aligned} p(lmid{x_i,Theta}) &= frac{p(x_i,l,Theta)}{p(x_i,Theta)} \ &= frac{alpha_{l}p_{l}(x_imid{	heta_{l}})}{sum_{l=1}^kalpha_{l}p_{l}(x_imid{	heta_{l}})} end{aligned}

故將其代回之前的式子,得到:

egin{aligned} 1-frac{1}{lambda}sum_{l=1}^{k}sum_{i=1}^{N}frac{alpha_{l}p_{l}(x_imid{	heta_{l}})}{sum_{l=1}^kalpha_{l}p_{l}(x_imid{	heta_{l}})}&=0 \ 1-frac{1}{lambda}sum_{i=1}^{N}sum_{l=1}^{k}frac{alpha_{l}p_{l}(x_imid{	heta_{l}})}{sum_{l=1}^kalpha_{l}p_{l}(x_imid{	heta_{l}})}&=0 \ 1-frac{1}{lambda}sum_{i=1}^{N}(1)&=0 \ 1-frac{N}{lambda}&=0 \ lambda&=N end{aligned}

故,

alpha_l=frac{1}{N}sum_{i=1}^{N}p(lmid{x_i,Theta^g})

在估計剩下的參數之前,先補充一下一會要用到的線性代數知識。

【知識點~Matrix Algebra】

矩陣的跡等於矩陣的主對角線上元素之和,且有以下性質

  • tr(A+B)=tr(A)+tr(B)
  • tr(AB)=tr(BA)
  • sum_ix_i^TAx_i=tr(AB) quad where quad B=sum_ix_ix_i^T

|A| 表示矩陣 A 的行列式,當 A 可逆時,有以下性質

  • |A^{-1}|=|A|^{-1}

a_{i,j} 表示矩陣 A 的第 i 行, j 列的元素,給出矩陣函數求導的一些公式:

  • frac{partial{x^TAx}}{partial{x}}=(A+A^T)x
  • frac{partiallog|A|}{partial{A}}=2A^{-1}-diag(A^{-1})
  • frac{partial{tr(AB)}}{partial{A}}=B+B^T-diag(B)

估計參數 	heta_l

對於 d 維高斯分布來說,參數 	heta=(mu,Sigma) ,且:

p_l(xmid{mu_l,Sigma_l})=frac{1}{(2pi)^{d/2}vertSigma_lvert^{1/2}}expig[-frac{1}{2}(x-mu_l)^TSigma_l^{-1}(x-mu_l)ig]

估計參數 mu_l

egin{aligned} &quadsum_{l=1}^ksum_{i=1}^NlogBig(p_l(x_imid{mu_l,Sigma_l})Big)p(lmid{x_i,Theta^g}) \ &= sum_{l=1}^ksum_{i=1}^NBig(log(2pi^{-d/2})+log|Sigma_l|^{-1/2}-frac{1}{2}(x_i-mu_l)^TSigma_l^{-1}(x_i-mu_l)Big)p(lmid{x_i,Theta^g}) \ end{aligned}

忽略其中的常數項(因為常數項在求導之後為 0 ),上式化簡為:

sum_{l=1}^ksum_{i=1}^NBig(-frac{1}{2}log|Sigma_l|-frac{1}{2}(x_i-mu_l)^TSigma_l^{-1}(x_i-mu_l)Big)p(lmid{x_i,Theta^g})

關於 mu_l 求導,得到:

sum_{l=1}^ksum_{i=1}^NBigg(-frac{1}{2}Big(Sigma_l^{-1}+(Sigma_l^{-1})^TBig)Big(x_i-mu_lBig)Big(-1Big)Bigg)p(lmid{x_i,Theta^g})

因為協方差矩陣 Sigma 為對稱陣,故 frac{1}{2}Big(Sigma_l^{-1}+(Sigma_l^{-1})^TBig)=Sigma_l^{-1} ,因此,上式繼續化簡為:

sum_{l=1}^ksum_{i=1}^NSigma_l^{-1}(x_i-mu_l)p(lmid{x_i,Theta^g})

令該式為 0 ,得到:

egin{aligned} sum_{l=1}^kSigma_l^{-1}sum_{i=1}^N{mu_l}p(lmid{x_i,Theta^g})&=sum_{l=1}^kSigma_l^{-1}sum_{i=1}^N{x_i}p(lmid{x_i,Theta^g}) \ sum_{i=1}^N{mu_l}p(lmid{x_i,Theta^g})&=sum_{i=1}^N{x_i}p(lmid{x_i,Theta^g}) \ mu_lsum_{i=1}^Np(lmid{x_i,Theta^g})&=sum_{i=1}^N{x_i}p(lmid{x_i,Theta^g}) \ mu_l&=frac{sum_{i=1}^N{x_i}p(lmid{x_i,Theta^g})}{sum_{i=1}^Np(lmid{x_i,Theta^g})} \ end{aligned}

估計參數 Sigma_l

忽略其中的常數項(因為常數項在求導之後為 0 ),上式化簡為:

egin{aligned} &quadsum_{l=1}^ksum_{i=1}^NlogBig(p_l(x_imid{mu_l,Sigma_l})Big)p(lmid{x_i,Theta^g}) \ &= sum_{l=1}^ksum_{i=1}^NBig(log(2pi^{-d/2})+log|Sigma_l|^{-1/2}-frac{1}{2}(x_i-mu_l)^TSigma_l^{-1}(x_i-mu_l)Big)p(lmid{x_i,Theta^g}) \ end{aligned}egin{aligned} &quadsum_{l=1}^ksum_{i=1}^NBig(-frac{1}{2}log|Sigma_l|-frac{1}{2}(x_i-mu_l)^TSigma_l^{-1}(x_i-mu_l)Big)p(lmid{x_i,Theta^g}) \ &= sum_{l=1}^ksum_{i=1}^NBig(frac{1}{2}log|Sigma_l^{-1}|p(lmid{x_i,Theta^g})-frac{1}{2}(x_i-mu_l)^TSigma_l^{-1}(x_i-mu_l)p(lmid{x_i,Theta^g})Big) \ &= sum_{l=1}^kBig(frac{1}{2}log|Sigma_l^{-1}|sum_{i=1}^{N}p(lmid{x_i,Theta^g})-frac{1}{2}Sigma_l^{-1}sum_{i=1}^{N}(x_i-mu_l)^T(I)(x_i-mu_l)p(lmid{x_i,Theta^g})Big) \ &= sum_{l=1}^kBigg(frac{1}{2}log|Sigma_l^{-1}|sum_{i=1}^{N}p(lmid{x_i,Theta^g})-frac{1}{2}trBig(Sigma_l^{-1}sum_{i=1}^{N}(x_i-mu_l)(x_i-mu_l)^Tp(lmid{x_i,Theta^g})Big)Bigg) \ end{aligned}

考慮方程 S(mu_l,Sigma_{l}^{-1}) 為上式中 sum_{l=1}^k 內部的式子,即:

S(mu_l,Sigma_{l}^{-1})= frac{1}{2}log|Sigma_l^{-1}|sum_{i=1}^{N}p(lmid{x_i,Theta^g})-frac{1}{2}trBig(Sigma_l^{-1}underbrace{sum_{i=1}^{N}(x_i-mu_l)(x_i-mu_l)^Tp(lmid{x_i,Theta^g})}_{M_l}Big)

S 關於 Sigma_l^{-1} 求導數,得到:

egin{aligned} frac{partial{S(mu_l,Sigma_l^{-1})}}{partial{Sigma_l^{-1}}} &= frac{1}{2}sum_{i=1}^{N}p(lmid{x_i,Theta^g})Big(2Sigma_l-diag(Sigma_l)Big)-frac{1}{2}frac{partial{tr(Sigma_l^{-1}M_l)}}{partialSigma_l^{-1}} \ &= frac{1}{2}sum_{i=1}^{N}p(lmid{x_i,Theta^g})Big(2Sigma_l-diag(Sigma_l)Big)-frac{1}{2}Big(2M_l-diag(M_l)Big) \ end{aligned}

令該導數值為 0 ,得到:

egin{aligned} sum_{i=1}^{N}p(lmid{x_i,Theta^g})Big(2Sigma_l-diag(Sigma_l)Big) &= Big(2M_l-diag(M_l)Big) \ 2Big(underbrace{sum_{i=1}^{N}Sigma_lp(lmid{x_i,Theta^g})-M_l}_{mathcal{A}}Big) &= diagBig(underbrace{sum_{i=1}^{N}Sigma_lp(lmid{x_i,Theta^g})-M_l}_{mathcal{A}}Big) \ 2mathcal{A} &= diag(mathcal{A}) end{aligned}

解得 mathcal{A}=0 ,即:

egin{aligned} sum_{i=1}^{N}Sigma_lp(lmid{x_i,Theta^g}) &= M_l \ sum_{i=1}^{N}Sigma_lp(lmid{x_i,Theta^g}) &= sum_{i=1}^{N}(x_i-mu_l)(x_i-mu_l)^Tp(lmid{x_i,Theta^g}) \ Sigma_lsum_{i=1}^{N}p(lmid{x_i,Theta^g}) &= sum_{i=1}^{N}(x_i-mu_l)(x_i-mu_l)^Tp(lmid{x_i,Theta^g}) \ Sigma_l &= frac{sum_{i=1}^{N}(x_i-mu_l)(x_i-mu_l)^Tp(lmid{x_i,Theta^g})}{sum_{i=1}^{N}p(lmid{x_i,Theta^g})} end{aligned}


推薦閱讀:

手把手教你快速構建自定義分類器
2017年歷史文章匯總|機器學習
大數據看190部國產片,年輕人的觀影口味發生了這些變化
通俗易懂說數據挖掘十大經典演算法
數據挖掘和網路爬蟲有什麼關聯區別?

TAG:機器學習 | 模式識別 | 數據挖掘 |