1. 前言
高斯混合模型是使用高斯分布對原始資料進行估計，其中高斯函數的均值 μ \mu μ和方差 σ \sigma σ以及各個高斯函數分量占的比例 α \alpha α這些參數是未知的。對于它們的求解是通過EM算法實作的。
高斯混合模型可以表述為如下的機率分布模型：
P ( y ∣ θ ) = ∑ k = 1 K α k ϕ ( y ∣ θ k ) P(y|\theta)=\sum_{k=1}^{K}\alpha_k\phi(y|\theta_k) P(y∣θ)=k=1∑Kαkϕ(y∣θk)
其中， α k \alpha_k αk是系數， α k ≥ 0 \alpha_k\ge0 αk≥0， ∑ k = 1 K α k = 1 \sum_{k=1}^K\alpha_k=1 ∑k=1Kαk=1； ϕ ( y ∣ θ k ) \phi(y|\theta_k) ϕ(y∣θk)是高斯分布密度， θ k = ( μ k , σ k 2 ) \theta_k=(\mu_k,\sigma_k^2) θk=(μk,σk2)
ϕ ( y ∣ θ k ) = 1 2 π σ k e x p ( − ( y − μ k ) 2 2 σ k 2 ) \phi(y|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k}exp(-\frac{(y-\mu_k)^2}{2\sigma_k^2}) ϕ(y∣θk)=2π
σk1exp(−2σk2(y−μk)2)
稱為第K個分模型。
2. 推導
2.1 确定隐藏參數
假設觀測資料 y 1 , y 2 , … , y N y_1,y_2,\dots,y_N y1,y2,…,yN是由高斯混合模型生成
P ( y ∣ θ ) = ∑ k = 1 K α k ϕ ( y ∣ θ k ) P(y|\theta)=\sum_{k=1}^{K}\alpha_k\phi(y|\theta_k) P(y∣θ)=k=1∑Kαkϕ(y∣θk)
其中，參數 θ = ( α 1 , α 2 , … , α k ; θ 1 , θ 2 , … , θ k ) \theta=(\alpha_1,\alpha_2,\dots,\alpha_k;\theta_1,\theta_2,\dots,\theta_k) θ=(α1,α2,…,αk;θ1,θ2,…,θk)。是高斯混合模型中直覺需要求解的未知變量。但是高斯混合模型中就隻有這兩類變量是未知的麼？
對于觀測資料 y j y_j yj， j = 1 , 2 , … , N j=1,2,\dots,N j=1,2,…,N，是這樣産生的：首先依據機率 α k \alpha_k αk選擇第 k k k個高斯分布模型 ϕ ( y ∣ θ k ) \phi(y|\theta_k) ϕ(y∣θk)；然後依據第 k k k個分模型的機率分布 ϕ ( y ∣ θ k ) \phi(y|\theta_k) ϕ(y∣θk)生成觀測資料 y j y_j yj。這個時候觀測資料 y j y_j yj， j = 1 , 2 , … , N j=1,2,\dots,N j=1,2,…,N，是已知的；但是呢，反映觀測資料 y j y_j yj具體是來自第 k k k個分模型的資料是未知的。 k = 1 , 2 , … , K k=1,2,\dots,K k=1,2,…,K，這裡以隐變量 γ j k \gamma_{jk} γjk表示，可以将其定義為：
γ j k = { 1 , 第 j 個觀測資料來自第 k 個分模型 0 , 否則 \gamma_{jk}= \begin{cases} 1, & \text{第$j$個觀測資料來自第$k$個分模型} \\ 0, & \text{否則} \end{cases} γjk={1,0,第j個觀測資料來自第k個分模型否則
其中， j = 1 , 2 , … , N ; k = 1 , 2 , … , K j=1,2,\dots,N;k=1,2,\dots,K j=1,2,…,N;k=1,2,…,K，顯然 γ j k \gamma_{jk} γjk是屬于典型的0-1随機變量。因而完整的觀測資料就可以寫為：
( y j , γ j 1 , γ j 2 , … , γ j K ) , j = 1 , 2 , … , N (y_j, \gamma_{j1},\gamma_{j2},\dots,\gamma_{jK}),j=1,2,\dots,N (yj,γj1,γj2,…,γjK),j=1,2,…,N
2.2 得到似然函數，确定Q函數
于是，可以寫出完全資料的似然函數：
P ( y , γ ∣ θ ) = ∏ j = 1 N P ( y j , γ j 1 , γ j 2 , … , γ j K ) = ∏ j = 1 N ∏ k = 1 K [ α k ϕ ( y j ∣ θ k ) ] γ j k = ∏ k = 1 K α k n k ∏ j = 1 N [ ϕ ( y j ∣ θ k ) ] γ j k = ∏ k = 1 K α k n k ∏ j = 1 N [ 1 2 π σ k e x p ( − ( y − μ k ) 2 2 σ k 2 ) ] γ j k \begin{aligned} P(y,\gamma|\theta) & = \prod_{j=1}^NP(y_j,\gamma_{j1},\gamma_{j2},\dots,\gamma_{jK}) \\ & = \prod_{j=1}^N\prod_{k=1}^K[\alpha_k\phi(y_j|\theta_k)]^{\gamma_{jk}}\\ & =\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\phi(y_j|\theta_k)]^{\gamma_{jk}} \\ & =\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\frac{1}{\sqrt{2\pi}\sigma_k}exp(-\frac{(y-\mu_k)^2}{2\sigma_k^2})]^{\gamma_{jk}} \end{aligned} P(y,γ∣θ)=j=1∏NP(yj,γj1,γj2,…,γjK)=j=1∏Nk=1∏K[αkϕ(yj∣θk)]γjk=k=1∏Kαknkj=1∏N[ϕ(yj∣θk)]γjk=k=1∏Kαknkj=1∏N[2π
σk1exp(−2σk2(y−μk)2)]γjk
其中， n k = ∑ j = 1 N γ j k n_k=\sum_{j=1}^{N}\gamma_{jk} nk=∑j=1Nγjk， N = ∑ k = 1 K n k N=\sum_{k=1}^{K}n_k N=∑k=1Knk
因而，完全資料的對數似然函數為：
l o g P ( y , γ ∣ θ ) = ∑ k = 1 K { n k l o g α k + ∑ j = 1 N γ j k [ l o g ( 1 2 π ) − l o g α k − 1 2 π ( y j − μ k ) 2 ] } logP(y,\gamma|\theta)=\sum_{k=1}^K \{n_klog\alpha_k+\sum_{j=1}^N\gamma_{jk}[log(\frac{1}{\sqrt{2\pi}})-log\alpha_k-\frac{1}{\sqrt{2\pi}}(y_j-\mu_k)^2]\} logP(y,γ∣θ)=k=1∑K{nklogαk+j=1∑Nγjk[log(2π
1)−logαk−2π
1(yj−μk)2]}
在上面得到了完全資料的對數似然函數，接下來就需要使用EM算法求解模型中的未知參數了。首先是EM算法的E步，也就是求解期望，得到Q函數。
Q ( θ , θ ( i ) ) = E [ l o g P ( y , γ ∣ θ ) ∣ y , θ ( i ) ] = E { ∑ k = 1 K { n k l o g α k + ∑ j = 1 N γ j k [ l o g ( 1 2 π ) − l o g α k − 1 2 π ( y j − μ k ) 2 ] } } = ∑ k = 1 K { ∑ j = 1 N ( E γ j k ) l o g α k + ∑ j = 1 N ( E γ j k ) [ l o g ( 1 2 π ) − l o g α k − 1 2 π ( y j − μ k ) 2 ] } \begin{aligned} Q(\theta,\theta^{(i)}) & = E[logP(y,\gamma|\theta)|y,\theta^{(i)}]\\ & = E\{\sum_{k=1}^K \{n_klog\alpha_k+\sum_{j=1}^N\gamma_{jk}[log(\frac{1}{\sqrt{2\pi}})-log\alpha_k-\frac{1}{\sqrt{2\pi}}(y_j-\mu_k)^2]\}\}\\ & =\sum_{k=1}^K \{\sum_{j=1}^N(E\gamma_{jk})log\alpha_k+\sum_{j=1}^N(E\gamma_{jk})[log(\frac{1}{\sqrt{2\pi}})-log\alpha_k-\frac{1}{\sqrt{2\pi}}(y_j-\mu_k)^2]\} \\ \end{aligned} Q(θ,θ(i))=E[logP(y,γ∣θ)∣y,θ(i)]=E{k=1∑K{nklogαk+j=1∑Nγjk[log(2π
1)−logαk−2π
1(yj−μk)2]}}=k=1∑K{j=1∑N(Eγjk)logαk+j=1∑N(Eγjk)[log(2π
1)−logαk−2π
1(yj−μk)2]}
這裡記 γ j k ^ = E ( γ j k ∣ y , θ ) \hat{\gamma_{jk}}=E(\gamma_{jk}|y,\theta) γjk^=E(γjk∣y,θ)，則：
γ j k ^ = α k ϕ ( y j ∣ θ k ) ∑ k = 1 K α k ϕ ( y j ∣ θ k ) \hat{\gamma_{jk}}=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^{K}\alpha_k\phi(y_j|\theta_k)} γjk^=∑k=1Kαkϕ(yj∣θk)αkϕ(yj∣θk)
γ j k ^ \hat{\gamma_{jk}} γjk^是在目前模型參數下第 j j j個觀測資料來自第 k k k個分模型的機率，成為分模型k對觀測資料 y j y_j yj的響應度。
将 γ j k ^ = E γ j k \hat{\gamma_{jk}}=E\gamma_{jk} γjk^=Eγjk及 n k = ∑ k = 1 N E γ j k n_k=\sum_{k=1}^{N}E\gamma_{jk} nk=∑k=1NEγjk，可得：
Q ( θ , θ ( i ) ) = ∑ k = 1 K { n k l o g α k + ∑ j = 1 N γ j k ^ [ l o g ( 1 2 π ) − l o g α k − 1 2 π ( y j − μ k ) 2 ] } Q(\theta,\theta^{(i)})=\sum_{k=1}^K \{n_klog\alpha_k+\sum_{j=1}^N\hat{\gamma_{jk}}[log(\frac{1}{\sqrt{2\pi}})-log\alpha_k-\frac{1}{\sqrt{2\pi}}(y_j-\mu_k)^2]\} Q(θ,θ(i))=k=1∑K{nklogαk+j=1∑Nγjk^[log(2π
1)−logαk−2π
1(yj−μk)2]}
2.3 EM算法M步
之後，便是使用EM算法的M步求 Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q(θ,θ(i))對 θ \theta θ的極大值，即求新一輪疊代的模型參數
θ i + 1 = a r g max θ Q ( θ , θ ( i ) ) \theta^{i+1}=arg\max_{\theta}Q(\theta,\theta^{(i)}) θi+1=argθmaxQ(θ,θ(i))
這裡用 μ k ^ , σ k 2 ^ , α k ^ \hat{\mu_k},\hat{\sigma_k^2},\hat{\alpha_k} μk^,σk2^,αk^來表示 θ i + 1 \theta^{i+1} θi+1中的各參數。對于參數 μ k ^ , σ k 2 ^ \hat{\mu_k},\hat{\sigma_k^2} μk^,σk2^可以通過上述的Q函數對應求導取極值得到；對于 α k ^ \hat{\alpha_k} αk^是在 ∑ k = 1 K = 1 \sum_{k=1}^K=1 ∑k=1K=1條件下求偏導數并令其為0得到。因而可以得到如下的參數更新公式
3. 高斯混合模型參數估計的EM算法
輸入：觀測資料 y 1 , y 2 , … , y N y_1,y_2,\dots,y_N y1,y2,…,yN，高斯混合模型
輸出：高斯混合模型中的參數
step1：取參數的初始值開始疊代
step2：E步，依據目前模型的參數，計算分模型 k k k對觀測資料 y j y_j yj的響應度
γ j k ^ = α k ϕ ( y j ∣ θ k ) ∑ k = 1 K α k ϕ ( y j ∣ θ k ) \hat{\gamma_{jk}}=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^{K}\alpha_k\phi(y_j|\theta_k)} γjk^=∑k=1Kαkϕ(yj∣θk)αkϕ(yj∣θk)
step3：M步，計算新一輪疊代的模型參數
step4：重複step2和step3，直到收斂。
後記：本篇部落客要講了高斯混合模型的推導以及演算的過程，沒有講解相關的實用例子，在後面的部落格中将會講到使用高斯混合模型實作自動化惡性良性腫瘤分割的任務。這裡将其記錄下來，留作後用。