当前位置:网站首页>高维高斯分布基础
高维高斯分布基础
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--打家劫舍问题
- 【历史上的今天】7 月 31 日:“缸中之脑”的提出者诞生;Wi-Fi 之父出生;USB 3.1 标准发布
- 南方科技大学:Xiaoying Tang | AADG:视网膜图像分割领域泛化的自动增强
- ECCV2022 Workshop | Multi-Object Tracking and Segmentation in Complex Environments
- 精心总结十三条建议,帮你创建更合适的MySQL索引
- Item 36: Specify std::launch::async if asynchronicity is essential.
- WAASAP WebClient UI页面标签的决定逻辑介绍
- 蓝图:杨辉三角排列
- leetcode:1648. 销售价值减少的颜色球【二分找边界】
- 继承的注意事项
猜你喜欢
RTL8762DK PWM(七)
机器学习初学者可以学哪些实战项目?
One line of code to solve CoreData managed object properties change in SwiftUI problem of animation effects
TFC CTF 2022 WEB Diamand WriteUp
欧拉系统(euleros):升级Mysql
MYSQL逻辑架构
南方科技大学:Xiaoying Tang | AADG:视网膜图像分割领域泛化的自动增强
现代企业架构框架1
Southern University of Science and Technology: Xiaoying Tang | AADG: Automatic Enhancement for Generalization in the Field of Retinal Image Segmentation
Modern Enterprise Architecture Framework 1
随机推荐
机器学习初学者可以学哪些实战项目?
LeetCode--The problem of robbery
[AMEX] LGBM Optuna American Express Credit Card Fraud Contest kaggle
RTL8762DK RTC(五)
【Cryptography/Cryptanalysis】Cryptanalysis method based on TMTO
Item 36: Specify std::launch::async if asynchronicity is essential.
Modern Enterprise Architecture Framework 1
In 2022, the latest eight Chongqing construction members (electrical construction workers) simulation question bank and answers
YOLO怎么入门?怎么实现自己的训练集?
OSF一分钟了解敏捷开发模式
Nmap 操作手册 - 完整版
RTL8762DK WDG (six)
WAASAP WebClient UI页面标签的决定逻辑介绍
[AMEX] LGBM Optuna美国运通信用卡欺诈赛 kaggle
Carefully organize 16 MySQL usage specifications to reduce problems by 80% and recommend sharing with the team
开源好用的 流程图绘制工具 drawio
OSD read SAP CRM One Order application log way of optimization
WebApi 打个Attribute 统一处理异常
【密码学/密码分析】基于TMTO的密码分析方法
Likou Binary Tree