Welcome to this fancy |

Rayinfos

园龄:2年8个月粉丝:3关注:0

多元高斯分布/高斯过程全解析

大纲

  • 公式推导
  • 参数估计
  • 高斯分布运算
  • 高斯分布性质
  • 高斯过程(Gaussian process)
  • 高斯混合模型

概念区分

  • 边缘分布(marginal distribution)和联合分布
  • 概率密度函数和概率分布函数

1. 多元高斯分布公式推导

首先我们知道一元高斯分布是:N(x|u,σ2)=12πσ2exp[12σ2(xu)2], 拓展到高维时:

N(x|u,Σ)=1(2π)D/21|Σ|1/2exp[12(xu)TΣ1(xu)]

其中,x 表示维度为 D 的向量,u 则是这些向量的平均值,Σ 表示所有向量 x 的协方差矩阵。

现在进行推导。为了简单起见,假设所有变量都是相互独立的,即对于概率分布函数 f(x0,x1,,xn)=f(x0)f(x1)...f(xn) 成立。

假设有很多变量 x=[x1x2],它们的均值为 u=[u1u2],方差为 σ=[σ1σ2]

由于 x1x2 是相互独立的,所以,x 的高斯分布函数可以表示为:

f(x)=f(x1,x2)=f(x1)f(x2)=12πσ12exp(12(x1u1σ1)2)×12πσ22exp(12(x2u2σ2)2)=1(2π)2/2(σ12σ22)1/2exp(12[(x1u1σ1)2+(x2u2σ2)2])

接下来,为了推出文章开篇的高维公式,我们要想办法得到协方差矩阵 Σ
对于二维的向量 x 而言,其协方差矩阵为:

Σ=[σ11σ12σ12σ22]=[σ12σ12σ21σ22]

协方差(Covariance)在概率论和统计学中用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况,即当两个变量是相同的情况。简单来讲,协方差就是衡量两个变量相关性的变量。当协方差为正时,两个变量呈正相关关系(同增同减);当协方差为负时,两个变量呈负相关关系(一增一减)。 而协方差矩阵,只是将所有变量的协方差关系用矩阵的形式表现出来而已。通过矩阵这一工具,可以更方便地进行数学运算。协方差公式为:

Cov(X,Y)=E(X,Y)E(X)E(Y)=i=1n(xix)(yiy)n1$$$X$$Y$$x$$y$$X$$x$$Y$$y$$x$$y$$X$$Y$$Z$$X$$Y$$$Cov(Z)=E[(XE[X])(YE[Y])]=[Cov(X,X)Cov(X,Y)Cov(Y,X)Cov(Y,Y)]

这样矩阵中之中每个元素 Σij=(ii)T(jj)1
X, Y两个变量独立时,Cov(X,Y)为0:

E(XY)=xyx×y×p(x,y)=xyx×y×px(x)×py(y)=xx×px(x)yy×py(y)=E(X)E(Y)

由于 x1x2 是相互独立的,所以 σ12=σ21=0。这样,Σ 退化成 [σ1200σ22]
Σ 的行列式 |Σ|=σ12σ22,代入公式 (4) 就可以得到:

f(x)=1(2π)2/2|Σ|1/2exp(12[(x1u1σ1)2+(x2u2σ2)2])

这样一来,我们已经推出了公式的左半部分,下面,开始处理右面的 exp 函数。
原始的高维高斯函数的 exp 函数为:exp[12(xu)TΣ1(xu)],根据前面算出来的 Σ,我们可以求出它的逆:Σ1=1σ12σ22[σ2200σ12]
接下来根据这个二维的例子,将原始的 exp() 展开:

exp[12(xu)TΣ1(xu)]=exp[12[x1u1   x2u2]1σ12σ22[σ2200σ12][x1u1x2u2]]=exp[12[x1u1   x2u2]1σ12σ22[σ22(x1u1)σ12(x2u2)]]=exp[12σ12σ22[σ22(x1u1)2+σ12(x2u2)2]]=exp[12[(x1u1)2σ12+(x2u2)2σ22]]

展开到最后,发现推出了原公式。说明原公式 N(x|u,Σ)=1(2π)D/21|Σ|1/2exp[12(xu)TΣ1(xu)]是成立的。


2. 参数估计

如果给定了很多数据点,并且知道它们服从某个高斯分布,我们要求高斯分布的参数(μΣ),估计模型参数的方法有很多,最常用的就是极大似然估计(MLE)。对于一维的高斯模型假如有m个数据点,则似然函数:

f(x1,x2,,xm)=i=1m12πσexp((xiμ~)22σ2)

取对数后求导,令导数为 0 得到似然方程。

lnfμ=1σ2i=1m(xiμ~)=0

lnfσ=mσ+1σ3i=1m(xiμ~)=0

得到 μ~=1mi=1mxiσ=1mi=1m(xiμ~)2

多维高斯分布时,假如有 mp 维向量 x,参数估计为:

μ~=1mi=1mx(i)Σ=1mi=1m(x(i)μ~)(x(i)μ~)T

在计算样本协方差矩阵时,我们要使用无偏估计,即将分母由 m 改为 (m1)


3. 高斯分布运算

3.1 一元高斯分布相乘

假设p(x1)=N(x|μ1,σ1),p(x2)=N(x|μ2,σ2)均是关于变量
x的分布,得两个高斯分布相乘仍为缩放的高斯分布:

p(x1)p(x2)=e12σ12(xμ1)2e12σ22(xμ2)2=e12σ12+σ22x22(μ1σ22+μ2σ12)x+constantσ12σ22

则高斯分布的参数:

μ=μ1σ22+μ2σ12σ12+σ22 ,  σ=σ12σ22σ12+σ22

上式可写为如下形式,从而推广至n个一维高斯分布相乘

μ=(μ1σ12+μ2σ22)σ2 ,  1σ2=1σ12+1σ22

新函数等价于正态分布 N(μ,σ2) 的概率密度函数乘以缩放因子 A 。其中,缩放因子A=e(μ1μ2)22(σ12+σ22)2π(σ12+σ22)

3.2 多元高斯分布相乘

μ=Σ(Σ11μ1+Σ21μ2)Σ=(Σ11+Σ21)1

3.3 高斯分布相加

两个高斯分布函数直接相加,很明显不是一个高斯函数。如果两个满足高斯分布的随机变量相加,那么他们的和还是一个高斯分布。具体的,如果 XN(μX,σX2), YN(μY,σY2)Z=X+Y 那么 ZN(μX+μY,σX2+σY2)

需要用到卷积运算:(fg)(n)=f(τ)g(nτ)dτ

FZ(z)=P(Zz)=P(X+Yz)=x+yzf(x,y)dxdy=dxzxf(x,y)dy=令u=y+xdxzf(x,ux)du=zduf(x,ux)dx

所以,Z的概率密度函数为:

fZ(z)=f(x,zx)dx

XY为独立随机变量时,Z的概率密度为fZ(z)=fY(zx)fX(x)dx

法二:使用特征函数证明
高斯分布的特征函数为:

φ(t)=exp(itμσ2t22)

所以,

φX+Y(t)=E(eit(X+Y))=φX(t)φY(t)=exp(itμXσX2t22)exp(itμYσY2t22)=exp(it(μX+μY)(σX2+σY2)t22).

3.4 高斯线性模型

{P(x)=N(x|μ,Λ1)P(y|x)=N(y|Ax+b,L1)P(y)=N(y|Aμ+b,L1+AΛ1AT)P(x|y)=N(x|μ+Λ1AT(L1+AΛ1AT)(yAμb),Λ1Λ1AT(L1+AΛ1AT)1AΛ1)

高斯线性系统推导如下:

p(x)=N(μ0,Σ0)y=Ax+b+ϵ, ϵN(0,Σy)p(y|x)=N(Ax+b,Σy)

yx 产生,在观测到 y 后可以对 x 进行更新(update): p(x|y)=N(μx|y,Σx|y),没观测到 y 可以对其预测(predict)P(y)
下面对μx|y,Σx|y进行计算。已知p(x,y)=p(y|x)p(x)p(y|x)p(x)的指数部分为:

xT(Σ01+ATΣy1A)x2xT(Σ01μ0+ATΣy1(yb))+(yb)TΣy1(yb)+μ0TΣ01μ0

通过配方可以得到:

Σx|y1=Σ01+ATΣy1Aμx|y=Σx|y(Σ01μ0+ATΣy1(yb))

下面对 p(y) 进行求解,我们知道 p(y)=p(y|x)p(x)dxp(x)p(y|x)=p(y)p(x|y) 通过上述的式子,如果对上式求积分或者配方会有些复杂。实际上,通过上式可以得到 p(x,y) 逆协方差矩阵:

Λ=[Σ01+ATΣy1AATΣy1Σy1AΣy1]

利用联合高斯分布的推断结论,可以得到:

μx|y=Σx|y(Σx|y1μ0Λ12(yμy))=Σx|y(Σ01μ0+ATΣy1Aμ0+ATΣy1(yμy))=Σx|y(Σ01μ0+ATΣy1(yb))

可以推知:μy=Aμ0+b, 再对Σy(这里Σy表示p(y)的协方差矩阵)进行计算:

Σ=[ΣxΣxyΣxyTΣy]=Λ1[IAT0I]Λ[I0AI]=[Σ0100Σy1]Λ1=[I0AI][Σ000Σy][IAT0I]=[Σ0Σ0ATAΣ0Σy+AΣ0AT]

因此 Σy=Σy+AΣ0ATp(y) 的分布参数如下:

μ=Aμ0+bΣ=Σy+AΣ0AT


4. 高斯分布性质

多元正态分布有4种等价的定义。

定义1--由标准正态随机向量线性组合得到

U=(U1,U2,,Uq) 为随机向量, U1,,Uq 独立服从标准正态。设 μp 维常数向量, Ap×q 维常数矩阵,则称 X=AU+μ 的分布为 p 元正态分布,或称 Xp 维正态随机向量,记作 XNp(μ,AA)

性质1--特征函数

在概率论中,任何随机变量的特征函数(ch.f)完全定义了它的概率分布。在实直线上,它由以下公式给出,其中X是任何具有该分布的随机变量: φX(t)=E[eitX]

φX(t)=E[eitX]=泰勒展开E(1+itX1t2X22!++(it)nXnn!)=E(1)+E(itX1)E(t2X22!)++E((it)nXnn!)=1+itE[X]一阶矩1t2E[X2]二阶矩2!++(it)nE[Xn]n阶矩n!)

k阶原点矩: E[Xk]  Ak=1ni=1nXik,k=1,2,
k阶中心矩: E[(XE(X))k]  Ak=1ni=1n(XiX¯)k,k=2,3,

可见特征函数包含了分布函数的所有矩(moment),也就是包含了分布函数的所有特征。
所以,特征函数其实是随机变量 X 的分布的另外一种描述方式。
假设某连续随机变量 X的概率密度函数为 f(x),那么可知:E(X)=+xf(x)dx,特征函数为:

φX(t)=E[eitX]=+eitxf(x)dx

特征函数把分布函数换到另外一个坐标系,也可以获得一些计算的好处

  • 假如我们不知道分布函数,但是通过实验算出了期望、方差、偏度、峰度等,那么可以用特征函数去代替分布函数
  • 两个分布函数的卷积 fg 通过特征函数更换坐标系后,可以变为更容易计算的乘法:φ(fg)=φ(f)φ(g)
  • 通过对 t 求导,可以简单求出各阶矩:φX(k)(0)=ikE[Xk]

由定义1得到的随机向量 X 的特征函数为

ΦX(t)=exp[itμ12tAAt]

其中 t=(t1,,tp)Rp
证明:首先考虑一维标准正态分布的特征函数为 ΦUi(ti)=exp[12ti2]
根据独立性有

ΦU(t)=exp[12i=1qti2]=exp[12tt]

进而根据 X 的定义得到

ΦX(t)=E[exp{itX}]=E[exp{it(AU+μ)}]=E[exp{itμ}]E[exp{itAU}]=E[exp{itμ}]E[exp{i(At)U}]

其中 E[exp{i(At)U}]ΦU(At) ,代入即得结论.

定义2--由特征函数定义

如果随机向量 X 的特征函数具有如下形式 ΦX(t)=exp[itμ12tΣt], 则称 X 服从 p 维正态分布,记作 XNp(μ,Σ)

性质2--正态随机向量任意线性变换仍服从正态分布

XNp(μ,Σ),令 Z=BX+d ,则 ZNs(Bμ+d,BΣB) ,其中 Bs×q 维矩阵,ds 维向量.

推论--子向量的均值与协方差:

X=[X(1)X(2)]rprNp(μ,Σ) ,将 μ, Σ 分为

μ=[μ(1)μ(2)]rpr,Σ=[Σ11Σ12Σ21Σ22]rpr

则有 X(1)Nr(μ(1),Σ11),X(2)Npr(μ(2),Σ22)
注意: Σ12Σ21 ,两者互为转置

性质3--多元正态 任意线性组合为一元正态

X=(X1,X2,,Xp)p 维随机向量,则 X 服从 p 元正态分布等价于对任意 p 维实向量, ξ=aX 是一维正态随机变量.
证明
当 X 为 p 元正态分布,由性质2知 ξ 为一维正态随机变量。
反之,如果对任意 aξ=aX 为一维正态随机变量,则 ξ 各阶矩存在,进而 X 的均值和协方差存在,分别设为 μ,Σ ,则

ξ=aXN(aμ,aΣa)

进而考察 X 的特征函数得到

ΦX(a)=exp[iaX]=exp[iξ]=Φξ(1)=exp[iaμ12aΣa]

刚好等于多元正态的特征函数,由特征函数与分布的一一对应得到结论.

定义3--任意线性组合为正态

如果 p 维随机向量 X 的任意线性组合均服从一元正态分布,则称 Xp正态随机向量.

性质4--联合密度函数

如果 XNp(μ,Σ)Σ>0 ,则 X联合密度函数

f(x)=1(2π)p/2|Σ|1/2exp[12(xμ)Σ1(xμ)]

定义4--密度函数

如果 p 维随机向量 X 的联合密度函数为

f(x)=1(2π)p/2|Σ|1/2exp[12(xμ)Σ1(xμ)]

则称 Xp 维正态随机向量.
注意:定义4要求 Σ>0 ,其他三个只要求 Σ0 ,一般也不考虑 X 为退化随机向量的情况.

5. 高斯条件分布和独立性

仅讨论 Σ0 (即半正定) 的情形

定理1--正态随机向量的独立性等价于协方差为0矩阵

定理2--条件分布

X=[X(1)X(2)]rprNp(μ,Σ)(Σ>0) ,则当 X(2)=x(2) 给定时, X(1) 的条件分布为

(X(1)X(2)=x(2))Nr(μ12,Σ112)

其中

μ12=μ(1)+Σ12Σ221(x(2)μ(2))Σ12=Σ11Σ12Σ221Σ21

证明:从回归的角度会比较容易理解,理论依据是,在均方意义下,线性回归的结果就是条件期望。将 X 中心化后做回归

X(1)μ(1)=β(X(2)μ(2))+ε

那么 β(x(2)μ(2)) 就是 X(1)μ(1) 的条件期望。现在假设对于每个变量,都有 n 个观测数据,并将回归形式改写为 Zt=βRt+ε ,那么利用最小二乘估计可以得到参数的估计量为 β=(RR)1RZ 。考虑当 n 充分大的情况下, (RR)1 对应了 Σ221RZ 对应了 Σ21 进而对 β 求转置后得到

X(1)=μ(1)+Σ12Σ221(X(2)μ(2))+ε

因此条件期望就是

μ12=μ(1)+Σ12Σ221(x(2)μ(2))

下面考虑条件方差的计算。做回归后得到的误差项 ε 是从 X(1) 中剔除了 X(2) 对其的影响,那么条件方差就应该等于误差项的方差,即

Σ12=Varε=Var(X(1)μ(1))Var[Σ12Σ221(X(2)μ(2))]=Σ11Σ12Σ221Σ22(Σ12Σ221)=Σ11Σ12Σ221Σ21

由此可以自然地得到下面的推论:
X(2)X(1)Σ12Σ221X(2) 相互独立
X(1)X(2)Σ21Σ111X(1) 相互独立
X(2)X(1)Npr(μ21,Σ21), 其中

μ21=μ(2)+Σ21Σ111(x(1)μ(1))Σ21=Σ22Σ21Σ111Σ12

问:如果是三个子向量,给定其中两个,求另一个的条件分布呢?
答:把给定的两个看做一个子向量就可以。

条件数字特征

就是刚刚推导的东西的定义

(1)条件期望(Conditional Expectation),回归系数(regression coefficient),偏相关系数(Partial correlation coefficient)
X=[X(1) X(2)]Np([μ(1) μ(2)],[Σ11Σ12 Σ21Σ22])

根据定理2有 (X(1)X(2)=x(2))Nr(μ12,Σ12),我们把

μ12=μ(1)+Σ12Σ221(x(2)μ(2))

称为条件期望(Conditional Expectation),记作 E(X(1)X(2)=x(2)) ;把 Σ12Σ221=defB 称为回归系数.

区分 E(X)E(XY)E(XY=y)
E(X):一个数
E(XY) :随机变量,关于 Y 的函数,没有固定的 y 值
E(XY=y) : y 的函数 f(y) ,对于给定的 y ,有唯一确定值与之对应

全期望公式(Law of total expectation)
X,Y 为离散型随机变量,下列期望和条件期望均存在,则

E(X)=E(E(XY))

Y 为连续型随机变量,则

E(X)=E(E(XY))=+E(XY=y) dFY(y)

若 Y 为离散型随机变量,则

E(X)=E(E(XY))=yE(XY=y)P(Y=y)

离散型的证明如下:

E(E(XY))=yE(XY=y)P(Y=y)=y(xxP(X=xY=y))P(Y=y)=yxxP(X=xY=y)P(Y=y)=yxxP(Y=yX=x)P(X=x)=xyxP(Y=yX=x)P(X=x)=xxP(X=x)(yP(Y=yX=x))=xxP(X=x)=E(X)

一个特殊情况:若 {Ai}i 是一个样本空间的有限集或可列集,则

E(X)=E(E(XY))=iE(XAi)P(Ai)

为了定义偏回归系数,将条件方差矩阵的元素具体表示为

Σ12=(σij)r×r(i,j=1,,r)

ρijr+1,,p=σijσiiσjj偏相关系数,即为 X(2)=(Xr+1,,Xp) 给定的条件下, Xi,Xj相关系数.

(2)全相关系数(了解)
Z=[XY]p1Np+1([μXμy],[ΣXXΣXyΣyXσyy]),则称

R=(yXΣXX1ΣXyσyy)1/2

YX=(X1,,Xp) 的全相关系数.

(3)最佳预测
X(1)=Y,g(x(2))=E(YX(2)=x(2)) ,则对任意函数 ϕ() ,可以证明

E[(Yg(x(2)))2]E[(Yϕ(x(2)))2]

也就是在均方准则下,条件期望是最优预测,证明方法就是加一项减一项,往证交叉项为0.

5.2 高斯边缘分布

如果联合分布 p(xa,xb) 是高斯分布,那么条件概率分布 p(xa|xb) 也是高斯分布,那么边缘概率分布 p(xa)=p(xa,xb) dxb 显然也是一个高斯分布。
我们主要研究联合分布的指数项二次型,这次考虑涉及到 xb 的项,结合条件高斯分布中对 z2=zTz=(xμx)TΣ1(xμx) 几何形式关于 xa,xb 的分解公式,可以得到:

12xbTΛbbxb+xbTm=12(xbΛbb1m)TΛbb(xbΛbb1m)+12mTΛbb1m

其中 m=ΛbbμbΛba(xaμa)ΛΣ1Λ为协方差矩阵的逆矩阵,又称为精度矩阵
上式中与 xb 相关的项转化为一个高斯分布的标准二次型,结合边缘概率公式需要积分:

exp{12(xbΛbb1m)TΛbb(xbΛbb1m)}dxb

上面只提出了关于 xb 的二项式,其最后一项 12mTΛbb1m 为和 xb 无关但和 xa 有关的项,结合前文提到的,除 xb 二次项以外的并和 xa 有关的项结合,得到:

(1)12[ΛbbμbΛba(xaμa)]TΛbb1[ΛbbμbΛba(xaμa)](2)12xaTΛaaxa+xaT(ΛaaΛabΛbb1Λba)xa(3)=12xaT(ΛaaΛabΛbb1Λbaxa)+xaT(ΛaaΛabΛbb1Λba)μa+b

b 为常数,是与 xa 无关的量,那么可以得到边缘概率的协方差矩阵:
Σa=(ΛaaΛabΛbb1Λba)1
均值为: μa=Σa(ΛaaΛabΛbb1Λba)μa
前文介绍过分块矩阵逆矩阵的恒等式,那么可以得出:
Σaa=(ΛaaΛabΛbb1Λba)1
最后可以得出边缘概率 p(xa) 的均值和协方差:
E[xa]=μacov[xa]=Σaa
边缘概率分布:
p(xa)=N(xa|μa,Σaa)

5.3 混合高斯分布

通过将更基本的概率分布(高斯分布)进行线性组合叠加,然后形式化为概率模型,被称为混合模型。高斯分布的线性组合可以给出相当复杂的概率密度形式。通过使用足够多的高斯分布,并且调节它们的均值和方差以及线性组合的系数、几乎所有的连续概率密度能够以任意的精度近似。考虑 K 个高斯概率密度的叠加,形式为:
p(x)=k=1Kπk N(x|μk,Σk)
称为混合高斯分布,每个高斯概率密度 N(x|μk,Σk) 被称为混个高斯分布的一个成分,并且有自己的均值和协方差 μkΣkπk 被称为混合系数,可以得到: k=1Kπk=1
根据概率的加和规则和乘积规则,边缘概率密度为: p(x)=k=1Kp(k)p(x|k)
这和上面的混合高斯分布公式是等价的,把 πk=p(k) 看成第k个成分的先验概率,把密度 N(x|μk,Σk)=p(x|k) 看成以k为条件的x的概率。
后验概率 p(k|x)=p(k)p(x|k)Σlp(l)p(x|l)=πkN(x|μk,Σk)ΣlπlN(x|μl,Σl)

π{π1,...,πK},μ{μ1,...,μK},Σ{Σ1,...,ΣK} ,对数似然函数为:

ln p(X|π,μ,Σ)=n=1Nln{xk=1Kπk N(xn|μk,Σk)}X={x1,...,xN}

因为该对数似然函数中对数里含有求和式,不能像一元高斯分布那样可以求得封闭的解析解,可以通过迭代数值优化方法以及期望最大化方法来求解。

6. 高斯过程

6.1 简介

高斯过程(Gaussian process, GP) 是一个概率统计学上的概念,更确切的说应该是随机过程(Stochastic process)中一个特殊例子。
在高斯过程中,连续输入空间中每个点都是与一个正态分布的随机变量相关联。此外,这些随机变量的每个有限集合都有一个多元正态分布。高斯过程的分布是所有那些(无限多个)随机变量的联合分布,正因如此,它是连续域(例如时间或空间)的分布。

【定义】 对于一个连续域 T (假设他是一个时间轴),如果我们在连续域上任选 n 个时刻: t1,t2,t3,...,tnT ,使得获得的一个 n 维向量 {ξ1,ξ2,ξ3,...,ξn} 都满足其是一个 n 维高斯分布,那么这个 {ξt} 就是一个高斯过程。

GP可以被mean和covariance function共同唯一决定其表达式,具体的,ξtGP(m(t),k(t,s))。因为我们知道一个高斯分布可以被mean和variance共同唯一决定,一个多元高斯分布可以对mean vector和covariance matrix共同唯一决定。mean需要描述每一个时间点 t 上的均值,但是这个时候就不能用向量了,因为是在连续域上的,维数是无限的,因此就应该定义成一个关于时刻 t 的函数: m(t) 。covariance function被称为核函数kernel,原因就是它捕捉了不同输入点之间的关系,并且反映在了之后样本的位置上。这样的话,就可以做到,利用点与点之间关系,以从输入的训练数据预测未知点的值。比如径向基函数 RBF:

k(s,t)=σ2exp(||st||22l2)

这里 σl 是径向基函数的超参数。st 表示高斯过程连续域上的两个不同的时间点, ||st||2 是一个二范式,简单点说就是 st 之间的距离,径向基函数输出的是一个标量,他代表的就是 st 两个时间点各自所代表的高斯分布之间的协方差值,很明显径向基函数是一个关于 st 距离负相关的函数,两个点距离越大,两个分布之间的协方差值越小,即相关性越小,反之越靠近的两个时间点对应的分布其协方差值就越大。

import numpy as np
def gaussian_kernel(x1, x2, l=1.0, sigma_f=1.0):
m, n = x1.shape[0], x2.shape[0]
dist_matrix = np.zeros((m, n), dtype=float)
for i in range(m):
for j in range(n):
dist_matrix[i][j] = np.sum((x1[i] - x2[j]) ** 2)
return sigma_f ** 2 * np.exp(- 0.5 / l ** 2 * dist_matrix)
train_X = np.array([1, 3, 7, 9]).reshape(-1, 1)#转换为4*1矩阵形式
print(gaussian_kernel(train_X, train_X)) #4*4矩阵,当 i=j 时,就是自身的方差

6.2 高斯过程回归

我们知道,高斯分布有一个很好的特性,那就是高斯分布的联合概率、边缘概率、条件概率仍然是满足高斯分布的,假设:
n 维的随机变量满足高斯分布: xN(μ,Σn×n)
那么如果我们把这个 n维的随机变量分成两部分: p 维的 xaq 维的 xb ,满足 n=q+p ,那么按照均值向量 μ 和协方差矩阵 Σ 的分块规则,就可以写成:

x=[xaxb]p+qμ=[μaμb]p+qΣ=[ΣaaΣabΣbaΣbb]

那么依据高斯分布的性质,我们知道下列条件分布依然是一个高维的高斯分布:

xb|xaN(μb|a,Σb|a)μb|a=ΣbaΣaa1(xaμa)+μbΣb|a=ΣbbΣbaΣaa1Σab

也就是说,设置了高斯过程的先验参数,一旦我们拿到一些观测值,那么就可以对高斯过程的均值函数和核函数进行修正,得到一个修正后的后验高斯过程,而更新后验参数的信息就来自于观测值。

构建先验分布

GP超参数的估计

将高斯过程对比高维高斯分布,我们把均值向量替换成均值函数,把协方差矩阵替换成核函数,就能够得到高斯过程基于观测值的后验过程的参数表达式:
我们的一组观测值,他们的时刻对应一个向量 X ,那么对应的值时另一个同纬度的向量的 Y ,假设有4组观测值,那就是 {(X[1],Y[1]),((X[2],Y[2])),((X[3],Y[3])),((X[4],Y[4]))}
那么余下的所有非观测点,在连续域上我们定义为 X ,值定义为 f(X)

首先,联合分布显然是满足无限维高斯分布的:

[Yf(X)]N([μ(X)μ(X)][k(X,X)k(X,X)k(X,X)k(X,X)])

从这个联合分布所派生出来的条件概率 f(X)|Y 同样也服从高斯分布,当然这里指的是无限维高斯分布,其实对比一下,把 Y 看作是 xa ,把 f(X) 看作是 xb ,直接类比条件分布的参数表达式:f(X)|YN(μ,k) 这里面的 μk 就是条件分布下的后验高斯过程的均值函数和核函数的形式。

类比我们就可以写成表达式:

μ=k(X,X)k(X,X)1(Yμ(X))+μ(X)k=k(X,X)k(X,X)k(X,X)1k(X,X)

以下例子中,高斯过程先验我们设置均值函数为 μ(X)=0 ,径向基函数 k(X,X)=σ2exp(||XX||22l2)X=[1,3,7,9] 的位置上设置一组观测值, YX取正弦的基础上加上一点高斯噪声。我们在四个观测点的基础上,来求高斯过程的后验。

import matplotlib.pyplot as plt
import numpy as np
#高斯核函数
def gaussian_kernel(x1, x2, l=0.5, sigma_f=0.2):
m, n = x1.shape[0], x2.shape[0]
dist_matrix = np.zeros((m, n), dtype=float)
for i in range(m):
for j in range(n):
dist_matrix[i][j] = np.sum((x1[i] - x2[j]) ** 2)
return sigma_f ** 2 * np.exp(- 0.5 / l ** 2 * dist_matrix)
#生成观测值,取sin函数没有别的用意,单纯就是为了计算出Y
def getY(X):
X = np.asarray(X)
Y = np.sin(X)*0.4 + np.random.normal(0, 0.05, size=X.shape)
return Y.tolist()
#根据观察点X,修正生成高斯过程新的均值和协方差
def update(X, X_star):
X = np.asarray(X)
X_star = np.asarray(X_star)
K_YY = gaussian_kernel(X, X) # K(X,X)
K_ff = gaussian_kernel(X_star, X_star) # K(X*, X*)
K_Yf = gaussian_kernel(X, X_star) # K(X, X*)
K_fY = K_Yf.T # K(X*, X) 协方差矩阵是对称的,因此分块互为转置
K_YY_inv = np.linalg.inv(K_YY + 1e-8 * np.eye(len(X))) # (N, N)
mu_star = K_fY.dot(K_YY_inv).dot(Y)
cov_star = K_ff - K_fY.dot(K_YY_inv).dot(K_Yf)
return mu_star, cov_star
f, ax = plt.subplots(2, 1, sharex=True,sharey=True)
#绘制高斯过程的先验
X_pre = np.arange(0, 10, 0.1)
mu_pre = np.array([0]*len(X_pre))
Y_pre = mu_pre
cov_pre = gaussian_kernel(X_pre, X_pre)
uncertainty = 1.96 * np.sqrt(np.diag(cov_pre))#取95%置信区间
ax[0].fill_between(X_pre, Y_pre + uncertainty,Y_pre - uncertainty, alpha=0.1)
ax[0].plot(X_pre, Y_pre, label="expection")
ax[0].legend()
#绘制基于观测值的高斯过程后验
X = np.array([1, 3, 7, 9]).reshape(-1, 1)#4*1矩阵
Y = getY(X)
X_star = np.arange(0, 10, 0.1).reshape(-1, 1)
mu_star, cov_star = update(X, X_star)
Y_star = mu_star.ravel()
uncertainty = 1.96 * np.sqrt(np.diag(cov_star))#取95%置信区间
ax[1].fill_between(X_star.ravel(), Y_star + uncertainty, Y_star - uncertainty, alpha=0.1)
ax[1].plot(X_star, Y_star, label="expection")
ax[1].scatter(X, Y, label="observation point", c="red", marker="x")
ax[1].legend()
plt.show()

在这里插入图片描述

  • 优点

    • (采用 RBF 作为协方差函数)具有平滑性质,能够拟合非线性数据
    • 高斯过程回归天然支持得到模型关于预测的不确定性(置信区间),直接输出关于预测点值的概率分布
    • 通过最大化边缘似然这一简洁的方式,高斯过程回归可以在不需要交叉验证的情况下给出比较好的正则化效果
  • 缺点

    • 高斯过程是一个非参数模型,每次的 inference 都需要对所有的数据点进行(矩阵求逆)。对于没有经过任何优化的高斯过程回归,n 个样本点时间复杂度大概是 O(n3) ,空间复杂度是 O(n2),在数据量大的时候高斯过程变得 intractable
    • 高斯过程回归中,先验是一个高斯过程,likelihood 也是高斯的,因此得到的后验仍是高斯过程。在 likelihood 不服从高斯分布的问题中(如分类),需要对得到的后验进行 approximate 使其仍为高斯过程
    • RBF 是最常用的协方差函数,但在实际中通常需要根据问题和数据的性质选择恰当的协方差函数

References

多元高斯分布完全解析 -知乎
高斯分布相乘、积分整理
多维高斯分布 -博客园
多维正态分布的最大似然估计 -博客园
高斯性质:多元统计分析第01讲--多元正态分布及参数估计(随机向量,多元正态分布定义,条件分布和独立性)-知乎
高斯条件分布:多元统计分析第02讲(条件分布,随机阵的正态分布,参数估计)
高斯边缘分布:深入理解高斯分布

Sum of normally distributed random variables
第三章·随机向量 ----概率论与数理统计
条件期望与全期望公式

汇总型
prml -gitbook
Gaussian Processes for Machine Learning.pdf

本文作者:Rayinfos

本文链接:https://www.cnblogs.com/rayinfos/p/17814095.html

版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。

posted @   Rayinfos  阅读(3579)  评论(1编辑  收藏  举报
点击右上角即可分享
微信分享提示
评论
收藏
关注
推荐
深色
回顶
收起
  1. 1 Singing Rib
  2. 2 Focus Rib
  3. 3 fossil Rib
  4. 4 damn 藤井風
  5. 5 きらり 藤井風
  6. 6 帰ろう 藤井風
Focus - Rib
00:00 / 00:00
An audio error has occurred, player will skip forward in 2 seconds.