<link href="https://****img.cn/public/favicon.ico" rel="SHORTCUT ICON">
<title>详解EM算法与混合高斯模型(Gaussian mixture model, GMM) - 林立民爱洗澡 - ****博客</title>
最近在看晓川老(shi)师(shu)的博士论文,接触了混合高斯模型(Gaussian mixture model, GMM)和EM(Expectation Maximization)算法,不禁被论文中庞大的数学公式所吓退。本文通过查阅相关资料,在复杂巧妙的推理公式中融入了自己的理解,详细梳理了混合高斯模型和EM算法。
1 单高斯模型(Gaussian single model, GSM)
简单回顾一下概率论讲过的高斯模型。
高斯模型是一种常用的变量分布模型,在数理统计领域有着广泛的应用(……好吧读了这么多年书没白费,教科书般的话语已植入骨髓)。一维高斯分布的概率密度函数如下:
(1)f(x)=12πσexp(−(x−μ)22σ2)f(x)=2πσ1exp(−2σ2(x−μ)2)f(x)=2πσ1exp(−2σ2(x−μ)2)(1)(1)
μμμ和σ2σ2σ2分别是高斯分布的均值和方差。
譬如将男生身高视为变量X, 假设男生的身高服从高斯分布,则X∼N(μ,σ2)X∼N(μ,σ2)X∼N(μ,σ2),女生亦如此。只是男女生身高分布可能具有不同的均值和方差。图1是从谷歌图库中搜索到的男女生身高分布图,来源不清,个人觉得男生的均值身高虚高……四个记号分别表示0~3σσσ准则。

图1 男女生身高分布差异
多维变量X=(x1,x2,...xn)X=(x1,x2,...xn)X=(x1,x2,...xn)的联合概率密度函数为:
(2)f(X)=1(2π)d/2∣Σ∣1/2exp[−12(X−u)TΣ−1(X−u)],X=(x1,x2...xn)f(X)=(2π)d/2∣Σ∣1/21exp[−21(X−u)TΣ−1(X−u)],X=(x1,x2...xn)f(X)=(2π)d/2∣Σ∣1/21exp[−21(X−u)TΣ−1(X−u)],X=(x1,x2...xn)(2)(2)
其中:
d:变量维度。对于二维高斯分布,有d=2;
u=(u1u2...un)u=⎝⎜⎜⎛u1u2...un⎠⎟⎟⎞u=⎝⎜⎜⎛u1u2...un⎠⎟⎟⎞:各维变量的均值;
ΣΣΣ:协方差矩阵,描述各维变量之间的相关度。对于二维高斯分布,有:
(3)Σ=[δ11δ12δ21δ22]Σ=[δ11δ21amp;δ12amp;δ22]Σ=[δ11δ21δ12δ22](3)(3)

图2 二维高斯数据分布
图2是二维高斯分布产生的数据示例,参数设定为:u=(00),Σ=[10.80.85]u=(00),Σ=[10.8amp;0.8amp;5]u=(00),Σ=[10.80.85]。关于二维高斯分布的参数设定对为高斯曲面的影响,这篇文章二维高斯分布(Two-dimensional Gaussian distribution)的参数分析有提及,主要是为了下文理解混合高斯分布做铺垫。服从二维高斯分布的数据主要集中在一个椭圆内部,服从三维的数据集中在一个椭球内部。
2 混合高斯模型(Gaussian mixture model, GMM)
2.1 为什么要有混合高斯模型
先来看一组数据。

图3 混合高斯分布所产生数据
如果我们假设这组数据是由某个高斯分布产生的,利用极大似然估计(后文还会提及)对这个高斯分布做参数估计,得到一个最佳的高斯分布模型如下。

图4 用单高斯模型对样本作分析不合理示意
有什么问题吗?一般来讲越靠近椭圆的中心样本出现的概率越大,这是由概率密度函数决定的,但是这个高斯分布的椭圆中心的样本量却极少。显然样本服从单高斯分布的假设并不合理。单高斯模型无法产生这样的样本。
实际上,这是用两个不同的高斯分布模型产生的数据。

图5 混合高斯模型对样本作分析示意
正当单高斯模型抓耳挠腮的时候,混合高斯模型就大摇大摆地进场了。它通过求解两个高斯模型,并通过一定的权重将两个高斯模型融合成一个模型,即最终的混合高斯模型。这个混合高斯模型可以产生这样的样本。
更一般化的描述为:假设混合高斯模型由K个高斯模型组成(即数据包含K个类),则GMM的概率密度函数如下:
(3)p(x)=∑k=1Kp(k)p(x∣k)=∑k=1KπkN(x∣uk,Σk)p(x)=k=1∑Kp(k)p(x∣k)=k=1∑KπkN(x∣uk,Σk)p(x)=k=1∑Kp(k)p(x∣k)=k=1∑KπkN(x∣uk,Σk)(3)(3)
其中,p(x∣k)=N(x∣uk,Σk)p(x∣k)=N(x∣uk,Σk)p(x∣k)=N(x∣uk,Σk)是第k个高斯模型的概率密度函数,可以看成选定第k个模型后,该模型产生x的概率;p(k)=πkp(k)=πkp(k)=πk是第k个高斯模型的权重,称作选择第k个模型的先验概率,且满足∑k=1Kπk=1k=1∑Kπk=1k=1∑Kπk=1。
所以,混合高斯模型并不是什么新奇的东西,它的本质就是融合几个单高斯模型,来使得模型更加复杂,从而产生更复杂的样本。理论上,如果某个混合高斯模型融合的高斯模型个数足够多,它们之间的权重设定得足够合理,这个混合模型可以拟合任意分布的样本。
2.2 直观上理解混合高斯模型
下面通过几张图片来帮助理解混合高斯模型。
首先从简单的一维混合高斯模型说起。

图6 一维混合高斯模型
在图6中,y1,y2和y3分别表示三个一维高斯模型,他们的参数设定如图所示。y4表示将三个模型的概率密度函数直接相加,注意的是这并不是一个混合高斯模型,因为不满足∑k=1Kπk=1k=1∑Kπk=1k=1∑Kπk=1的条件。而y5和y6分别是由三个相同的高斯模型融合生成的不同混合模型。由此可见,调整权重将极大影响混合模型的概率密度函数曲线。另一方面也可以直观地理解混合高斯模型可以更好地拟合样本的原因:它有更复杂更多变的概率密度函数曲线。理论上,混合高斯模型的概率密度函数曲线可以是任意形状的非线性函数。
下面再给出一个二维空间3个高斯模型混合的例子。

(a) 3个类别高斯分布截面轮廓线

(b) 混合高斯分布截面轮廓线

© 二维混合高斯分布概率密度函数图
图7 二维混合高斯模型
(a) 图表示的是3个高斯模型的截面轮廓图,3个模型的权重系数已在图中注明,由截面轮廓图可知3个模型之间存在很多重叠区域。其实这也正是混合高斯模型所希望的。因为如果它们之间的重叠区域较少,那么生成的混合高斯模型一般较为简单,难以生成较为复杂的样本。
设定好了3个高斯模型和它们之间的权重系数之后,就可以确定二维混合高斯分布的概率密度函数曲面,如图©所示。图(b)是对于图©概率密度曲面的截面轮廓线。从图7也可以看出,整个混合高斯分布曲面相对比于单高斯分布曲面已经异常复杂。实际上,通过调整混合高斯分布的系数(π,μ,Σ)(π,μ,Σ)(π,μ,Σ),可以使得图©的概率密度曲面去拟合任意的三维曲面,从而采样生成所需要的数据样本。
3 极大似然估计(Maximum Likehood Estimate, MLE)(最大化对数似然函数)
● 最大化对数似然函数(log-likelihood function)的意义
首先直观化地解释一下最大化对数似然函数要解决的是什么问题。
假设我们采样得到一组样本ytytyt,而且我们知道变量Y服从高斯分布(本文只提及高斯分布,其他变量分布模型类似),数学形式表示为Y∼N(μ,Σ)Y∼N(μ,Σ)Y∼N(μ,Σ)。采样的样本如图8所示,我们的目的就是找到一个合适的高斯分布(也就是确定高斯分布的参数μ,Σμ,Σμ,Σ),使得这个高斯分布能产生这组样本的可能性尽可能大。

图8 最大化似然函数的意义
那怎么找到这个合适的高斯分布呢(在图8中的表示就是1~4哪个分布较为合适)?这时候似然函数就闪亮登场了。
似然函数数学化:设有样本集Y=y1,y2...yNY=y1,y2...yNY=y1,y2...yN。p(yn∣μ,Σ)p(yn∣μ,Σ)p(yn∣μ,Σ)是高斯分布的概率分布函数,表示变量Y=ynY=ynY=yn的概率。假设样本的抽样是独立的,那么我们同时抽到这N个样本的概率是抽到每个样本概率的乘积,也就是样本集Y的联合概率。此联合概率即为似然函数:
(4)L(μ,Σ)=L(y1,y2...yN;μ,Σ)=∏n=1Np(yn;μ,Σ)L(μ,Σ)=L(y1,y2...yN;μ,Σ)=n=1∏Np(yn;μ,Σ)L(μ,Σ)=L(y1,y2...yN;μ,Σ)=n=1∏Np(yn;μ,Σ)(4)(4)
对式子(4)进行求导并令导数为0(即最大化似然函数,一般还会先转化为对数似然函数再最大化),所求出的参数就是最佳的高斯分布对应的参数。
所以最大化似然函数的意义就是:通过使得样本集的联合概率最大来对参数进行估计,从而选择最佳的分布模型。
对于图8产生的样本用最大化似然函数的方法,最终可以得到序号1对应的高斯分布模型是最佳的模型。
4 EM算法(最大化Q函数)
4.1 为什么要有EM算法(EM算法与极大似然估计分别适用于什么问题)
● 尝试用极大似然估计的方法来解GMM模型
解GMM模型,实际上就是确定GMM模型的参数(μ,Σ,π)(μ,Σ,π)(μ,Σ,π),使得由这组参数确定的GMM模型最有可能产生采样的样本。
先试试看用极大似然估计的方法来解GMM模型会出现什么样的问题。
如第3小节所述,要利用极大似然估计求解模型最重要的一步就是求出似然函数,即样本集出现的联合概率。而对于混合高斯模型,如何求解某个样本ytytyt的概率?显然我们得先知道这个样本来源于哪一类高斯模型,然后求这个高斯模型生成这个样本的概率p(yt)p(yt)p(yt)。
但是问题来了:我们只有样本。不知道样本到底来源于哪一类的高斯模型。那么如何求解样本的生成概率p(yt)p(yt)p(yt)?
先引入一个隐变量γγγ。它是一个K维二值随机变量,在它的K维取值中只有某个特定的元素γkγkγk的取值为1,其它元素的取值为0。实际上,隐变量描述的就是:每一次采样,选择第k个高斯模型的概率,故有:
(5)p(γk=1)=πkp(γk=1)=πkp(γk=1)=πk(5)(5)
当给定了γγγ的一个特定的值之后(也就是知道了这个样本从哪一个高斯模型进行采样),可以得到样本y的条件分布是一个高斯分布,满足:
(6)p(y∣γk=1)=N(y∣μk,Σk)p(y∣γk=1)=N(y∣μk,Σk)p(y∣γk=1)=N(y∣μk,Σk)(6)(6)
而实际上,每个样本到底是从这K个高斯模型中哪个模型进行采样的,是都有可能的。故样本y的概率为:
(7)p(y)=∑γp(γ)p(y∣γ)=∑k=1KπkN(y∣μk,Σk)p(y)=∑γp(γ)p(y∣γ)=k=1∑KπkN(y∣μk,Σk)p(y)=∑γp(γ)p(y∣γ)=k=1∑KπkN(y∣μk,Σk)(7)(7)
样本集Y(n个样本点)的联合概率为:
(8)L(μ,Σ,π)=L(y1,y2...yN;μ,Σ,π)=∏n=1Np(yn;μ,Σ,π)=∏n=1N∑k=1KπkN(yn∣μk,Σk)L(μ,Σ,π)=L(y1,y2...yN;μ,Σ,π)=n=1∏Np(yn;μ,Σ,π)=n=1∏Nk=1∑KπkN(yn∣μk,Σk)L(μ,Σ,π)=L(y1,y2...yN;μ,Σ,π)=n=1∏Np(yn;μ,Σ,π)=n=1∏Nk=1∑KπkN(yn∣μk,Σk)(8)(8)
对数似然函数表示为:
(9)lnL(μ,Σ,π)=∑n=1Nln∑k=1KπkN(yn∣μk,Σk)lnL(μ,Σ,π)=n=1∑Nlnk=1∑KπkN(yn∣μk,Σk)lnL(μ,Σ,π)=n=1∑Nlnk=1∑KπkN(yn∣μk,Σk)(9)(9)
好了,然后求导,令导数为0,得到模型参数(μ,Σ,π)(μ,Σ,π)(μ,Σ,π)。
貌似问题已经解决了,喜大普奔。
然而仔细观察可以发现,对数似然函数里面,对数里面还有求和。实际上没有办法通过求导的方法来求这个对数似然函数的最大值。
MLE(极大似然估计)略显沮丧。这时候EM算法走过来,安慰着说:兄弟别哭,老哥帮你。
● 极大似然估计与EM算法适用问题分析
下面先阐述一下极大似然估计与EM算法分别适用于解决什么样的问题。

图9 极大似然估计适用问题

图10 EM算法适用问题
如果我们已经清楚了某个变量服从的高斯分布,而且通过采样得到了这个变量的样本数据,想求高斯分布的参数,这时候极大似然估计可以胜任这个任务;而如果我们要求解的是一个混合模型,只知道混合模型中各个类的分布模型(譬如都是高斯分布)和对应的采样数据,而不知道这些采样数据分别来源于哪一类(隐变量),那这时候就可以借鉴EM算法。EM算法可以用于解决数据缺失的参数估计问题(隐变量的存在实际上就是数据缺失问题,缺失了各个样本来源于哪一类的记录)。
下面将介绍EM算法的两个步骤:E-step(expectation-step,期望步)和M-step(Maximization-step,最大化步);
4.2 E-step
我们现有样本集Y=(y1,y2...yT)Y=(y1,y2...yT)Y=(y1,y2...yT),通过隐变量γt,kγt,kγt,k(表示ytytyt这个样本来源于第k个模型)的引入,可以将数据展开成完全数据:
(yt,γt,1,γt,2...γt,K),t=1,2...T(yt,γt,1,γt,2...γt,K),t=1,2...T(yt,γt,1,γt,2...γt,K),t=1,2...T
所谓的完全数据,就是不缺失的数据。只有样本集Y=(y1,y2...yT)Y=(y1,y2...yT)Y=(y1,y2...yT)的数据是不完整的,存在信息缺失的。若ytytyt由第1类采样而来,则有γt,1=1,γt,2=0...γt,K=0γt,1=1,γt,2=0...γt,K=0γt,1=1,γt,2=0...γt,K=0,表示为(yt,1,0,...0)(yt,1,0,...0)(yt,1,0,...0)。
所以要求能采到这组数据的可能性,需要分两步走:①第t个样本由哪一类产生?②如果第t个样本由第k类产生,那么第k类产生第t个样本的概率为多少?
综合考虑上面两步,有了完全数据的似然函数:
(10)p(y,γ∣μ,Σ,π)=∏t=1Tp(yt,γt,1,γt,2...γt,K∣μ,Σ,π) =∏t=1T∏k=1K(πkN(yt;μk,Σk))γt,k =∏k=1Kπk∑t=1Tγt,k∏t=1T(N(yt;μk,Σk))γt,kp(y,γ∣μ,Σ,π)=t=1∏Tp(yt,γt,1,γt,2...γt,K∣μ,Σ,π) =t=1∏Tk=1∏K(πkN(yt;μk,Σk))γt,k =k=1∏Kπk∑t=1Tγt,kt=1∏T(N(yt;μk,Σk))γt,kp(y,γ∣μ,Σ,π)=t=1∏Tp(yt,γt,1,γt,2...γt,K∣μ,Σ,π) =t=1∏Tk=1∏K(πkN(yt;μk,Σk))γt,k =k=1∏Kπk∑t=1Tγt,kt=1∏T(N(yt;μk,Σk))γt,k(10)(10)
第1个等号到第2个等号的理解:若ytytyt由第1类采样而来,则有γt,1=1,γt,2=0...γt,K=0γt,1=1,γt,2=0...γt,K=0γt,1=1,γt,2=0...γt,K=0,
(11)p(yt,γt,1,γt,2...γt,K∣μ,Σ,π)=∏k=1K(πkN(yt;μk,Σk))γt,k =(π1N(yt;μ1,Σ1))γt,1(π2N(yt;μ2,Σ2))γt,2...(πKN(yt;μK,ΣK))γt,K =(π1N(yt;μ1,Σ1))1(π2N(yt;μ2,Σ2))0...(πKN(yt;μK,ΣK))0 =π1N(yt;μ1,Σ1)p(yt,γt,1,γt,2...γt,K∣μ,Σ,π)=k=1∏K(πkN(yt;μk,Σk))γt,k =(π1N(yt;μ1,Σ1))γt,1(π2N(yt;μ2,Σ2))γt,2...(πKN(yt;μK,ΣK))γt,K =(π1N(yt;μ1,Σ1))1(π2N(yt;μ2,Σ2))0...(πKN(yt;μK,ΣK))0 =π1N(yt;μ1,Σ1)p(yt,γt,1,γt,2...γt,K∣μ,Σ,π)=k=1∏K(πkN(yt;μk,Σk))γt,k =(π1N(yt;μ1,Σ1))γt,1(π2N(yt;μ2,Σ2))γt,2...(πKN(yt;μK,ΣK))γt,K =(π1N(yt;μ1,Σ1))1(π2N(yt;μ2,Σ2))0...(πKN(yt;μK,ΣK))0 =π1N(yt;μ1,Σ1)(11)(11)
注意式子(11)与式子(7)的差别:如果求p(yt)p(yt)p(yt)则需要考虑ytytyt有可能来源于k个类;而如果求的是p(yt,γt,1,γt,2...γt,K)p(yt,γt,1,γt,2...γt,K)p(yt,γt,1,γt,2...γt,K)则已经限定了ytytyt只会来源于某个类。
第2个等式到第3个等式的理解:先交换累乘符号。由于πkπkπk与t无关,故而可以从内部的累乘符号中提取出来。
实际上完全数据的似然函数描述的就是采集到这些样本的可能性。
完全数据的对数似然函数为:
(12)lnp(y,γ∣μ,Σ,π)=∑k=1K(∑t=1Tγt,k)lnπk+∑t=1Tγt,k(−ln(2π)−12ln∣Σk∣−12(yt−μt)T(Σk)−1(yt−μt))lnp(y,γ∣μ,Σ,π)=k=1∑K(t=1∑Tγt,k)lnπk+t=1∑Tγt,k(−ln(2π)−21ln∣Σk∣−21(yt−μt)T(Σk)−1(yt−μt))lnp(y,γ∣μ,Σ,π)=k=1∑K(t=1∑Tγt,k)lnπk+t=1∑Tγt,k(−ln(2π)−21ln∣Σk∣−21(yt−μt)(12)