Loading [MathJax]/jax/output/HTML-CSS/jax.js

上一篇中我们着重介绍了对于离散数据的生成模型,紧接上一篇,本篇我们介绍对于连续数据的生成模型。好吧,废话我们就不多说了,直接进入正文。

MLE for MVN


Basics about MVN

谈到连续分布,我们很自然地就会想到高斯分布,从小学到现在,印象中第一个走入我脑海中的看着比较高端大气上档次的就是Gaussian分布了。这次我们的重点也会完全集中在Gaussian分布之了,在正式讨论之前,我们先介绍一些关于Gaussian分布的基础知识。

D维空间中,MVN(Multivariate Normal)多变量正态分布的概率分布函数具有如下形式:

N(x|μ,Σ)1(2π)D/2det(Σ)1/2exp[12(xμ)TΣ1(xμ)]

上式中的指数部分是xμ之间的Mahalanobis距离。为了更好地理解这个量,我们对Σ做特征值分解,即Σ=UΛUT,其中U是一正交阵,满足UTU=I,而Λ是特征值矩阵。

通过特征值分解,我们有:

Σ1=UTΛ1U1=UΛ1UT=Di=11λiuiuTi

其中,uiU的第i列。因此Mahalanobis距离可被改写为:

(xμ)TΣ1(xμ)=(xμ)T(Di=11λiuiuTi)(xμ)=Di=11λi(xμ)TuiuTi(xμ)=Di=1y2ixi

其中,yiuTi(xμ)。另2维空间中的椭圆方程为:

y21λ1+y22λ2=1

因此我们可知Gaussian概率密度的等高线沿着椭圆分布,特征向量决定椭圆的朝向,而特征值则决定椭圆有多“椭”。一般来说,如果我们将坐标系移动μ,然后按U旋转,此时的欧拉距离即为Mahalanobis距离。

MLE for MVN

以下我们给出MVN参数的MLE(极大似然估计)的证明:

Theorem 3.1 若我们获取的N个独立同分布的样本xi N(x|μ,Σ),则关于μ以及 Σ的极大似然分布如下: MLE for Gaussian

我们不加证明地给出如下公式组:

(bTa)a=b(aTAa)a=(A+AT)aAtr(BA)=BTAlog|A|=ATtr(ABC)=tr(CAB)=tr(BCA)

最后一个等式称为迹的循环置换性质(cyclic permutation property)。利用这个性质,我们使用trace trick可以得到下式:

xTAx=tr(xTAx)=tr(xxTA)=tr(AxxT)

证明:

对数似然函数为:

l(μ,Σ)=logp(D|μ,Σ)=N2log|Λ|12Ni=1(xiμ)TΛ(xiμ)

其中,Λ=Σ1为精度矩阵。令yi=xiμ并利用链式法则有:

μ(xiμ)TΣ1(xiμ)=yiyTiΣ1yiyiμ=(ΣT+Σ1)yi

即:

μl(μ,Σ)=12Ni=12Σ1(xiμ)=0

故有:

ˆμ=1NNi=1xi=ˉx

即最大似然均值即为经验均值。

利用trace trick我们重写对数似然函数为:

l(Λ)=N2log|Λ|12itr[(xiμ)(xiμ)TΛ]=N2log|Λ|12tr[SuΛ]

其中,SuNi=1(xiμ)(xiμ)T,业界尊称其为分散度矩阵(Scatter Matrix),以后我们聊LDA的时候会再次碰到。对Λ求偏导有:

l(Λ)Λ=N2ΛT12STu=0

得:

ˆΣ=1NNi=1(xiμ)(xiμ)T

证毕。

Gaussian Discriminant Analysis


我们在上一篇中提到了Naive Bayes方法,其实质无非是估计在每一类下特定的样本出现的概率,进而我们可以把该特定样本分配给概率值最大的那个类。而对于连续数据而言,其实质其实也是一样的,每一个MVN(我们可以看做一类或者一个Component)都可能生成一些数据,我们估计在每一个Component下生成特定样本的概率,然后把该特定样本分配给概率值最大的那个Component即可。即我们可以定义如下的条件分布:

p(x|y=c,θ)=N(x|μc,Σc)

上述模型即为高斯判别分析(Gaussian Discriminant Analysis,GDA)(注意,该模型为生成模型,而不是判别模型)。如果Σc是对角阵,即所有的特征都是独立的时,该模型等同于Naive Bayes.

QDA

在上式中带入高斯密度函数的定义,则有:

p(y=c|x,θ)=πc|2πΣc|12exp[12(xμc)TΣ1C(xμc)]cπc|2πΣc|12exp[12(xμc)TΣ1c(xμc)]

上式中π为各个Component的先验概率分布。根据上式得到的模型则称为Quadratic Discriminant Analysis(QDA).以下给出在2类以及3类情形下可能的决策边界形状,如下图所示:

Decision Boundary

Linear Discriminant Analysis(LDA)

当各个Gaussian Component的协方差矩阵相同时,此时我们有:

p(y=c|x,θ)πcexp[μTcΣ1x12xTΣ1x12μTcΣ1μc]=exp[μTcΣ1x12μTcΣ1μc+logπc]exp[12xTΣ1x]

上式中exp[12xTΣ1x]是独立于c的,分子分母相除抵消到此项。

令:

γc=12μTcΣ1μc+logπcβc=Σ1μc

于是有:

p(y=c|x,θ)=eβTcx+γcceβTcx+γc=S(η)c

其中η=[βT1x+γ1,...,βTCx+γc],S为softmax函数(类似于max函数,故得此名),定义如下:

Softmax Function

若将ηc除以一个常数(temperature),当T0时,我们有:

Softmax Function

换句话说,当温度很低时,分布集中在概率最大的那个状态上,而当温度高时,所有的状态呈现均匀分布。

NOTE:该术语来自于统计物理学,在统计物理学中,人们更倾向于使用波尔兹曼分布(Boltzmann distribution)一词。

关于式16有一个有趣的性质,即对该式取log,我们则会得到一个关于x的线性方程。因此对于任意两类之间的决策边界将会是一条直线,据此该模型也被称为线性判别分析(Linear Discriminant Analysis,LDA)。而且对于二分类问题,我们可以得到:

p(y=c|x,θ)=p(y=c|x,θ)

βTcx+γc=βTcx+γc

xT(βcβc)=γcγc

Two-class LDA

为了加深对以上等式的理解,对于二分类的情况,我们做如下说明:

p(y=1|x,θ)=eβT1x+γ1eβT1x+γ1+eβT0x+γ0=11+e(β0β1)Tx+(γ0γ1)=sigm((β1β0)Tx+(γ1γ0))

其中,sigm(η)代表Sigmoid函数

现有:

γ1γ0=12μT1Σ1μ1+12μT0Σ1μ0+log(π1/π0)=12(μ1μ0)TΣ1(μ1+μ0)+log(π1/π0)

因此若我们另:

ω=β1β0=Σ1(μ1μ0)x0=12(μ1+μ0)(μ1μ0)log(π1/π0)(μ1μ0)TΣ1(μ1μ0)

则有ωTx0=(γ1γ0),即:

p(y=1|x,θ)=sigm(ωT(xx0))

因此最后的决策规则很简单:将x平移x0,然后投影到ω上,通过结果是正还是负决定它到底属于哪一类。

Two class LDA

Σ=σ2I时,ωμ1μ0同向。这时我们只需要判断投影点离μ1μ0中的那个点近。当它们的先验概率π1=π0时,投影点位于其中点;当π1>π0时,则x0越趋近于μ0,直线的更大部分先验地属于类1;反之亦然。1

Inference in joint Gaussian distributions


给定联合概率分布p(x1,x2),如果我们能够计算边际概率分布p(x1)以及条件概率分布p(x1|x2)想必是极好的而且是及有用的。以下我们仅给出结论,下式表明如果两变量符合联合高斯分布,则它们的边际分布以及条件分布也都是高斯分布

Theorem 3.22 假定x=(x1,x2)服从联合高斯分布,且参数如下: μ=(μ1μ2) Σ=(Σ11Σ12Σ21Σ22) 则我们可以得到如下边际概率分布: p(x1)=N(x1|μ1,Σ11)p(x2)=N(x2|μ2,Σ22) 另其后验条件分布为: p(x1|x2)=N(x1|μ1|2,Σ1|2) μ1|2=μ1+Σ12Σ122(x2μ2)=μ1Λ111Λ12(x2μ2)=Σ1|2(Λ11μ1Λ12(x2μ2)) Σ1|2=Σ11Σ12Σ122Σ21=Λ111

Linear Gaussian Systems


给定两变量,xy.令xRDx为一隐含变量,yRDy为关于x的包含噪声的观察值。此外,我们假定存在如下prior和likelihood: p(x)=N(x|μx,Σx)p(y|x)=N(y|Ax+b,Σy) 上式即称为Linear Gaussian System。此时我们有:

Theorem 3.33 给定一Linear Gaussian System.其后验分布p(x|y)具有如下形式: p(x|y)=N(x|μx|y,Σx|y)Σ1x|y=Σ1x+ATΣ1yAμx|y=Σx|y[ATΣ1y(yb)+Σ1xμx]

Inferring an unknown vector from noisy measurements

下面我们举一个简单的例子以进一步说明Linear Gaussian System: 现有N个观测向量,yi N(x,Σy),prior服从高斯分布x N(μ0,Σ0).令A=I,b=0,此外,我们采用ˉy作为我们的有效估计值,其精度为NΣ1y,我们有: p(x|y1,...,yN)=N(x|μN,ΣN)Σ1N=Σ10+NΣ1yμN=ΣN(Σ1y(Nˉy)+Σ10μ0) 为了更直观地解释以上模型,如下图所示:

Radar Blips

我们可将x视为2维空间中一个物体的真实位置(但我们并不知道),例如一枚导弹或者一架飞机,yi则是我们的观测值(含噪声),可以视为雷达上的一些点。当我们得到越来越多的点时,我们就能够更好地进行定位。4

现假定我们有多个测量设备,且我们想利用多个设备的观测值进行估计,这种方法称为sensor fusion.如果我们有具有不同方差的多组观测值,那么posterior将会是它们的加权平均。如下图所示:

Sensor Fusion

我们采用不带任何信息的关于x的先验分布,即p(x)=N(μ0,Σ0)=N(0,1010I2),我们得到2个观测值,y1 N(x,Σy,1,y2 N(x,Σy,2),我们需要计算p(x|y1,y2).

如上图(a),我们设定Σy,1=Σy,2=0.01I2,因此两个传感器都相当可靠,posterior mean即位于两个观测值中间;如上图(b),我们设定Σy,1=0.05I2Σy,2=0.01I2,因此传感器2比传感器1可靠,此时posterior mean更靠近于y2;如上图(c),我们有: Σy,1=0.01(10111),Σy,2=0.01(11110)

从上式我们不难看出,传感器1在第2个分量上更可靠,传感器2在第1个分量上也更可靠。此时,posterior mean采用传感器1的第二分量以及传感器1的第二分量。

NOTE:当sensor测量精度未知时,我们则需要它们的测量精度也进行估计。

The Wishart Distribution


在多变量统计学中,Wishart分布是继高斯分布后最重要且最有用的模型。 ------Press

既然Press他老人家都说了Wishart分布很重要,而且我们下一部分会用到它,那么我们就必须得介绍介绍它了。(它主要被用来Model关于协方差矩阵的不确定性)

Wishart Distribution

Wishart概率密度函数具有如下形式: Wi(Λ|S,ν)=1ZWi|Λ|(νD1)/2exp(12tr(ΛS1))

其中,ν为自由度,S为缩放矩阵。其归一项具有如下形式:(很恐怖,对吧!) ZWi=2νD/2ΓD(ν/2)|S|ν/2

其中,ΓD(a)为多变量Gamma函数: ΓD(x)=πD(D1)/4Di=1Γ(x+(1i)/2)

于是有Γ1(a)=Γ(a)且有: ΓD(ν0/2)=Di=1Γ(ν0+1i2)

仅当ν>D1时归一项存在

其实Wishart分布和Gaussian分布是有联系的。具体而言,令xi N(0,Σ),则离散度矩阵S=Ni=1xixTi服从Wishart分布,且S Wi(Σ,1)。于是有E(S)=NΣ.

更一般地,我们可以证明Wi(S,ν)的mean和mode具有如下形式: mean=νS,mode=(νD1)S

仅当ν>D+1时mode存在

D=1时,Wishart分布退化为Gamma分布,且有: Wi(λ|s1,v)=Ga(λ|ν2,s2)

Inverse Wishart Distribution

Σ1 Wi(S,ν),则Σ IW(S1,ν+D+1),其中IW为逆Wishart分布。当ν>D1S0时,我们有: IW(Σ|S,ν)=1ZIW|Σ|(ν+D+1)/2exp(12tr(S1Σ1))ZIW=|S|ν/22νD/2ΓD(ν/2)

此外我们可以证明逆Wishart分布具有如下性质: mean=S1νD1,mode=S1ν+D+1

D=1时,它们退化为逆Gamma分布: IW(σ2|S1,ν)=IG(σ2|ν/2,S/2)

Inferring the parameters of an MVN5


到目前为止,我们已经讨论了在θ=(μ,Σ)已知的条件下如何inference,现我们讨论一下如何对参数本身进行估计。假定xi N(μ,Σ) for i=1:N.本节主要分为两个部分:

  • 计算p(μ|D,Σ);
  • 计算p(Σ|D,μ).

Posterior distribution of μ

我们已经就如何计算μ的极大似然估计值进行了讨论,现我们讨论如何计算其posterior.

其likelihood具有如下形式: p(D|μ)=N(ˉx|μ,1NΣ)

为了简便起见,我们采用共轭先验分布,即高斯。特别地,若p(μ)=N(μ|m0,V0).此时我们可以根据之前Linear Gaussian System的结论得到(和我们之前提到的雷达的例子雷同): p(μ|D,Σ)=N(μ|mN,VN)V1N=V10+NΣ1mN=VN(Σ1(Nˉx)+V10m0)

我们可以通过设定V0=I提供一个不带任何信息的先验,此时我们有p(μ|D,Σ)=N(ˉx,1NΣ),即和MLE得到的结果相同。

Posterior distribution of Σ

现我们讨论如何计算p(Σ|D,μ),其likelihood具有如下形式:6 p(D|μ,Σ)|Σ|N2exp(12tr(SμΣ1))

之前我们提到过如果采用共轭先验能够减少计算的复杂度,而此likelihood的共轭先验就是我们之前提到的非常恐怖的逆Wishart分布,即:

IW(Σ|S10,ν0)|Σ|(ν0+D+1)/2exp(12tr(S0Σ1))

上式中N0=ν+D+1控制着先验的强度,和likelihood中的N的作用基本相同。

将先验和likelihood相乘我们得到如下posterior:

p(Σ|D,μ)|Σ|N2exp(12tr(SμΣ1))|Σ|(ν0+D+1)/2exp(12tr(S0Σ1))=|Σ|N+(ν0+D+1)2exp(12tr(Σ1(Sμ+S0)))=IW(Σ|SN,νN)νN=ν0+NS1N=S0+Sμ

总而言之,从上式我们可以看到,posterior vN的强度为ν0N;posterior离散度矩阵是先验离散度矩阵S0和数据离散度矩阵Sμ之和。

MAP Estimation

根据我们之前得到的关于Σ的极大似然估计值,即:

MLE for Covariance

从上式我们可以看出矩阵的rank为min(N,D)。若N小于D,该矩阵不是full rank的,因此不可逆。另尽管N可能大于D,ˆΣ也可能是ill-conditioned(接近奇异)。

为了解决上述问题,我们可以采用posterior mean或mode.我们可以证明Σ的MAP估计值如下:(使用我们推导MLE时所用的技巧,其实并不难,亲证下式正确):

ˆΣMAP=SNνN+D+1=S0+SμN0+Nμ

当我们采用improper uniform prior,即S0=0,N0=0时,我们即得MLE估计值。

D/N较大时,选择一个包含信息的合适的prior就相当必要了。令μ=ˉx,故有Sμ=Sˉx,此时MAP估计值可被重写为prior mode和MLE的convex combination.令Σ0S0N0为prior mode,则有:

ˆΣMAP=S0+SˉxN0+N=N0N0+NS0N0+NN0+NSN=λΣ0+(1λ)ˆΣmle

其中λ=N0N0+N,控制着向prior shrinkage的程度。对于λ而言,我们可以通过交叉验证设置其值。7

而对于先验的协方差矩阵S0,一般采用如下prior:S0=diag(ˆΣmle).因此,我们有:

S_0

由上式我们可以看出,对角线的元素和极大似然估计值相等,而非对角线则趋近于0.因此该技巧也被称为shrinkage estimation,or regularized estimation.

  1. LDA算法可能导致overfitting,具体解决方法请参阅Machine Learning:a probabilistic perspective一书4.2.*部分
  2. 具体证明请参考ML:APP一书4.3.4.3一节
  3. 证明请参考ML:APP一书4.4.3节
  4. 另外一种方法为Kalman Filter Algorithm
  5. 本部分部分参考Regularized Gaussian Covariance Estimation
  6. 参见MLE for Gaussian部分
  7. 其他方法见ML:APP一书4.6.2.1节
Comment

本文为机器学习系列第二篇,主要研究一下离散数据生成模型。我们会介绍两种模型,一为Dirichlet-multinomial model,一为朴素贝叶斯分类器(Naive Bayes Classfier,NBC),最后我们介绍一下关于互信息的基本知识以及将其用于特征选择。

Read More

苹果的味道 © qingyuanxingsi Powered by Pelican and Twitter Bootstrap. Icons by Font Awesome and Font Awesome More