摘自《统计学习方法》 李航著 清华大学出版社
链接：
EM介绍请点击这里

EM算法在高斯混合模型中的应用

EM算法的一个重要应用是高斯混合模型的参数估计。高斯混合模型应用广泛，在许多情况下，EM算法是学习高斯混合模型的有效方法。想详细理解单高斯模型和混合高斯模型请点击这里

高斯混合模型

定义

高斯混合模型是指具有如下形式的概率分布模型：

P(y|Θ)=∑Kk=1αϕ(y|Θk) 公式（1）

ϕ(y|Θk)=12π√σkexp−(y−μk)22σ2k 公式（2）

---

其中，α是系数，α≥0，∑Kk=1αk=1；ϕ(y|Θk)是高斯分布密度，Θ=(μk,σk)。公式（2）我们称为第K个分模型。

---

一般混合模型可以由任意概率分布密度代替公式（2）中的高斯分布密度，我们只介绍最常用的高斯混合模型。

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

假设观测数据y1,y2,...,yN由高斯混合模型生成，

P(y|Θ)=∑Kk=1ϕ(y|Θk) 公式（3）

---

其中，Θ=(α1,α2,...,αk;θ1,θ2,...,θk)。我们用EM算法估计高斯混合模型的参数θ。

---

明确隐变量，写出完全数据的对数似然函数

可以设想观测数据yi， j=1,2,...,N，是这样产生的：首先依概率αk选择第k个高斯分布分模型ϕ(y|θk)；然后依第k个分模型的概率分布ϕ(y|θk)生成观测数据yj。这是观测数据yj，j=1,2,...,N，是已只的；反映观测数据yj来自第k个分模型的数据是未知的，k=1,2,...,K，以隐变量γjk表示，其定义如下：

γk={1第j个观察来自第k个分模型0 否则

j=1,2,...,N;k=1,2,....,K 公式（4）

---

其中γjk是0-1随机变量。

---

有了观测数据yj及未观测数据γjk，那么完全数据是

(yj,γj1,γj2,...,γjK)，j=1,2,...,N

于是可以写出完全数据的似然函数：

P(y,γ|θ)=∏j=1NP(yj,γj1,γj2,...,γjk|θ)

=∏Kk=1∏Nj=1[αkϕ(yj|θk)]γjk 

=∏Kk=1αnk∏Nj=1[ϕ(yj|θk)]γjk  

=∏αnkk∏Nj=1[12π√σkexp(−(yj−μk)22σ2k)]γjk  

---

式中，nk=∑Nj=1γjk，∑Kk=1nk=N
公式注释：
上述第二个等号为什么成立？我们思考一下。在上节EM算法介绍中，我们在硬币例题介绍了

P(y|θ)=∑zP(y,z|θ)=∑zP(z|θ)P(y|z,θ)

现在，我们来对应一下。P(z|θ)对应高斯混合模型的P(γjk|θ)，也就是αk。在第k个模型的情况下，观测值yi属于第k个模型的概率是ϕ(yj|θk)，也就是ϕ(y|Θk)=12π√σkexp−(y−μk)22σ2k，然而不只是由k这种模型，还有1~K个模型，所以是每个模型算出概率后相乘。注意γjk的值不是1就是0，所以当yi不属于第k类的时候，任何数的0次方为1。当当yi属于第k类的时候，γjk的值为1，任何数的1次方还是其本身。

---

那么，完全数据的对数似然函数为:

logP(y,γ|θ)=∑Kk=1nk{logαk+∑Nj=1γjk[log(12π√)−logσk−12σ2k(yj−μk)2]}

EM算法的E步：确定Q函数

Q(θ,θ(i))=E(logP(y,γ|θ)|y,θ(i))

=E{∑Kk=1{nklogαk+∑Nj=1γjk[log12π√−logσk−12σ2k(yj−μk)2]}}

=∑Kk=1{∑Nj=1(Eγjk)logαk+∑Nj=1(Eγrk)[log(1(√2π))−logσk−12σ2(yj−μk)2]}

公式（4）

---

从第2个等号到第3个等号这里涉及到一些期望的知识，以图片的形式补充到下面。

资料来源：http://www.360doc.com/content/13/1124/03/9482_331690142.shtml

---

这里需要计算E(γjk|y,θ)，记为γ‘jk。

γjk‘=E(γjk|y,θ)=P(γrk=1|y,θ)

=P(γ=1,yj|θ)∑Kk=1P(γjk=1,yj|θ) 

=P(yj|γjk=1,θ)P(γjk=1|θ)∑Kk=1P(yj|γjk=1,θ)P(γjk=1|θ)  

=αkϕ(yj|θk)∑Kk=1αkϕ(yj|θk)j=1, 2,…, N; k = 1, 2, …, K

γ‘jk是在当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据yj的响应度。
将γ‘jk=Eγjk及nk=∑Nj=1Eγjk代入公式（4）（确定Q函数的公式），即可获得

Q(θ,θ(i))=∑Kk=1{nklogαk+∑Nk=1γ‘jk[log(12π√)−logσk−12σ2k(yj−μk)2]} 公式（5）

---

确定EM算法的M步

迭代的M步是求函数Q(θ,θ(i))对θ的极大值，即求新一轮迭代的模型参数：

θ(i+1)=argmaxQ(θ,θ(i))

用μk^,σ2k^及αk^，k=1,2,...,K，表示θ(i+1)的各个参数。求μk^，σ2k^只需将公式（5）分别对μk^，σ2k^求偏导数并令其为0，即可得到；求αk^是在∑Kk=1αk=1条件下求偏导并令其为0得到的。其结果如下所示：

μk^=∑Nj=1γjkyj^∑Nj=1γjk^k=1,2,...,K公式（6）

σ2k^=∑Nj=1γjk^(yi−uk)2∑Nj=1γjk^k=1,2,...,K   公式（7）

αk^=nkN=∑Nj=1γjk^Nk=1,2,...,K   公式（7）

重复以上计算，直到对数似然函数值不再有明显变化为止。

---

算法总结

输入：观测数据y1,y2,...,yN，高斯混合模型；
输出：高斯混合模型参数；
步骤：
（1）取参数的初始值开始迭代
（2）E步：依据当前模型参数，计算分模型k对观测数据yi的响应度

γk^=αkϕ(yj|θk)∑Kk=1αkϕ(yj|θk)j=1,2,...,N; k=1,2,...,K

（3）M步：计算新一轮迭代的模型参数:

μk^=∑Nj=1γjkyj^∑Nj=1γjk^k=1,2,...,K

σ2k^=∑Nj=1γjk^(yi−uk)2∑Nj=1γjk^k=1,2,...,K   

αk^=nkN=∑Nj=1γjk^N k=1,2,...,K 

（4）重复第（2）步和第（3）步，直到收敛。