EM算法推导以及在高斯混合模型中的应用（详细）

文章目录

EM算法的导出
EM算法的收敛性
EM算法在高斯混合模型中的应用
- 高斯混合模型
- 高斯混合模型参数估计的EM算法
- - 算法步骤
  - 算法推导
  - - 以一维情形为例p=1

参考教材：《统计学习方法》第二版——李航

EM算法的导出

EM算法通过不断求解下界函数的极大值来逼近求解对数似然函数的最大。

大致推导：
设一个含有隐变量Z的概率模型，对应有观测数据Y和待求解参数 θ \theta θ。现在的目标是最大化关于参数 θ \theta θ的对数似然函数，即：
L ( θ ) = l o g P ( Y ∣ θ ) = l o g ∑ Z P ( Y , Z ∣ θ ) = l o g ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) (1) L(\theta)=logP(Y|\theta)=log\sum_ZP(Y,Z|\theta)=log(\sum_Z P(Y|Z,\theta)P(Z|\theta)) \tag{1} L(θ)=logP(Y∣θ)=logZ∑P(Y,Z∣θ)=log(Z∑P(Y∣Z,θ)P(Z∣θ))(1)
极大化上式的困难在于最右边的式子中含有和（Z离散情况是和，连续情况下是积分）的对数形式。
EM算法考虑构造(1)中右侧式子的下界来逼近它，从而在不断的迭代逼近中达到一个较好的效果。
假设 θ i \theta_i θi是在第i次迭代得到的参数估计值。希望新的估计值 θ \theta θ能使得 L ( θ ) > L ( θ i ) L(\theta)>L(\theta_i) L(θ)>L(θi).考虑：
L ( θ ) − L ( θ i ) = l o g ( ∑ Z P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) − l o g P ( Y ∣ θ i ) = l o g ( ∑ Z P ( Z ∣ Y , θ i ) P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) P ( Z ∣ Y , θ i ) P ( Y ∣ θ i ) ) ≥ ∑ Z P ( Z ∣ Y , θ i ) l o g P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) P ( Z ∣ Y , θ i ) P ( Y ∣ θ i ) (2) \begin{aligned} L(\theta)-L(\theta_i)&=log(\sum_Z P(Y|Z,\theta)P(Z|\theta))-logP(Y|\theta_i) \\ &=log(\sum_Z P(Z|Y,\theta_i)\frac{P(Y|Z,\theta)P(Z|\theta))}{P(Z|Y,\theta_i)P(Y|\theta_i)}) \tag{2} \\ &\ge \sum_Z P(Z|Y,\theta_i)log \frac{P(Y|Z,\theta)P(Z|\theta))}{P(Z|Y,\theta_i)P(Y|\theta_i)} \\ \end{aligned} L(θ)−L(θi)=log(Z∑P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θi)=log(Z∑P(Z∣Y,θi)P(Z∣Y,θi)P(Y∣θi)P(Y∣Z,θ)P(Z∣θ)))≥Z∑P(Z∣Y,θi)logP(Z∣Y,θi)P(Y∣θi)P(Y∣Z,θ)P(Z∣θ))(2)

(2)中最后一个不等号利用到了对数函数的凹性质（或Jensen inequality），这也是(2)中第一个等式到第二个等式的推导中利用 P ( Z ∣ Y , θ i ) P(Z|Y,\theta_i) P(Z∣Y,θi)进行等价变换的原因，因为它是Z的条件概率，关于Z求和值一定为1且非负。这样就找到了对数似然函数 L ( θ ) L(\theta) L(θ)的下界：
B ( θ , θ i ) = L ( θ i ) + ∑ Z P ( Z ∣ Y , θ i ) l o g P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ i ) P ( Y ∣ θ i ) B(\theta,\theta_i)=L(\theta_i)+\sum_Z P(Z|Y,\theta_i)log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta_i)P(Y|\theta_i)} B(θ,θi)=L(θi)+∑ZP(Z∣Y,θi)logP(Z∣Y,θi)P(Y∣θi)P(Y∣Z,θ)P(Z∣θ)
并且利用条件概率计算公式可得 B ( θ i , θ i ) = L ( θ i ) B(\theta_i,\theta_i)=L(\theta_i) B(θi,θi)=L(θi)，这说明 B ( θ , θ i ) B(\theta,\theta_i) B(θ,θi)是 L ( θ ) L(\theta) L(θ)的一个不错的下界，至少在 θ i \theta_i θi上二者是相等的。也就是说最大化 B ( θ , θ i ) B(\theta,\theta_i) B(θ,θi)得到的最坏的结果也不会比当前迭代中得到的 θ i \theta_i θi更差。

Q函数的推导
Q函数其实已经求出来了，其实就是 B ( θ , θ i ) B(\theta,\theta_i) B(θ,θi)，但是为了求解方便需要进一步简化。
省去 B ( θ , θ i ) B(\theta,\theta_i) B(θ,θi)中与 θ \theta θ无关的项可得：
Q ( θ , θ i ) = ∑ Z P ( Z ∣ Y , θ i ) l o g ( P ( Y ∣ Z , θ ) P ( Z ∣ θ ) ) = ∑ Z P ( Z ∣ Y , θ i ) l o g P ( Y , Z ∣ θ ) = E [ l o g P ( y , Z ∣ θ ) ∣ y , θ i ] (3) \begin{aligned} Q(\theta,\theta_i)&=\sum_Z P(Z|Y,\theta_i)log(P(Y|Z,\theta)P(Z|\theta)) \\ &=\sum_Z P(Z|Y,\theta_i)logP(Y,Z|\theta) \\ &=E[logP(y,Z|\theta)|y,\theta_i] \tag{3} \end{aligned} Q(θ,θi)=Z∑P(Z∣Y,θi)log(P(Y∣Z,θ)P(Z∣θ))=Z∑P(Z∣Y,θi)logP(Y,Z∣θ)=E[logP(y,Z∣θ)∣y,θi](3)
(3)最后一个等号表明Q函数的本质是给定观测y和当前迭代步求得的参数下，对数似然函数的条件期望，实际的计算仍然常常采用第一个等号.
注：(3)这里第一个等式到第二个等式的推导与(1)中最后一个等号有关，在GMM模型中这一点使得(Y,Z)联合概率的计算非常方便。

（ending)

EM算法的收敛性

（待续

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代表y来自第k个分模型的概率； ∑ α k = 1 ; ϕ ( y ∣ θ k ) \sum\alpha_k=1; \phi(y|\theta_k) ∑αk=1;ϕ(y∣θk)是高斯分布密度， θ k = ( μ k , σ k 2 ) \theta_k=(\mu_k,\sigma^2_k) θk=(μk,σk2)

注：一般的混合模型可以是任意的分布密度，不局限于高斯密度。那么这类模型的参数EM算法求解过程都大致相同。

高斯混合模型参数估计的EM算法

算法步骤

输入：观测数据，高斯混合模型初始化参数
输出：高斯混合模型最终参数
(1) E步：求Q函数: Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q(θ,θ(i))
(2) M步： θ ( i + 1 ) = a r g m a x θ Q ( θ , θ ( i ) ) \theta^{(i+1)} =argmax_{\theta}Q(\theta,\theta^{(i)}) θ(i+1)=argmaxθQ(θ,θ(i))
(3) 重复(1),(2)步直至收敛。

算法推导

以一维情形为例p=1

观测数据 y i y_i yi的产生模式：首先依分布列 α k , k = 1 , 2 , . . . , K \alpha_k,k=1,2,...,K αk,k=1,2,...,K选择第k个高斯分布模型 ϕ ( y ∣ θ k ) \phi(y|\theta_k) ϕ(y∣θk)，由此生成观测数据 y i y_i yi。
已知：观测 y j , j = 1 , 2 , . . . , N . y_j, j=1,2,...,N. yj,j=1,2,...,N.
未知： θ = ( α 1 , . . . , α K , θ 1 , . . . , θ K ) . \theta=(\alpha_1,...,\alpha_K, \theta_1,...,\theta_K). θ=(α1,...,αK,θ1,...,θK).
设反映观测数据 y j y_j yj来自第k个分模型的潜在变量(0-1离散随机变量，即两点分布)为 γ j k \gamma_{jk} γjk,在 j = 1 , 2 , . . . , N j=1,2,...,N j=1,2,...,N之间独立， ∑ k γ j k = 1 , ∑ k P ( γ j k = 1 ) = ∑ k α k = 1 \sum_k\gamma_{jk}=1,\sum_kP(\gamma_{jk}=1)=\sum_k\alpha_k=1 ∑kγjk=1,∑kP(γjk=1)=∑kαk=1.
完全数据为: ( y j , γ j 1 , γ j 2 , . . . , γ j K ) , j = 1 , 2 , . . . , N (y_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jK}),j=1,2,...,N (yj,γj1,γj2,...,γjK),j=1,2,...,N
完全数据的似然函数：
P ( y , γ ∣ θ ) = ∏ j = 1 N P ( y j , γ j 1 , γ j 2 , . . . , γ j K ∣ θ ) = ∏ k = 1 K ∏ j = 1 N [ α k ϕ ( y j ∣ θ k ) ] γ j k = ∏ k = 1 K α k n k ∏ j = 1 N [ ϕ ( y j ∣ θ k ) ] γ j k (4) \begin{aligned} P(y,\gamma|\theta)&=\prod_{j=1}^{N}P(y_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jK}|\theta)\\ &=\prod_{k=1}^{K}\prod_{j=1}^{N}[\alpha_k \phi(y_j|\theta_k)]^{\gamma_{jk}} \tag{4} \\ &=\prod_{k=1}^{K}\alpha_k^{n_k}\prod_{j=1}^{N}[ \phi(y_j|\theta_k)]^{\gamma_{jk}} \\ \end{aligned} P(y,γ∣θ)=j=1∏NP(yj,γj1,γj2,...,γjK∣θ)=k=1∏Kj=1∏N[αkϕ(yj∣θk)]γjk=k=1∏Kαknkj=1∏N[ϕ(yj∣θk)]γjk(4)
其中， n k = ∑ j γ j k ; ∑ k n k = N n_k=\sum_j \gamma_{jk} ; \sum_kn_k=N nk=∑jγjk;∑knk=N
ϕ ( y j ∣ θ k ) = 1 2 π σ k e x p ( − ( y i − μ k ) 2 2 σ k 2 ) \phi(y_j|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k}exp(-\frac{(y_i-\mu_k)^2}{2\sigma^2_k}) ϕ(yj∣θk)=2π σk1exp(−2σk2(yi−μk)2)
对数似然函数：
l o g P ( y , γ ∣ θ ) = ∑ k = 1 K { n k l o g α k + ∑ j = 1 N γ j k l o g [ ϕ ( y j ∣ θ k ) ] } logP(y,\gamma|\theta)=\sum_{k=1}^{K}\{{n_k}log\alpha_k+\sum_{j=1}^{N}\gamma_{jk}log[ \phi(y_j|\theta_k)]\} logP(y,γ∣θ)=∑k=1K{nklogαk+∑j=1Nγjklog[ϕ(yj∣θk)]}

E步求Q函数：(书上的推导公式)
这时记第i次迭代得到的参数估计值为 θ ( i ) \theta^{(i)} θ(i).
E [ l o g P ( y , γ ∣ θ ) ∣ y , θ ( i ) ] = ∑ k = 1 K { E ( n k ) l o g α k + ∑ j = 1 N E ( γ j k ) l o g [ ϕ ( y j ∣ θ k ) ] } (5) E[logP(y,\gamma|\theta)|y,\theta^{(i)}]=\sum_{k=1}^{K}\{{E(n_k)}log\alpha_k+\sum_{j=1}^{N}E(\gamma_{jk})log[ \phi(y_j|\theta_k)] \} \tag{5} E[logP(y,γ∣θ)∣y,θ(i)]=k=1∑K{E(nk)logαk+j=1∑NE(γjk)log[ϕ(yj∣θk)]}(5)
注意此处以及以下的E(.)表示条件期望，记 γ ^ j k = E ( γ j k ∣ y , θ ( i ) ) \hat{\gamma}_{jk}=E(\gamma_{jk}|y,\theta^{(i)}) γ^jk=E(γjk∣y,θ(i))。由(4)中的推导知， n k n_k nk与 γ k \gamma_k γk有关， E ( n k ) = ∑ j E ( γ j k ∣ y , θ ( i ) ) E(n_k)=\sum_jE(\gamma_{jk}|y,\theta^{(i)}) E(nk)=∑jE(γjk∣y,θ(i))，
γ ^ j k = E ( γ j k ∣ y , θ ( i ) ) = P ( γ j k = 1 ∣ y , θ ( i ) ) = P ( γ j k = 1 , y j ∣ θ ( i ) ) ∑ k P ( γ j k = 1 , y j ∣ θ ( i ) ) = P ( y j ∣ γ j k = 1 , θ ( i ) ) P ( γ j k = 1 ∣ θ ( i ) ) ∑ k P ( y j ∣ γ j k = 1 , θ ( i ) ) P ( γ j k = 1 ∣ θ ( i ) ) = α k ( i ) ϕ ( y j ∣ θ k ( i ) ) ∑ k α k ( i ) ϕ ( y j ∣ θ k ( i ) ) (6) \begin{aligned} \hat{\gamma}_{jk}&=E(\gamma_{jk}|y,\theta^{(i)}) \\ &=P(\gamma_{jk}=1|y,\theta^{(i)}) \tag{6} \\ &=\frac{P(\gamma_{jk}=1,y_j|\theta^{(i)})}{\sum_kP(\gamma_{jk}=1,y_j|\theta^{(i)})} \\ &=\frac{P(y_j|\gamma_{jk}=1,\theta^{(i)})P(\gamma_{jk}=1|\theta^{(i)})}{\sum_kP(y_j|\gamma_{jk}=1,\theta^{(i)})P(\gamma_{jk}=1|\theta^{(i)})} \\ &=\frac{\alpha^{(i)}_{k} \phi(y_j|\theta_k^{(i)})}{\sum_k\alpha^{(i)}_{k} \phi(y_j|\theta_k^{(i)})} \\ \end{aligned} γ^jk=E(γjk∣y,θ(i))=P(γjk=1∣y,θ(i))=∑kP(γjk=1,yj∣θ(i))P(γjk=1,yj∣θ(i))=∑kP(yj∣γjk=1,θ(i))P(γjk=1∣θ(i))P(yj∣γjk=1,θ(i))P(γjk=1∣θ(i))=∑kαk(i)ϕ(yj∣θk(i))αk(i)ϕ(yj∣θk(i))(6)

(6)中第2个等号到第3个等号的推导：
P ( γ j k = 1 ∣ y , θ ( i ) ) = P ( γ j k = 1 , y j , θ ( i ) ) P ( θ ( i ) ) P ( θ ( i ) ) P ( y j , θ ( i ) ) = P ( γ j k = 1 , y j ∣ θ ( i ) ) P ( y j ∣ θ ( i ) ) (7) P(\gamma_{jk}=1|y,\theta^{(i)})=\frac{P(\gamma_{jk}=1,y_j,\theta^{(i)})}{P(\theta^{(i)})} \frac{P(\theta^{(i)})}{P(y_j,\theta^{(i)})}=\frac{P(\gamma_{jk}=1,y_j|\theta^{(i)} )}{P(y_j|\theta^{(i)})} \tag{7} P(γjk=1∣y,θ(i))=P(θ(i))P(γjk=1,yj,θ(i))P(yj,θ(i))P(θ(i))=P(yj∣θ(i))P(γjk=1,yj∣θ(i))(7)
P ( y j ∣ θ ( i ) ) = ∑ k P ( γ j k = 1 , y j ∣ θ ( i ) ) (8) P(y_j|\theta^{(i)})=\sum_kP(\gamma_{jk}=1,y_j|\theta^{(i)}) \tag{8} P(yj∣θ(i))=k∑P(γjk=1,yj∣θ(i))(8)
根据条件分布定义进行推导就OK。

γ ^ j k \hat{\gamma}_{jk} γ^jk的含义在书中被描述为“在当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据 y j y_j yj的响应度”。结合(5)式不难理解，其中既有“先验信息”——概率 α k \alpha_k αk，还包含“当前的信息”——观测数据的密度。

把 γ ^ j k = α k ( i ) ϕ ( y j ∣ θ k ( i ) ) ∑ k α k ( i ) ϕ ( y j ∣ θ k ( i ) ) \hat{\gamma}_{jk}=\frac{\alpha^{(i)}_{k} \phi(y_j|\theta_k^{(i)})}{\sum_k\alpha^{(i)}_{k} \phi(y_j|\theta_k^{(i)})} γ^jk=∑kαk(i)ϕ(yj∣θk(i))αk(i)ϕ(yj∣θk(i))带入(5)即得到Q函数：
Q ( θ , θ ( i ) ) = E [ l o g P ( y , γ ∣ θ ) ∣ y , θ ( i ) ] = ∑ k = 1 K { ∑ j = 1 N γ ^ j k l o g α k + ∑ j = 1 N γ ^ j k l o g [ ϕ ( y j ∣ θ k ) ] } Q(\theta,\theta^{(i)})=E[logP(y,\gamma|\theta)|y,\theta^{(i)}]=\sum_{k=1}^{K} \{ \sum_{j=1}^{N}\hat{\gamma}_{jk} log\alpha_k+\sum_{j=1}^{N}\hat{\gamma}_{jk}log[ \phi(y_j|\theta_k)] \} Q(θ,θ(i))=E[logP(y,γ∣θ)∣y,θ(i)]=k=1∑K{j=1∑Nγ^jklogαk+j=1∑Nγ^jklog[ϕ(yj∣θk)]}
其中， l o g ϕ ( y j ∣ θ k ) = l o g 1 2 π − l o g σ k − ( y i − μ k ) 2 2 σ k 2 log\phi(y_j|\theta_k)=log\frac{1}{\sqrt{2\pi}}-log\sigma_k-\frac{(y_i-\mu_k)^2}{2\sigma^2_k} logϕ(yj∣θk)=log2π 1−logσk−2σk2(yi−μk)2.

M步：最大化Q函数
求偏导，令其为0即可解得新一次迭代的参数值：

∂ Q ∂ μ k = ∑ j γ ^ j k y j − μ k σ k 2 = 0 , μ ^ k = ∑ j γ ^ j k y j ∑ j γ ^ j k \frac{\partial Q}{\partial \mu_k}=\sum_j \hat{\gamma}_{jk}\frac{y_j-\mu_k}{\sigma^2_k } =0, \hat{\mu}_k=\frac{\sum_j \hat{\gamma}_{jk}y_j}{\sum_j \hat{\gamma}_{jk}} ∂μk∂Q=∑jγ^jkσk2yj−μk=0,μ^k=∑jγ^jk∑jγ^jkyj
∂ Q ∂ σ k 2 = ∑ j γ ^ j k ( − 1 σ k 2 1 2 σ k 2 + ( y i − μ k ) 2 2 σ k 4 ) = 0 , σ ^ k 2 = ∑ j γ ^ j k ( y i − μ k ) 2 ∑ j γ ^ j k \frac{\partial Q}{\partial \sigma^2_k}=\sum_j\hat{\gamma}_{jk}(-\frac{1}{\sqrt{\sigma^2_k}}\frac{1}{2\sqrt{\sigma^2_k}}+\frac{(y_i-\mu_k)^2}{2\sigma^4_k})=0, \hat{\sigma}_k^2=\frac{\sum_j\hat{\gamma}_{jk}(y_i-\mu_k)^2 }{\sum_j\hat{\gamma}_{jk}} ∂σk2∂Q=∑jγ^jk(−σk2 12σk2 1+2σk4(yi−μk)2)=0,σ^k2=∑jγ^jk∑jγ^jk(yi−μk)2
记 n k = ∑ j = 1 N γ ^ j k n_k=\sum_{j=1}^{N}\hat{\gamma}_{jk} nk=∑j=1Nγ^jk
∂ Q ∂ α i = ∂ ∂ α i ∑ k = 1 K n k l o g α k = ∂ ∂ α i ( ∑ k ≠ 1 n k l o g α k + n 1 l o g ( 1 − ∑ k ≠ 1 α k ) ) = n i α i − n 1 α 1 = 0 \frac{\partial Q}{\partial \alpha_i}=\frac{\partial }{\partial \alpha_i}\sum_{k=1}^{K} n_klog\alpha_k =\frac{\partial }{\partial \alpha_i}(\sum_{k\ne 1} n_klog\alpha_k+n_1 log(1-\sum_{k\ne 1}\alpha_k ))=\frac{n_i}{\alpha_i}-\frac{n_1}{\alpha_1}=0 ∂αi∂Q=∂αi∂∑k=1Knklogαk=∂αi∂(∑k=1nklogαk+n1log(1−∑k=1αk))=αini−α1n1=0
n i α 1 = n 1 α i n_i \alpha_1=n_1\alpha_i niα1=n1αi, N α 1 = n 1 N \alpha_1=n_1 Nα1=n1, α ^ 1 = n 1 N \hat{\alpha}_1 =\frac{n_1}{N} α^1=Nn1
同理， α ^ k = n k N \hat{\alpha}_k =\frac{n_k}{N} α^k=Nnk

代码部分等待更新

补充多维情形
E步相似，M步也相似。M步具体的参数更新式子如下：
γ ^ j k = γ ^ j k = α k ( i ) ϕ ( Y j ∣ θ k ( i ) ) ∑ k α k ( i ) ϕ ( Y j ∣ θ k ( i ) ) \hat{\gamma}_{jk}=\hat{\gamma}_{jk}=\frac{\alpha^{(i)}_{k} \phi(Y_j|\theta_k^{(i)})}{\sum_k\alpha^{(i)}_{k} \phi(Y_j|\theta_k^{(i)})} γ^jk=γ^jk=∑kαk(i)ϕ(Yj∣θk(i))αk(i)ϕ(Yj∣θk(i))，此处 θ k ( i ) = ( μ k , Σ k ) \theta^{(i)}_k=(\mu_k,\Sigma_k) θk(i)=(μk,Σk)
1、 μ ^ k = ∑ j γ ^ j k Y j ∑ j γ ^ j k \hat{\mu}_k=\frac{\sum_j \hat{\gamma}_{jk}Y_j}{\sum_j \hat{\gamma}_{jk}} μ^k=∑jγ^jk∑jγ^jkYj
2、 σ ^ k 2 = ∑ j γ ^ j k ( Y j − μ k ) ( Y j − μ k ) T ∑ j γ ^ j k \hat{\sigma}_k^2 =\frac{\sum_j\hat{\gamma}_{jk}(Y_j-\mu_k)(Y_j-\mu_k)^T}{\sum_j\hat{\gamma}_{jk}} σ^k2=∑jγ^jk∑jγ^jk(Yj−μk)(Yj−μk)T
3、 α ^ k = 1 N ∑ j γ ^ j k \hat{\alpha}_k =\frac{1}{N}\sum_j \hat{\gamma }_{jk} α^k=N1∑jγ^jk