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，直到收敛。
后记：本篇博客主要讲了高斯混合模型的推导以及演算的过程，没有讲解相关的实用例子，在后面的博客中将会讲到使用高斯混合模型实现自动化肿瘤分割的任务。这里将其记录下来，留作后用。