EM 算法和高斯混合模型

《统计学习方法》对 EM 算法和 GMM 的介绍都比较详细，但有不少地方绕弯路，而且记号稍显混乱；《西瓜书》对 EM 算法和 GMM 的介绍都非常简洁，有许多细节被略去；B 站上的白板推导系列里也有 EM 算法和 GMM 详细推导。

为了这次数据挖掘小组作业，在这里整理一下 EM 和 GMM。

EM 算法是解决含有隐变量情况下的极大似然估计算法

极大似然估计

已知独立同分布数据集 X={x1,...,xn}，即 X∼i.i.dp(x;Θ)。它们的联合概率密度记为 p(x;Θ)，其中 x=(x1,...,xn)。对参数 Θ 的极大似然估计化归为求解如下问题 Θ^MLE=arg⁡maxθlnp(x;Θ)=arg⁡maxθ∑i=1nlnp(xi;Θ) 如果含有隐变量 z∼i.i.dq(z;Θ)，隐含数据集记为 Z={z1,...,zn}，它们的联合概率密度记为 q(z;Θ)，其中 z=(z1,...,zn)，此时 MLE 变为求解如下问题 Θ^MLE=arg⁡maxΘlnp(x;Θ)=arg⁡maxΘln∫zp(x,z;Θ)dz=arg⁡maxΘln∫zp(z;Θ)p(x|z;Θ)dz=arg⁡maxΘ∫zp(z;Θ)p(x|z;Θ)dz 因隐含数据集 Z={z1,...,zn} 无法观测到，联合概率密度 p(z;Θ) 无法计算，无法求解上述问题

对数似然函数的下界

在有隐变量 Z 的情况下，为了最大化对数似然函数 lnp(x;Θ)，先对它做一些变形

假设 q(z) 是隐变量 Z 服从的概率分布，q(z) 是对应的联合概率密度，则有 lnp(x;Θ)=ln∫zp(x,z;Θ)dz=ln∫zp(x,z;Θ)q(z)q(z)dz 然后用 Jensen 不等式进行放缩 lnp(x;Θ)=ln∫zp(x,z;Θ)q(z)q(z)dz≥∫zlnp(x,z;Θ)q(z)q(z)dz 这就得到了对数似然函数 lnp(x;Θ) 的下界 lnp(x;Θ)≥∫zlnp(x,z;Θ)q(z)q(z)dz

由 Jensen 不等式的性质知，等号成立当且仅当 p(x,z;Θ)=Cq(z) 两边对 z 积分可得 ∫zp(x,z;Θ)dz=C∫zq(z)dz 因为 q(z) 是概率密度函数，积分是 1；p(x,z;Θ) 对 z 积分是边缘概率密度 p(x;Θ)。由此可知 C=p(x;Θ)，所以 q(z)=p(x,z;Θ)C=p(x,z;Θ)p(x;Θ)=p(z|x;Θ) 即 q(z) 是观测到数据集 x 后的后验分布。这也启发了 EM 算法选取 q(z) 的原因。

EM 算法的导出

我们把真实参数记为 Θ，EM 算法第 t 步迭代时得到的参数近似值记为 Θ(t)

EM 算法把第 t 步迭代时以当前参数估计值 Θ(t) 的 z 的后验分布作为 q(z)，即 q(z)=p(z|x;Θ(t))，这样就得到了第 t 步迭代时对数似然的下界，记为 B(Θ,Θ(t))lnp(x;Θ)≥∫zlnp(x,z;Θ)p(z|x;Θ(t))p(z|x;Θ(t))dz=B(Θ,Θ(t)) 为了最大化对数似然函数，EM 算法对它的下界 B(Θ,Θ(t)) 进行最大化，求解得到的参数作为下一步迭代的参数 Θ(t+1)Θ(t+1)=arg⁡maxΘB(Θ,Θ(t))=arg⁡maxΘ∫zlnp(x,z;Θ)p(z|x;Θ(t))p(z|x;Θ(t))dz 展开对数项后，注意到 p(z|x;Θ(t)) 关于 Θ 是常数，上面的迭代公式可化为 Θ(t+1)=arg⁡maxΘB(Θ,Θ(t))=arg⁡maxΘ∫zlnp(x,z;Θ)p(z|x;Θ(t))dz=arg⁡maxΘQ(Θ,Θ(t)) 这就是 EM 算法最终的迭代公式，通常把最后的积分项记为 Q(Θ,Θ(t))，称为 Q 函数

• 经过这一系列的迭代，每一步选取的 Θ 都使下界 B(Θ,Θ(t)) 得到增大，迫使对数似然函数也被增大。再由单调有界原理可以证明 EM 算法将会收敛到一个局部极大值点。

下图是《统计学习方法》的插图，解释了 EM 算法的迭代过程

如果将积分项写成数学期望，则迭代公式还可写作 Θ(t+1)=arg⁡maxΘEz|x;Θt[lnp(x,z;Θ)] 可以看到这个迭代公式分两步，第一步求期望，第二步求最大化，这也是 EM 算法 (Expectation-Maximum) 名字的由来

总结一下，EM 算法分为 E 步（求期望步）和 M 步（最大化步）：

1. 使用当前步估计的参数 Θ(t)，求对数似然函数关于 z 的期望 Ez|x;Θt[lnp(x,z;Θ)]
2. 最大化对数似然函数的期望，求得参数 Θ(t+1) 作为下一步迭代的参数

EM 算法的收敛性

正如前文所言，EM 算法的迭代使对数似然的下界 B(Θ,Θ(t)) 不断增大，从而迫使对数似然函数增大。因此只需要证明 lnp(x;Θ(t)) 关于迭代次数 t 单调递增，再用单调有界原理即可证明算法收敛。

下面证明 lnp(x;Θ(t)) 单调递增： lnp(x;Θ)=lnp(x,z;Θ)−lnp(z|x;Θ) 两边对概率分布 p(z|x;Θ(t)) 取期望，注意到 lnp(x;Θ) 和 z 无关 lnp(x;Θ)=∫zlnp(x,z;Θ)p(z|x;Θ(t))dz−∫zlnp(z|x;Θ)p(z|x;Θ(t))dz 等式右边第一项就是 Q(Θ,Θ(t))，第二项我们记为 H(Θ,Θ(t))

作差比大小： lnp(x;Θ(t+1))−lnp(x;Θ(t))=[Q(Θ(t+1),Θ(t))−Q(Θ(t),Θ(t))]+[H(Θ(t),Θ(t))−H(Θ(t+1),Θ(t))] 因为 Θ(t+1)=arg⁡maxΘQ(Θ,Θ(t))，自然有 Q(Θ(t+1),Θ(t))≥Q(Θ(t),Θ(t))，所以等号右边第一项非负

再来看右边第二项 H(Θ(t),Θ(t))−H(Θ(t+1),Θ(t))=−∫zlnp(z|x;Θ(t+1))p(z|x;Θ(t))p(z|x;Θ(t))dz≥−ln∫zp(z|x;Θ(t+1))p(z|x;Θ(t))p(z|x;Θ(t))dz=−ln1=0 这说明第二项也是非负的，所以我们证明了 lnp(x;Θ(t+1))−lnp(x;Θ(t))≥0 再由单调有界原理可知，lnp(x;Θ(t)) 最终会收敛到某一值

更进一步可以证明 EM 算法得到的参数估计序列 Θ(t) 会收敛到某个稳定点 Θ^，这里不再深入。

• EM 算法容易受初值影响，只能收敛到局部极大值点而不能收敛到全局最大值点

高斯混合模型（GMM）

用高斯分布对数据集建模

混合高斯模型 (Gaussian Mixture Model) 使用 K 个高斯分布的加权和对数据集进行建模，认为观测数据分两步生成：

1. 由离散型随机变量 z 选择数据属于哪个成分，z 的概率分布是 P(z=k)=αk；
2. 假设选出的是第 k 个成分，在第 k 个高斯分布中采样得到一个观测数据，重复 n 次这样的操作得到整个数据集 x={x1,..,xn}

如果数据的生成过程看作上述两步的话，整个数据集服从的分布就是高斯分布的加权和 pM(x;Θ)=∑zp(x,z)=∑k=1Kp(x|z)P(z=k)=∑k=1Kαkϕ(x;μk,σk2) 其中 ϕ 是高斯分布的概率密度函数，Θ 是所有参数，即 Θ=(α,μ,σ2)。μk,σk2 是第 k 个高斯分布成分的均值和方差，即 ϕ(x;μk,σk2)=12πσke−(x−μk)22σk2

用 EM 算法进行参数估计

概率分布 pM(x;Θ) 中需要估计的参数是

α=(α1,...,αK),μ=(μ1,...,μK),σ2=(σ12,...,σK2)

由于随机变量 z 无法被观测到，只能通过 EM 算法进行求解：

E-Step

记概率分布的真实参数 Θ=(α,μ,σ2), 记算法第 t 步迭代时估计的参数近似值为 Θ(t)=(μ(t),σ(t),α(t))

先写出完整数据 (X,Z) 的似然函数 p(x,z;Θ)=∏k=1K∏i=1n[αkϕ(xi;μk,σk2)]zik=∏k=1K∏i=1n[αk12πσke−(xi−μk)22σk2]zik 这里 zik 是 0-1 随机变量，即示性函数 zik=I(zi=k)，表示第 i 个数据 xi 是否来自第 k 个高斯成分 zik={1,zi=k0,zi≠k 再写出完整数据 (X,Z) 的对数似然函数 lnp(x,z;Θ)=∑k=1K∑i=1nzik[lnαk−ln2π−lnσk−(xi−μk)22σk2] 再写出 Q 函数 Q(Θ,Θ(t))=Ez|x;Θ(t)[lnp(x,z;Θ)]=∑k=1K∑i=1nEz|x;Θ(t)[zik][lnαk−ln2π−lnσk−(xi−μk)22σk2] 为了方便以下叙述，这里先计算 Ez|x;Θ(t)[zik]Ez|x;Θ(t)[zik]=P(zi=k|xi;Θ(t))=P(zi=k,xi;Θ(t))p(xi;Θ(t))=P(zi=k,xi;Θ(t))∑k=1KP(xi,zi=k;Θ(t))=p(xi|zi=k;Θ(t))P(zi=k;Θ(t))∑k=1Kp(xi|zi=k;Θ(t))P(zi=k;Θ(t))=αk(t)ϕ(xi;Θk(t))∑k=1Kαk(t)ϕ(xi;Θk(t))

M-Step

为求解 Θ(t+1)=arg⁡maxΘQ(Θ,Θ(t))，只需分别令 Q 函数对 μk,σk2,αk 偏导数为 0

求导时需注意到 Ez|x;Θ(t)[zik] 的表达式与 μk,σk2,αk 都无关

先令 Q 对 μk 的偏导数为 0 ∂Q(Θ,Θ(t))∂μk=∑i=1nEz|x;Θ(t)[zik]xi−μkσk2=0 化简得到 μk(t+1)=∑i=1nEz|x;Θ(t)[zik]xi∑i=1nEz|x;Θ(t)[zik]

再令 Q 对 σk2 的偏导数为 0 ∂Q(Θ,Θ(t))∂σk2=∑i=1nEz|x;Θ(t)[zik](xi−μk)2−σk22σk4=0 化简得到 σk2(t+1)=∑i=1nEz|x;Θ(t)[zik](xi−μk)2∑i=1nEz|x;Θ(t)[zik] 最后求解 αk(t+1)，注意需要满足约束条件 ∑k=1Kαk=1，因此用拉格朗日乘子法求解 L=Q(Θ,Θ(t))+λ(∑k=1Kαk−1)

∂L∂αk=λ+1αk∑i=1nEz|x;Θ(t)[zik]=0

化简得 αk(t+1)=−∑i=1nEz|x;Θ(t)[zik]λ 加上约束条件 ∑k=1Kαk=1 可得 −∑k=1K∑i=1nEz|x;Θ(t)[zik]=λ 注意到 ∑k=1KEz|x;Θ(t)[zik]=1，所以 λ=−nαk(t+1)=∑i=1nEz|x;Θ(t)[zik]n 这样就得到了 GMM 最终的迭代公式 Ez|x;Θ(t)[zik]=αk(t)ϕ(xi;Θk(t))∑k=1Kαk(t)ϕ(xi;Θk(t))μk(t+1)=∑i=1nEz|x;Θ(t)[zik]xi∑i=1nEz|x;Θ(t)[zik]σk2(t+1)=∑i=1nEz|x;Θ(t)[zik](xi−μk)2∑i=1nEz|x;Θ(t)[zik]αk(t+1)=∑i=1nEz|x;Θ(t)[zik]n 模型训练结束后可用来进行聚类，因为 zik 是指示数据 xi 是否来自于第 k 个高斯成分的随机变量，可以使用 zik 的数学期望作为数据 xi 的类标记 λi=arg⁡maxkEz|x;Θ^[zik] 下面是西瓜书给出的高斯混合聚类算法伪代码，和本文记号有一点区别

周志华. 机器学习. 北京：清华大学出版社，2016. Print.

李航. 统计学习方法（第 2 版）. 清华大学出版社，2019. Print.

一篇推导比较详细的知乎

白板推导系列视频