在上一篇中我们着重介绍了对于离散数据的生成模型,紧接上一篇,本篇我们介绍对于连续数据的生成模型。好吧,废话我们就不多说了,直接进入正文。
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=U−TΛ−1U−1=UΛ−1UT=D∑i=11λiuiuTi
其中,ui是U的第i列。因此Mahalanobis距离可被改写为:
(x−μ)TΣ−1(x−μ)=(x−μ)T(D∑i=11λiuiuTi)(x−μ)=D∑i=11λi(x−μ)TuiuTi(x−μ)=D∑i=1y2ixi
其中,yi≜uTi(x−μ)。另2维空间中的椭圆方程为:
y21λ1+y22λ2=1
因此我们可知Gaussian概率密度的等高线沿着椭圆分布,特征向量决定椭圆的朝向,而特征值则决定椭圆有多“椭”
。一般来说,如果我们将坐标系移动μ,然后按U旋转,此时的欧拉距离即为Mahalanobis距离。
MLE for MVN
以下我们给出MVN参数的MLE(极大似然估计)的证明:
Theorem 3.1 若我们获取的N个独立同分布的样本xi∼ N(x|μ,Σ),则关于μ以及 Σ的极大似然分布如下:
我们不加证明地给出如下公式组:
∂(bTa)∂a=b∂(aTAa)∂a=(A+AT)a∂∂Atr(BA)=BT∂∂Alog|A|=A−Ttr(ABC)=tr(CAB)=tr(BCA)
最后一个等式称为迹的循环置换性质(cyclic permutation property)。利用这个性质,我们使用trace trick
可以得到下式:
xTAx=tr(xTAx)=tr(xxTA)=tr(AxxT)
证明:
对数似然函数为:
l(μ,Σ)=logp(D|μ,Σ)=N2log|Λ|−12N∑i=1(xi−μ)TΛ(xi−μ)
其中,Λ=Σ−1为精度矩阵。令yi=xi−μ并利用链式法则有:
∂∂μ(xi−μ)TΣ−1(xi−μ)=∂∂yiyTiΣ−1yi∂yi∂μ=−(Σ−T+Σ−1)yi
即:
∂∂μl(μ,Σ)=−12N∑i=1−2Σ−1(xi−μ)=0
故有:
ˆμ=1NN∑i=1xi=ˉx
即最大似然均值即为经验均值。
利用trace trick我们重写对数似然函数为:
l(Λ)=N2log|Λ|−12∑itr[(xi−μ)(xi−μ)TΛ]=N2log|Λ|−12tr[SuΛ]
其中,Su≜∑Ni=1(xi−μ)(xi−μ)T,业界尊称其为分散度矩阵(Scatter Matrix
),以后我们聊LDA的时候会再次碰到。对Λ求偏导有:
∂l(Λ)∂Λ=N2Λ−T−12STu=0
得:
ˆΣ=1NN∑i=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类情形下可能的决策边界形状,如下图所示:
Linear Discriminant Analysis(LDA)
当各个Gaussian Component的协方差矩阵相同时,此时我们有:
p(y=c|x,θ)∝πcexp[μTcΣ−1x−12xTΣ−1x−12μTcΣ−1μc]=exp[μTcΣ−1x−12μ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+γc∑c′eβTc′x+γc′=S(η)c
其中η=[βT1x+γ1,...,βTCx+γc],S为softmax函数(类似于max函数,故得此名),定义如下:
若将ηc除以一个常数(temperature),当T→0时,我们有:
换句话说,当温度很低时,分布集中在概率最大的那个状态上,而当温度高时,所有的状态呈现均匀分布。
NOTE:该术语来自于统计物理学,在统计物理学中,人们更倾向于使用波尔兹曼分布(Boltzmann distribution)一词。
关于式16
有一个有趣的性质,即对该式取log,我们则会得到一个关于x的线性方程。因此对于任意两类之间的决策边界将会是一条直线,据此该模型也被称为线性判别分析(Linear Discriminant Analysis,LDA)。而且对于二分类问题,我们可以得到:
p(y=c|x,θ)=p(y=c′|x,θ)
βTcx+γc=βTc′x+γ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(x−x0))
因此最后的决策规则很简单:将x平移x0,然后投影到ω上,通过结果是正还是负决定它到底属于哪一类。
当Σ=σ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
给定两变量,x和y.令x∈RDx为一隐含变量,y∈RDy为关于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(y−b)+Σ−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) 为了更直观地解释以上模型,如下图所示:
我们可将x视为2维空间中一个物体的真实位置(但我们并不知道),例如一枚导弹或者一架飞机,yi则是我们的观测值(含噪声),可以视为雷达上的一些点。当我们得到越来越多的点时,我们就能够更好地进行定位。4
现假定我们有多个测量设备,且我们想利用多个设备的观测值进行估计,这种方法称为sensor fusion
.如果我们有具有不同方差的多组观测值,那么posterior将会是它们的加权平均。如下图所示:
我们采用不带任何信息的关于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|Λ|(ν−D−1)/2exp(−12tr(ΛS−1))
其中,ν为自由度,S为缩放矩阵。其归一项具有如下形式:(很恐怖,对吧!
)
ZWi=2νD/2ΓD(ν/2)|S|ν/2
其中,ΓD(a)为多变量Gamma函数: ΓD(x)=πD(D−1)/4D∏i=1Γ(x+(1−i)/2)
于是有Γ1(a)=Γ(a)且有: ΓD(ν0/2)=D∏i=1Γ(ν0+1−i2)
仅当ν>D−1时归一项存在。
其实Wishart分布和Gaussian分布是有联系的。具体而言,令xi∼ N(0,Σ),则离散度矩阵S=∑Ni=1xixTi服从Wishart分布,且S∼ Wi(Σ,1)。于是有E(S)=NΣ.
更一般地,我们可以证明Wi(S,ν)的mean和mode具有如下形式: mean=νS,mode=(ν−D−1)S
仅当ν>D+1时mode存在。
当D=1时,Wishart分布退化为Gamma分布,且有: Wi(λ|s−1,v)=Ga(λ|ν2,s2)
Inverse Wishart Distribution
若Σ−1∼ Wi(S,ν),则Σ∼ IW(S−1,ν+D+1),其中IW为逆Wishart分布。当ν>D−1且S≻0时,我们有: IW(Σ|S,ν)=1ZIW|Σ|−(ν+D+1)/2exp(−12tr(S−1Σ−1))ZIW=|S|−ν/22νD/2ΓD(ν/2)
此外我们可以证明逆Wishart分布具有如下性质: mean=S−1ν−D−1,mode=S−1ν+D+1
当D=1时,它们退化为逆Gamma分布: IW(σ2|S−1,ν)=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)V−1N=V−10+NΣ−1mN=VN(Σ−1(Nˉx)+V−10m0)
我们可以通过设定V0=∞I提供一个不带任何信息的先验,此时我们有p(μ|D,Σ)=N(ˉx,1NΣ),即和MLE得到的结果相同。
Posterior distribution of Σ
现我们讨论如何计算p(Σ|D,μ),其likelihood具有如下形式:6 p(D|μ,Σ)∝|Σ|−N2exp(−12tr(SμΣ−1))
之前我们提到过如果采用共轭先验能够减少计算的复杂度,而此likelihood的共轭先验就是我们之前提到的非常恐怖的逆Wishart分布,即:
IW(Σ|S−10,ν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+NS−1N=S0+Sμ
总而言之,从上式我们可以看到,posterior vN的强度为ν0加N;posterior离散度矩阵是先验离散度矩阵S0和数据离散度矩阵Sμ之和。
MAP Estimation
根据我们之前得到的关于Σ的极大似然估计值,即:
从上式我们可以看出矩阵的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.令Σ0≜S0N0为prior mode,则有:
ˆΣMAP=S0+SˉxN0+N=N0N0+NS0N0+NN0+NSN=λΣ0+(1−λ)ˆΣmle
其中λ=N0N0+N,控制着向prior shrinkage
的程度。对于λ而言,我们可以通过交叉验证设置其值。7
而对于先验的协方差矩阵S0,一般采用如下prior:S0=diag(ˆΣmle).因此,我们有:
由上式我们可以看出,对角线的元素和极大似然估计值相等,而非对角线则趋近于0.因此该技巧也被称为shrinkage estimation,or regularized estimation.
- LDA算法可能导致overfitting,具体解决方法请参阅Machine Learning:a probabilistic perspective一书4.2.*部分 ↩
- 具体证明请参考ML:APP一书4.3.4.3一节 ↩
- 证明请参考ML:APP一书4.4.3节 ↩
- 另外一种方法为Kalman Filter Algorithm ↩
- 本部分部分参考Regularized Gaussian Covariance Estimation ↩
- 参见MLE for Gaussian部分 ↩
- 其他方法见ML:APP一书4.6.2.1节 ↩