本文从高斯混合模型出发,引出EM算法,最后回归到利用EM算法求解高斯混合模型。理论部分力求详尽不留证明疑点,所以略显冗长。实验部分给出了生成高斯混合分布样本和利用EM算法求解高斯混合模型的matlab代码。
理论部分
高斯混合模型(GMM)
顾名思义,高斯混合模型就是由多个高斯分布混合构成的模型。 K K K高斯混合分布的概率密度为:
p ( x ) = ∑ k = 1 K ϕ k N ( x ∣ μ k , Σ k ) . p(\mathbf{x})=\sum_{k=1}^K \phi_k\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k). p(x)=k=1∑KϕkN(x∣μk,Σk).
这里, ∑ k = 1 K ϕ k = 1 \sum_{k=1}^{K}\phi_k=1 ∑k=1Kϕk=1为混合系数,
N ( x ∣ μ , Σ ) = 1 ( 2 π ) D / 2 1 ∣ Σ ∣ 1 / 2 exp { − 1 2 ( x − μ ) T Σ − 1 ( x − μ ) } \mathcal{N}(\mathbf{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})=\frac{1}{(2\pi)^{D/2}}\frac{1}{|\boldsymbol{\Sigma}|^{1/2}}\exp\left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\} N(x∣μ,Σ)=(2π)D/21∣Σ∣1/21exp{
−21(x−μ)TΣ−1(x−μ)}
为 D D D维高斯分布,其中 μ \boldsymbol{\mu} μ为其均值向量, Σ \boldsymbol{\Sigma} Σ为其协方差矩阵。
直观来看,高斯混合分布可以看做下面分步过程的整合:
第一步,以 ϕ k \phi_k ϕk概率选择第 k k k个高斯模型;
第二步,利用第 k k k个高斯模型生成一个样本 x \mathbf{x} x。
为了讨论方便,记 θ = { ϕ 1 , … , ϕ K , μ 1 , … , μ K , Σ 1 , … , Σ K } \theta=\{\phi_1,\dots,\phi_K,\boldsymbol{\mu}_1,\dots,\boldsymbol{\mu}_K,\boldsymbol{\Sigma}_1,\dots,\boldsymbol{\Sigma}_K\} θ={
ϕ1,…,ϕK,μ1,…,μK,Σ1,…,ΣK}为高斯混合分布的参数集合。现在要解决的问题是,对于给定服从高斯混合分布的独立同分布样本集 X = { x 1 , … , x n } \mathbf{X}=\{\mathbf{x}_1, \dots, \mathbf{x}_n\} X={
x1,…,xn},最大化其对数似然函数:
max θ ln p ( X ∣ θ ) = ∑ i = 1 n ln ∑ k = 1 K ϕ k N ( x i ∣ μ k , Σ k ) . \max_{\theta}\ln p(\mathbf{X}|\theta)=\sum_{i=1}^n\ln\sum_{k=1}^K\phi_k\mathcal{N}(\mathbf{x}_i|\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k). θmaxlnp(X∣θ)=i=1∑nlnk=1∑KϕkN(xi∣μk,Σk).
这里, p ( X ∣ θ ) = ∏ i = 1 n p ( x i ∣ θ ) p(\mathbf{X}|\theta)=\prod\limits_{i=1}^np(\mathbf{x}_i|\theta) p(X∣θ)=i=1∏np(xi∣θ)。由于 ln \ln ln函数里的求和项,我们无法直接求得闭式解,而利用EM算法可以得到一个局部最优数值解。
EM算法
从分步角度来看,求解高斯混合模型的难点在于,我们不知道一个样本 x i \mathbf{x}_i xi具体是由 K K K个高斯模型中的哪一个生成的。所以,对于第 i i i个样本 x i \mathbf{x}_i xi来说,我们构造一个隐变量 z i z_i zi用来表示 x i \mathbf{x}_i xi来自于哪个高斯模型。也即, z i = k z_i=k zi=k当且仅当 x i \mathbf{x}_i xi来自于第 k k k个高斯模型。注意,虽然这个隐变量的取值是客观确定的,但对我们来说是不可见,因此仍将其看作随机变量。记 Z = { z 1 , … , z n } \mathbf{Z}=\{z_1,\dots,z_n\} Z={
z1,…,zn},对于固定的模型参数 θ ˉ \bar{\theta} θˉ,下面的不等式给出了EM算法的框架。
ln p ( X ∣ θ ˉ ) = ∑ Z p ( Z ∣ X , θ ˉ ) ln p ( X ∣ θ ˉ ) = ∑ Z p ( Z ∣ X , θ ˉ ) ln p ( X , Z ∣ θ ˉ ) p ( Z ∣ X , θ ˉ ) ≤ max θ ∑ Z p ( Z ∣ X , θ ˉ ) ln p ( X , Z ∣ θ ) p ( Z ∣ X , θ ˉ ) ( 1 ) ≤ ln ∑ Z p ( X , Z ∣ θ m a x ) p ( Z ∣ X , θ ˉ ) p ( Z ∣ X , θ ˉ ) ( J e n s e n 不 等 式 ) = ln p ( X ∣ θ m a x ) \begin{aligned} \ln p(\mathbf{X}|\bar{\theta})&=\sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\bar{\theta})\ln p(\mathbf{X}|\bar{\theta})\\ &=\sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\bar{\theta})\ln\frac{p(\mathbf{X},\mathbf{Z}|\bar{\theta})}{p(\mathbf{Z}|\mathbf{X},\bar{\theta})} \\ &\leq\max_\theta \sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\bar{\theta})\ln\frac{p(\mathbf{X},\mathbf{Z}|\theta)}{p(\mathbf{Z}|\mathbf{X},\bar{\theta})} & (1)\\ &\leq \ln \sum_{\mathbf{Z}}\frac{p(\mathbf{X},\mathbf{Z}|\theta_{max})}{p(\mathbf{Z}|\mathbf{X},\bar{\theta})}p(\mathbf{Z}|\mathbf{X},\bar{\theta}) & (Jensen不等式)\\ &=\ln p(\mathbf{X}|\theta_{max}) \end{aligned} lnp(X∣θˉ)=Z∑p(Z∣X,θˉ)lnp(X∣θˉ)=Z∑p(Z∣X,θˉ)lnp(Z∣X,θˉ)p(X,Z∣θˉ)≤θmaxZ∑p(Z∣X,θˉ)lnp(Z∣X,θˉ)p(X,Z∣