当前位置:网站首页>高维高斯分布基础
高维高斯分布基础
2022-08-01 01:08:00 【Adenialzz】
高维高斯分布基础
多位高斯分布的几何理解
多维高斯分布表达式为:
p ( x ∣ μ , Σ ) = 1 ( 2 π ) p / 2 ∣ Σ ∣ 1 / 2 e − 1 2 ( x − μ ) T Σ − 1 ( x − μ ) p(x|\mu,\Sigma)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}e^{-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)} p(x∣μ,Σ)=(2π)p/2∣Σ∣1/21e−21(x−μ)TΣ−1(x−μ)
其中 x , μ ∈ R p , Σ ∈ R p × p x,\mu\in\mathbb{R}^{p},\Sigma\in\mathbb{R}^{p\times p} x,μ∈Rp,Σ∈Rp×p , Σ \Sigma Σ 为协方差矩阵,一般而言是半正定矩阵,这里我们强化一下条件,只考虑正定矩阵。
首先我们处理指数上的数字,指数上的数字可以记为 x x x 和 μ \mu μ 之间的马氏距离。
马氏距离(Mahalanobis Distance)是度量学习中一种常用的距离指标,同欧氏距离、曼哈顿距离、汉明距离等一样被用作评定数据之间的相似度指标。但却可以应对高维线性分布的数据中各维度间非独立同分布的问题。
两个向量 x \bf{x} x 和 y \mathbf{y}{} y 之间的马氏距离为:
D M ( x , y ) = ( x − y ) T Σ − 1 ( x − y ) ) D_M(\bf{x},\bf{y})=\sqrt{(\bf{x}-\bf{y})^T\Sigma^{-1}(\bf{x}-\bf{y}))} DM(x,y)=(x−y)TΣ−1(x−y))
其中 Σ \Sigma Σ 是多维随机变量的协方差矩阵, μ \mu μ 为样本均值,如果协方差矩阵是单位向量( Σ = I \Sigma=I Σ=I),也就是各维度独立同分布,马氏距离就变成了欧氏距离。关于马氏距离,详见:马氏距离(Mahalanobis Distance)。
对于对称的协方差矩阵可进行特征值分解, Σ = U Λ U T = ( u 1 , u 2 , ⋯ , u p ) d i a g ( λ i ) ( u 1 , u 2 , ⋯ , u p ) T = ∑ i = 1 p u i λ i u i T \Sigma=U\Lambda U^{T}=(u_{1},u_{2},\cdots,u_{p})diag(\lambda_{i})(u_{1},u_{2},\cdots,u_{p})^{T}=\sum\limits _{i=1}^{p}u_{i}\lambda_{i}u_{i}^{T} Σ=UΛUT=(u1,u2,⋯,up)diag(λi)(u1,u2,⋯,up)T=i=1∑puiλiuiT ,于是:
Σ − 1 = ∑ i = 1 p u i 1 λ i u i T \Sigma^{-1}=\sum\limits _{i=1}^{p}u_{i}\frac{1}{\lambda_{i}}u_{i}^{T} Σ−1=i=1∑puiλi1uiT
Δ = ( x − μ ) T Σ − 1 ( x − μ ) = ∑ i = 1 p ( x − μ ) T u i 1 λ i u i T ( x − μ ) = ∑ i = 1 p y i 2 λ i \Delta=(x-\mu)^{T}\Sigma^{-1}(x-\mu)=\sum\limits _{i=1}^{p}(x-\mu)^{T}u_{i}\frac{1}{\lambda_{i}}u_{i}^{T}(x-\mu)=\sum\limits _{i=1}^{p}\frac{y_{i}^{2}}{\lambda_{i}} Δ=(x−μ)TΣ−1(x−μ)=i=1∑p(x−μ)Tuiλi1uiT(x−μ)=i=1∑pλiyi2
我们注意到 y i y_{i} yi 是 x − μ x-\mu x−μ 在特征向量 u i u_{i} ui 上的投影长度,因此上式子就是 Δ \Delta Δ 取不同值时的同心椭圆。例如,在维度 P = 2 P=2 P=2 时,取 Δ = 1 \Delta=1 Δ=1,则有: y 1 2 λ 1 + y 2 2 λ 2 = 1 \frac{y_1^2}{\lambda_1}+\frac{y_2^2}{\lambda_2}=1 λ1y12+λ2y22=1 ,明显就是椭圆的曲线方程。
多维高斯模型的限制
下面我们看多维高斯模型在实际应用时的两个限制:
- 参数 Σ , μ \Sigma,\mu Σ,μ 的自由度为 O ( p 2 ) O(p^{2}) O(p2) 对于维度很高的数据其自由度太高。
- 解决方案:高自由度的来源是 Σ \Sigma Σ 有 p ( p + 1 ) 2 \frac{p(p+1)}{2} 2p(p+1) 个自由参数,可以假设其是对角矩阵,甚至假设其对角线上的元素都相同,此时称为各向同性的高斯分布。前一种的算法有 Factor Analysis,后一种有概率 PCA (p-PCA) 。
- 第二个问题是单个高斯分布是单峰的,对有多个峰的数据分布不能得到好的结果。
- 解决方案:使用多个单高斯模型组合得到高斯混合模型 GMM。
多维高斯分布的边缘概率和条件概率
对于高斯模型的线性变换,有这样一个定理(暂不证明):
定理:已知 x ∼ N ( μ , Σ ) , y ∼ A x + b x\sim\mathcal{N}(\mu,\Sigma), y\sim Ax+b x∼N(μ,Σ),y∼Ax+b, x ∈ R p , y ∈ R p x\in\mathbb{R}^p,y\in\mathbb{R}^p x∈Rp,y∈Rp,那么 y ∼ N ( A μ + b , A Σ A T ) , Σ ∈ R p × p , A ∈ R 1 × p y\sim\mathcal{N}(A\mu+b, A\Sigma A^T),\ \ \Sigma \in \mathbb{R}^{p\times p },\ \ A\in\mathbb{R}^{1\times p} y∼N(Aμ+b,AΣAT), Σ∈Rp×p, A∈R1×p。
我们将 p p p 维样本数据拆分为 m + n m+n m+n 维: x = ( x 1 , x 2 , ⋯ , x p ) T = ( x a , m × 1 , x b , n × 1 ) T x=(x_1, x_2,\cdots,x_p)^T=(x_{a,m\times 1}, x_{b,n\times1})^T x=(x1,x2,⋯,xp)T=(xa,m×1,xb,n×1)T 。
对应的高斯模型的参数也进行拆分:均值 μ = ( μ a , m × 1 , μ b , n × 1 ) T \mu=(\mu_{a,m\times1}, \mu_{b,n\times1})^T μ=(μa,m×1,μb,n×1)T,协方差矩阵 Σ = ( Σ a a Σ a b Σ b a Σ b b ) \Sigma=\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix} Σ=(ΣaaΣbaΣabΣbb),已知 x ∼ N ( μ , Σ ) x\sim\mathcal{N}(\mu,\Sigma) x∼N(μ,Σ)。
下面介绍如何求多维高斯分布的边缘概率和条件概率 p ( x a ) , p ( x b ) , p ( x a ∣ x b ) , p ( x b ∣ x a ) p(x_a),p(x_b),p(x_a|x_b),p(x_b|x_a) p(xa),p(xb),p(xa∣xb),p(xb∣xa) 。
求边缘概率
构造 x a = ( I m × m O m × n ) ( x a x b ) x_a=\begin{pmatrix}{I}_{m\times m}&{O}_{m\times n}\end{pmatrix}\begin{pmatrix}x_a\\x_b\end{pmatrix} xa=(Im×mOm×n)(xaxb),其中 I , O I,O I,O 分别是单位矩阵和零矩阵,代入上述定理中得到:
E [ x a ] = ( I O ) ( μ a μ b ) = μ a V a r [ x a ] = ( I O ) ( Σ a a Σ a b Σ b a Σ b b ) ( I O ) = Σ a a \mathbb{E}[x_a]=\begin{pmatrix}{I}&{O}\end{pmatrix}\begin{pmatrix}\mu_a\\\mu_b\end{pmatrix}=\mu_a\\ Var[x_a]=\begin{pmatrix}{I}&{O}\end{pmatrix}\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}\begin{pmatrix}{I}\\{O}\end{pmatrix}=\Sigma_{aa} E[xa]=(IO)(μaμb)=μaVar[xa]=(IO)(ΣaaΣbaΣabΣbb)(IO)=Σaa
所以 x a ∼ N ( μ a , Σ a a ) x_a\sim\mathcal{N}(\mu_a,\Sigma_{aa}) xa∼N(μa,Σaa), 边缘概率 p ( x a ) p(x_a) p(xa) 就得到了。 类似的, x b ∼ N ( μ b , Σ b b ) x_b\sim\mathcal{N}(\mu_b,\Sigma_{bb}) xb∼N(μb,Σbb)。
求条件概率
对于两个条件概率,通常是用配方法(如 PRML 的证明),这里我们用一种构造法。首先引入三个量,令:
x b ⋅ a = x b − Σ b a Σ a a − 1 x a μ b ⋅ a = μ b − Σ b a Σ a a − 1 μ a Σ b b ⋅ a = Σ b b − Σ b a Σ a a − 1 Σ a b x_{b\cdot a}=x_b-\Sigma_{ba}\Sigma_{aa}^{-1}x_a\\ \mu_{b\cdot a}=\mu_b-\Sigma_{ba}\Sigma_{aa}^{-1}\mu_a\\ \Sigma_{bb\cdot a}=\Sigma_{bb}-\Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab} xb⋅a=xb−ΣbaΣaa−1xaμb⋅a=μb−ΣbaΣaa−1μaΣbb⋅a=Σbb−ΣbaΣaa−1Σab
特别的,最后一个式子叫做 Σ a a \Sigma_{aa} Σaa 的 Schur Complementary。可以看到:
x b ⋅ a = ( − Σ b a Σ a a − 1 I n × n ) ( x a x b ) x_{b\cdot a}=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&{I}_{n\times n}\end{pmatrix}\begin{pmatrix}x_a\\x_b\end{pmatrix} xb⋅a=(−ΣbaΣaa−1In×n)(xaxb)
再有定理,有:
E [ x b ⋅ a ] = ( − Σ b a Σ a a − 1 I n × n ) ( μ a μ b ) = μ b ⋅ a V a r [ x b ⋅ a ] = ( − Σ b a Σ a a − 1 I n × n ) ( Σ a a Σ a b Σ b a Σ b b ) ( − Σ a a − 1 Σ b a T I n × n ) = Σ b b ⋅ a \mathbb{E}[x_{b\cdot a}]=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&\mathbb{I}_{n\times n}\end{pmatrix}\begin{pmatrix}\mu_a\\\mu_b\end{pmatrix}=\mu_{b\cdot a}\\ Var[x_{b\cdot a}]=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&{I}_{n\times n}\end{pmatrix}\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}\begin{pmatrix}-\Sigma_{aa}^{-1}\Sigma_{ba}^T\\{I}_{n\times n}\end{pmatrix}=\Sigma_{bb\cdot a} E[xb⋅a]=(−ΣbaΣaa−1In×n)(μaμb)=μb⋅aVar[xb⋅a]=(−ΣbaΣaa−1In×n)(ΣaaΣbaΣabΣbb)(−Σaa−1ΣbaTIn×n)=Σbb⋅a
则对于我们构造的 x b ⋅ a x_{b\cdot a} xb⋅a ,有 x b ⋅ a ∼ N ( μ b ⋅ a , Σ b b ⋅ a ) x_{b\cdot a}\sim\mathcal{N}(\mu_{b\cdot a},\Sigma_{bb\cdot a}) xb⋅a∼N(μb⋅a,Σbb⋅a) ,这里可以看到最初这个构造的设计中,核心的构造就是 x b ⋅ a x_{b\cdot a} xb⋅a ,而 μ b ⋅ a , Σ b b ⋅ a \mu_{b\cdot a},\ \Sigma_{bb\cdot a} μb⋅a, Σbb⋅a 只是两个记号,在这种构造的推导下来表示一下均值和方差。
由我们最初的构造,有 x b = x b ⋅ a + Σ b a Σ a a − 1 x a x_b=x_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}x_a xb=xb⋅a+ΣbaΣaa−1xa ,再由定理:
E [ x b ∣ x a ] = μ b ⋅ a + Σ b a Σ a a − 1 x a \mathbb{E}[x_b|x_a]=\mu_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}x_a E[xb∣xa]=μb⋅a+ΣbaΣaa−1xa
V a r [ x b ∣ x a ] = Σ b b ⋅ a Var[x_b|x_a]=\Sigma_{bb\cdot a} Var[xb∣xa]=Σbb⋅a
所以 x b ∣ x a ∼ N ( μ b ⋅ a + Σ b a Σ a a − 1 x a , Σ b b ⋅ a ) x_b|x_a\sim \mathcal{N}(\mu_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}x_a,\Sigma_{bb\cdot a}) xb∣xa∼N(μb⋅a+ΣbaΣaa−1xa,Σbb⋅a) 。类似的, x a ∣ x b ∼ N ( μ a ⋅ b + Σ a b Σ b b − 1 x b , Σ a a ⋅ b ) x_a|x_b\sim \mathcal{N}(\mu_{a\cdot b}+\Sigma_{ab}\Sigma_{bb}^{-1}x_b,\Sigma_{aa\cdot b}) xa∣xb∼N(μa⋅b+ΣabΣbb−1xb,Σaa⋅b) 。
根据边缘概率和条件概率求联合概率
已知: p ( x ) = N ( μ , Λ − 1 ) , p ( y ∣ x ) = N ( A x + b , L − 1 ) p(x)=\mathcal{N}(\mu,\Lambda^{-1}),p(y|x)=\mathcal{N}(Ax+b,L^{-1}) p(x)=N(μ,Λ−1),p(y∣x)=N(Ax+b,L−1),求解: p ( y ) , p ( x ∣ y ) p(y),p(x|y) p(y),p(x∣y)。
- 这种类型的问题在线性高斯模型、PCA降维等机器学习模型中经常出现。
- 这里的 Λ , L \Lambda, L Λ,L 称为精度矩阵,它们是协方差矩阵的逆。
解:令 y = A x + b + ϵ , ϵ ∼ N ( 0 , L − 1 ) y=Ax+b+\epsilon,\epsilon\sim\mathcal{N}(0,L^{-1}) y=Ax+b+ϵ,ϵ∼N(0,L−1),且 ϵ \epsilon ϵ 与 x x x 相互独立,还是根据上节的定理,有
E [ y ] = E [ A x + b + ϵ ] = A μ + b V a r [ y ] = A Λ − 1 A T + L − 1 \mathbb{E}[y]=\mathbb{E}[Ax+b+\epsilon]=A\mu+b\\ Var[y]=A \Lambda^{-1}A^T+L^{-1} E[y]=E[Ax+b+ϵ]=Aμ+bVar[y]=AΛ−1AT+L−1
此时,就已经得到 y ∼ N ( A μ + b , L − 1 + A Λ − 1 A T ) y\sim\mathcal{N}(A\mu+b,L^{-1}+A\Lambda^{-1}A^T) y∼N(Aμ+b,L−1+AΛ−1AT) ,即 p ( y ) = N ( A μ + b , L − 1 + A Λ − 1 A T ) p(y)=\mathcal{N}(A\mu+b,L^{-1}+A\Lambda^{-1}A^T) p(y)=N(Aμ+b,L−1+AΛ−1AT) 。
因此:
V a r [ x ∣ y ] = Λ − 1 − Λ − 1 A T ( L − 1 + A Λ − 1 A T ) − 1 A Λ − 1 Var[x|y]=\Lambda^{-1}-\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}A\Lambda^{-1} Var[x∣y]=Λ−1−Λ−1AT(L−1+AΛ−1AT)−1AΛ−1
引入 z = ( x y ) z=\begin{pmatrix}x\\y\end{pmatrix} z=(xy),我们可以得到 C o v [ x , y ] = E [ ( x − E [ x ] ) ( y − E [ y ] ) T ] Cov[x,y]=\mathbb{E}[(x-\mathbb{E}[x])(y-\mathbb{E}[y])^T] Cov[x,y]=E[(x−E[x])(y−E[y])T]。对于这个协方差可以直接计算:
C o v ( x , y ) = E [ ( x − E [ x ] ) ( y − E [ y ] ) T ] = E [ ( x − μ ) ( A x + b − A μ − b + ϵ ) T ] = E [ ( x − μ ) ( A x − A μ + ϵ ) T ] = E [ ( x − μ ) ( A x − A μ ) T + ( x − μ ) ϵ T ] = E [ ( x − μ ) ( A x − A μ ) T ] = E [ ( x − μ ) ( x − μ ) T ] A T = V a r [ x ] A T = Λ − 1 A T \begin{align} Cov(x,y)&=\mathbb{E}[(x-\mathbb{E}[x])(y-\mathbb{E}[y])^T]\\ &=\mathbb{E}[(x-\mu)(Ax+b-A\mu-b+\epsilon)^T]\\ &=\mathbb{E}[(x-\mu)(Ax-A\mu+\epsilon)^T]\\ &=\mathbb{E}[(x-\mu)(Ax-A\mu)^T+(x-\mu)\epsilon^T]\\ &=\mathbb{E}[(x-\mu)(Ax-A\mu)^T]\\ &=\mathbb{E}[(x-\mu)(x-\mu)^T]A^T\\ &=Var[x]A^T\\ &=\Lambda^{-1}A^T \end{align} Cov(x,y)=E[(x−E[x])(y−E[y])T]=E[(x−μ)(Ax+b−Aμ−b+ϵ)T]=E[(x−μ)(Ax−Aμ+ϵ)T]=E[(x−μ)(Ax−Aμ)T+(x−μ)ϵT]=E[(x−μ)(Ax−Aμ)T]=E[(x−μ)(x−μ)T]AT=Var[x]AT=Λ−1AT
注意到协方差矩阵的对称性,所以 p ( z ) = N ( μ A μ + b ) , ( Λ − 1 Λ − 1 A T A Λ − 1 L − 1 + A Λ − 1 A T ) ) p(z)=\mathcal{N}\begin{pmatrix}\mu\\A\mu+b\end{pmatrix},\begin{pmatrix}\Lambda^{-1}&\Lambda^{-1}A^T\\A\Lambda^{-1}&L^{-1}+A\Lambda^{-1}A^T\end{pmatrix}) p(z)=N(μAμ+b),(Λ−1AΛ−1Λ−1ATL−1+AΛ−1AT))。根据上一节的公式,我们可以得到:
E [ x ∣ y ] = μ + Λ − 1 A T ( L − 1 + A Λ − 1 A T ) − 1 ( y − A μ − b ) \mathbb{E}[x|y]=\mu+\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}(y-A\mu-b) E[x∣y]=μ+Λ−1AT(L−1+AΛ−1AT)−1(y−Aμ−b)
V a r [ x ∣ y ] = Λ − 1 − Λ − 1 A T ( L − 1 + A Λ − 1 A T ) − 1 A Λ − 1 Var[x|y]=\Lambda^{-1}-\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}A\Lambda^{-1} Var[x∣y]=Λ−1−Λ−1AT(L−1+AΛ−1AT)−1AΛ−1
故得到: p ( x ∣ y ) = N ( μ + Λ − 1 A T ( L − 1 + A Λ − 1 A T ) − 1 ( y − A μ − b ) , Λ − 1 − Λ − 1 A T ( L − 1 + A Λ − 1 A T ) − 1 A Λ − 1 ) p(x|y)=\mathcal{N}(\mu+\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}(y-A\mu-b),\Lambda^{-1}-\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}A\Lambda^{-1}) p(x∣y)=N(μ+Λ−1AT(L−1+AΛ−1AT)−1(y−Aμ−b),Λ−1−Λ−1AT(L−1+AΛ−1AT)−1AΛ−1) 。
Ref
边栏推荐
- leetcode:1562. 查找大小为 M 的最新分组【模拟 + 端点记录 + 范围合并】
- WindowInsetsControllerCompat简单使用
- Exam preparation plan
- Qlib quantitative source analysis: qlib/qlib/contrib/model/GBDT py
- Carefully organize 16 MySQL usage specifications to reduce problems by 80% and recommend sharing with the team
- 2022年最新重庆建筑八大员(电气施工员)模拟题库及答案
- NIO编程
- GDB 源码分析系列文章五:动态库延迟断点实现机制
- Matlab/ArcGIS processing GPM global monthly precipitation data
- MYSQL经典面试题
猜你喜欢
leetcode: 1648. Color ball with decreasing sales value [Boundary find by two points]
RTL8762DK 点灯/LED(三)
【历史上的今天】7 月 31 日:“缸中之脑”的提出者诞生;Wi-Fi 之父出生;USB 3.1 标准发布
SC7A20(士兰微-加速度传感器)示例
Kyoto University: Masaki Waga | Dynamic Masking for Reinforcement Learning in Black Box Environments
TCP协议详解
zeno使用方法笔记
RTL8762DK UART (two)
Matlab/Arcgis processing nc data
MYSQL逻辑架构
随机推荐
从零造键盘的键盘超级喜欢,IT人最爱
Pylint检查规则中文版
TFC CTF 2022 WEB Diamand WriteUp
Force buckle 2326, 197
NIO编程
一体化步进电机在无人机自动机场的应用
Carefully summarize thirteen suggestions to help you create more suitable MySQL indexes
YOLO怎么入门?怎么实现自己的训练集?
MVCC总结
Nmap Operation Manual - Full Version
Carefully summarize thirteen suggestions to help you create more suitable MySQL indexes
Luogu P3373: Segment tree
二叉树遍历非递归程序 -- 使用栈模拟系统栈
虚继承的原理
Application of integrated stepper motor in UAV automatic airport
继承和友元,静态成员的关系
Web3.0:构建 NFT 市场(一)
机器学习应该如何入门?
MYSQL-批量插入数据
pycaret source code analysis: download dataset\Lib\site-packages\pycaret\datasets.py