主成分分析-交互固定效应基础:协方差矩阵的几何意义

发布时间:2022-01-08 阅读 258

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

   

作者: 周翔宇 ( 中山大学 )
E-Mail: zhouxy98@mail2.sysu.edu.cn

备注: 本文翻译自 Vincent Spruyt 推文A geometric interpretation of the covariance matrix,并结合了若干相关文献介绍协方差矩阵的几何意义。特此致谢!


目录
[toc]


1. 引言

1.1 你将会在本篇推文中看到什么?

我们会通过数据的线性变换来从几何意义的角度直观理解协方差矩阵。具体来说,我们将首先回顾如何描述一组一维数据,随后介绍如何用协方差矩阵描述二维乃至多维数据。过程中适当的地方会用Stata作图举例。为了尽量兼顾数学上的严谨性,涉及到的若干证明我会附在文中供大家参考。

1.2 协方差矩阵的特征值分解有什么用?

协方差矩阵 ( Covariance Matrix ) 可以很好地刻画随机变量之间的相关性,在金融工程、统计和机器学习等领域都有非常广泛的应用。

协方差矩阵的特征值分解在计量经济学和具体的实证研究中也应用广泛,其中主成分分析(Principal Component Analysis)极具代表性。实证研究中,变量数过多会给后续的处理造成极大阻碍。有时我们拿到的数据集包含成百上千个变量,但其中绝大多数变量常常并没有对数据集的信息做出太多贡献。面对这种情况,PCA 可以通过正交变换“合成”少数几个线性不相关的新变量,使得它们能够描述原数据中的绝大部分信息。这样我们就可以舍弃“低效率”的原数据集,从而更高效地解决问题。

简单来讲,PCA 的步骤如下:

  1. 计算原数据集的协方差矩阵
  2. 计算该协方差矩阵的特征向量和特征值
  3. 每个特征向量的特征值表示沿给定特征向量所解释的方差量,将所有特征向量按“解释能力”从大到小排序,即按照特征值的降序排列特征向量
  4. 根据我们想要在模型中“捕捉”多少方差来选择保留我们感兴趣的前2-3个或更多的特征向量及对应的特征值(也就是前2-3个主成分),同时舍弃余下全部特征值
  5. 用主成分分析构造回归模型时,把保留下来的各主成分作为新自变量代替原来的自变量做回归分析。这样,我们以原数据集的极少一部分信息为代价,换来了描述效率的极大提高

了解了协方差矩阵特征值分解的巨大威力后,我们从回顾方差开始本文的第二部分。

2. 从方差开始

对于一组给定的数据,该怎么描述它的“形状”?这个问题对于一维数据并不难回答:只要我们知道了样本均值和方差,就可以对它的分布有一个大体的认知。下面我们从方差开始。

方差是一组一维数据离散程度的度量。若样本量为 N , 均值为 x¯ , 则样本方差的无偏估计为:

以服从一维正态分布的数据为例:

正态分布密度函数图 Stata代码如下:

  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
//Fig01:不同标准差的正态分布密度函数图twoway function y=normalden(x), range(-4 4)  /// 	|| function y=normalden(x,0.8),  ///		range(-4 4) lpattern(dash)  ///	|| function y=normalden(x,1.2),  ///		range(-4 4) lpattern(dash_dot)  ///	|| ,title("Normal Densities with different variances")  ///		legend(label(1 "N(0,1^2)") label(2 "N(0,0.8^2)")  ///		label(3 "N(0,1.2^2)"))              

可以看到,方差很好地描述了一维数据的特征:方差越大,数据分布越分散;方差越小,数据分布越集中!

但是,拿到这样一组二维数据时,方差的局限便显现出来了: 直观上看,图中的点呈现出长条形,整体向左上方和右下方延伸,而在右上和左下这个方向上数据较为集中。这意味着 x 和 y 存在比较明显的负相关关系:x 越大,y 一般越小。

但是,如果沿用一维数据的分析思路,计算水平方向的方差 σx2 和垂直方向的方差 σy2 只能描述数据沿着坐标轴方向的离散程度,并不能很好地体现两个分量的相关关系。看来我们必须找到某种新的方法才能更好地描述这组数据。

一组二维数据
一组二维数据

生成上图利用了对白数据进行线性变换的方法。相关细节会在第 4 节展开

Stata代码如下:

  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
//Fig02:生成一组二维数据并绘图clearset obs 1000  //设置观察值N=1000set seed 2022  //设置随机数种子,令所生成数据可重现gen x = rnormal()  //生成服从标准正态分布的一组xset obs 1000  //设置观察值N=1000set seed 1204  //设置seed,令所生成数据可重现gen y = rnormal()  //生成服从标准正态分布的一组y
set matsize 5000 //设置矩阵尺寸以便存入数据mkmat x y, matrix(data) //将x、y作为列向量存入矩阵datamatrix whitedata = data' //转置data得到矩阵whitedata以便后续操作
matrix S = (1,0\0,3) //输入尺度矩阵matrix R = (cos(_pi/4), -sin(_pi/4)\sin(_pi/4),cos(_pi/4)) //输入旋转矩阵matrix T = R*S //得到线性变化矩阵T
matrix newdata1 = T*whitedata //将x、y分别放大后整体逆时针旋转pi/4matrix newdata1T = newdata1' //取处理过的newdata1的转置,命名为newdata1Tsvmat newdata1T, names(col) //将newdata1T中的列提出,作为新的变量r1,r2
rename r1 x1 //将提取出来的变量r1重命名为x1rename r2 y1 //将提取出来的变量r2重命名为y1twoway scatter y1 x1, xlabel(-16 -8 0 8 16) msize(tiny) //作图

3.描述高维数据

书接上回,当遇到方差无法很好地描述的高维数据时,协方差矩阵就派上用场了。

3.1 协方差矩阵的定义

在统计学与概率论中,协方差矩阵是一个矩阵,其 (i,j) 位置的元素是第 i 个与第 j 个随机变量之间的协方差。这是从标量随机变量到高维度随机向量的自然推广。

假设 X 是由 2 个随机变量组成的列向量 ,

且 μi 是 Xi 的期望值,即:

那么,协方差矩阵的第 (i,j) 项被定义为如下形式:

而协方差矩阵为:

即:

又因为:

对于经过标准化的二维数据集:

易证明其协方差矩阵

这意味着它必定是一个对称矩阵,非对角线上的是协方差,对角线上的元素是相应随机变量与自己的协方差,也就是其方差。

3.2 几个例子

几个例子 Stata代码如下 ( 数据线性变换工作与上图代码类似。为保持简洁,这里只展示作图部分 ) :

  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
*snip*twoway scatter y1 x1, xlabel(-16 -8 0 8 16) ylabel(-10 -5 0 5 10) msize(tiny) saving(pic1)  //画出第一幅散点图并保存,下同twoway scatter y2 x2, xlabel(-16 -8 0 8 16) ylabel(-10 -5 0 5 10) msize(tiny) saving(pic2)twoway scatter y3 x3, xlabel(-16 -8 0 8 16) ylabel(-10 -5 0 5 10) msize(tiny) saving(pic3)twoway scatter y4 x4, xlabel(-16 -8 0 8 16) ylabel(-10 -5 0 5 10) msize(tiny) saving(pic4)graph combine pic1.gph pic2.gph pic3.gph pic4.gph, imargin(small) col(2)  //将四幅图以2*2的格式呈现在同一张图中

四幅图所对应数据的协方差矩阵分别为:

从上图可以直观地看到,协方差矩阵对角线上的元素描述了数据沿着坐标轴方向的离散情况,非对角线上的元素则描述了数据在对角线方向上的分散情况。

4. 寻找“最大方差”的方向

现在我们已经知道,协方差矩阵定义了数据在不同方向上的离散情况,但是协方差矩阵本身似乎还不够直观。以上图[1]为例,数据明显朝着左上方和右下方延伸,我们能否通过它的协方差矩阵找到某个指向数据最大分散方向的向量,且其模(长度)等于数据在该方向上的“方差”呢?也就是如何从原始数据 D 的协方差矩阵 Σ 得到这个向量?

要实现上面的想法,可以把数据降维处理。以上图[1]为例,问题可以转化为:找到某个单位向量 v=[a,b] ,使原始数据沿着 v 的方向分布最分散。

4.1 数据的降维处理

我们将二维数据 D 投影到这个一维单位向量 v 上,用得到的投影值构成一组新的一维数据 D1 。也就是将 v 与数据集 D 中的每个数据求内积得到一组新的数据 D1 ,即:

其中 D 为按列排列的数据集:

此时 D1 为一组一维数据,其中的元素是原始数据在 v 方向上的投影值:

4.2 求最大方差与特征向量

使新数据 D1 方差最大的 v 就是所求。本例中, D 是经过标准化的二维数据,新数据 D1 的方差为:

将(1)式代入上式可得:

现在,问题转化成了在 vv=1 的约束下,求 Vn 的极值。下面用拉格朗日乘子法来求解(参考GitHub:瑞利商与极值计算),定义拉格朗日函数:

对 v 求梯度,并令其值为0:

神奇的事情发生了,协方差矩阵 Σ特征向量 vi正好能使 Vn 取得极值,且极值为对应的特征值 λi !同时,由于 Σ=Σ ,可以证明其特征值不等的特征向量两两正交。 (证明参见 知乎:对称矩阵特征向量正交)

现在我们知道:协方差矩阵的最大特征向量总是指向数据“最大方差”的方向,且它的模等于相应的最大特征值,第二大特征向量(对于 2×2 的协方差矩阵,第二大特征向量即是其最小特征向量)总是指向数据“第二大方差”的方向,且与第一大特征向量正交。

最大特征向量
最大特征向量

上图中数据集的协方差矩阵为:

对应的最大、最小特征向量分别为:

在图中标出两个模为对应特征值的特征向量,可以直观地感受到“指向最大方差方向”和“指向最小方差方向”的含义。

注意到 vmax 与 vmax 垂直

Stata代码如下:

  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line

//Fig04:生成一组二维数据clearset obs 1000 //设置观察值N=1000set seed 2022 //设置seed,令所生成数据可重现gen x = rnormal() //生成服从标准正态分布的一组x
set obs 1000 //设置观察值N=1000set seed 1204 //设置seed,令所生成数据可重现gen y = rnormal() //生成服从标准正态分布的一组y
*白数据旋转缩放后绘图set matsize 5000 //设置矩阵尺寸以便存入数据mkmat x y, matrix(data) //将x、y作为列向量存入矩阵datamatrix whitedata = data' //转置data得到矩阵whitedata以便后续操作
matrix S = (1,0\0,3) //输入尺度矩阵matrix R = (cos(_pi/4), -sin(_pi/4)\sin(_pi/4),cos(_pi/4)) //输入旋转矩阵
matrix T = R*S
matrix newdata1 = T*whitedata //将x、y分别放大后整体逆时针旋转pi/4matrix newdata1T = newdata1' //取处理过的newdata1的转置,命名为newdata1Tsvmat newdata1T, names(col) //将newdata1T中的列提出,作为新的变量r1,r2
rename r1 x1 //将提取出来的变量r1重命名为x1rename r2 y1 //将提取出来的变量r2重命名为y1
gen a = -9*2^0.5/2 in 1 //特征向量坐标replace a = 9*2^0.5/2 in 2
gen v_max = 9*2^0.5/2 in 1 //特征向量坐标replace v_max = -9*2^0.5/2 in 2
gen b = -2^0.5/2 in 1 //特征向量坐标replace b = 2^0.5/2 in 2
gen v_min = -2^0.5/2 in 1 //特征向量坐标replace v_min = 2^0.5/2 in 2
twoway (scatter y1 x1, msize(tiny)) ///(line v_max a, lwidth(medthick) lcolor(orange)) ///(line v_min b , lwidth(medthick) lcolor(green)),xlabel(-9 0 9) ylabel(-7 0 7)

5.协方差矩阵的特征值分解与线性变换

接下来让我们进入较为技术化的部分,探究对白数据 (White data) 进行的线性变换矩阵与新数据的协方差矩阵之间的关系。

5.1 矩阵的线性变换

首先介绍 ( 数据 ) 矩阵的线性变换:用Stata生成一组二维白数据 D, 它的协方差矩阵为2*2单位矩阵。如下图所示:

"随机生成的一组二维白数据"
"随机生成的一组二维白数据"

Stata代码如下:

  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
//Fig05:生成一组二维白数据clearset obs 1000  //设置观察值N=1000set seed 2022  //设置随机数种子,令所生成数据可重现gen x = rnormal()  //生成服从标准正态分布的一组x
set obs 1000 //设置观察值N=1000set seed 1204 //设置随机数种子,令所生成数据可重现gen y = rnormal() //生成服从标准正态分布的一组ytwoway scatter y x, msize(tiny)

可以直观地看到,D 中 x 和 y 没有明显的相关性。但是利用线性变换矩阵 T=RS ,可以将 D 变换各种各样形态各异的数据。其中 T 为两个旋转矩阵 R 和放缩矩阵 S 的乘积。下面具体来看:

线性变换 Ds=SD 意味着在 x 方向上将数据 D 变为原来的 σx 倍,在 y 方向上将数据 D 变为原来的 σy 倍。

考察线性变换对协方差矩阵的影响:白数据 D 的协方差矩阵为:

新数据 Ds 的协方差矩阵为:

即放缩矩阵 S 与新数据 Ds 的协方差矩阵 Σs 存在如下关系:

例如:

此时,Ds 的协方差矩阵为:

如图,用Stata实现该线性变换,可以看到,数据在 X 轴方向上被“拉长了”。

Stata代码如下:

  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
//Fig06:缩放效果图clearset obs 1000  //设置观察值N=1000set seed 2022  //设置seed,令所生成数据可重现gen x = rnormal()  //生成服从标准正态分布的一组x
set obs 1000 //设置观察值N=1000set seed 1204 //设置seed,令所生成数据可重现gen y = rnormal() //生成服从标准正态分布的一组yset matsize 5000 //设置矩阵尺寸以便存入数据mkmat x y, matrix(data) //将x、y作为列向量存入矩阵datamatrix whitedata = data' //转置data得到矩阵whitedata以便后续操作
matrix S = (4,0\0,1) //输入尺度矩阵matrix T = S
matrix newdata1 = T*whitedata //进行线性变换matrix newdata1T = newdata1' //取处理过的newdata1的转置,命名为newdata1Tsvmat newdata1T, names(col) //将newdata1T中的列提出,作为新的变量r1,r2drop x y
rename r1 x //将提取出来的变量r1重命名为xrename r2 y //将提取出来的变量r2重命名为y
twoway scatter y x, xlabel(-16 -8 0 8 16) /// ylabel(-10 -5 0 5 10) msize(tiny) //画出散点图

线性变换 Dr=RD 意味着将数据整体绕原点逆时针旋转 θ

例如:

现在对上面通过拉长得到的 Ds 进一步进行旋转变换:

用Stata实现该线性变换,新数据 Dnew 如图:

"缩放、旋转后的数据"
"缩放、旋转后的数据"

Stata代码如下:

  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
  • ounter(line
//Fig07clearset obs 1000  //设置观察值N=1000set seed 2022  //设置seed,令所生成数据可重现gen x = rnormal()  //生成服从标准正态分布的一组x
set obs 1000 //设置观察值N=1000set seed 1204 //设置seed,令所生成数据可重现gen y = rnormal() //生成服从标准正态分布的一组y
*白数据旋转缩放后绘图set matsize 5000 //设置矩阵尺寸以便存入数据mkmat x y, matrix(data) //将x、y作为列向量存入矩阵datamatrix whitedata = data' //转置data得到矩阵whitedata以便后续操作
matrix S = (4,0\0,1) //输入缩放矩阵matrix R = (cos(_pi/3), -sin(_pi/3)\sin(_pi/3),cos(_pi/3)) //输入旋转矩阵
matrix T = R*S
matrix newdata1 = T*whitedata //将x、y分别放大后整体逆时针旋转pi/3matrix newdata1T = newdata1' //取处理过的newdata1的转置,命名为newdata1Tsvmat newdata1T, names(col) //将newdata1T中的列提出,作为新的变量r1,r2drop x yrename r1 x //将提取出来的变量r1重命名为xrename r2 y //将提取出来的变量r2重命名为y
twoway scatter y x, xlabel(-16 -8 0 8 16) msize(tiny)

对于放缩矩阵,我们得到了式 (3) ,但是,当加入旋转矩阵 R 后,线性变换矩阵 T=RS 与新数据的协方差矩阵 Σnew 的关系如何呢?下面我们尝试运用矩阵的特征值分解将协方差矩阵 Σnew 分解为旋转矩阵和缩放矩阵的乘积。

5.2 协方差矩阵的特征值分解

对于由白数据 D 经线性变换得到的数据集:

我们取它的单位化后的特征向量 v,有:

其中 λ 为对应的特征值。

保留 (4) 中协方差矩阵 Σnew 的每一个特征向量和特征值,并将特征向量 vi 作为列向量组成一个新矩阵 V;同时将对应的特征值 λi 放入对角矩阵 L 中。具体来说,在二维数据的情形下:

与 (4) 对应,容易得到:

(5) 式两边同乘 V1,得到:

注意到我们在上一小节中已经证明了协方差矩阵的最大特征向量指向方差最大的方向,且不同的特征向量两两正交。因为这里的 v单位化后的特征向量,可以将矩阵 V 看作旋转矩阵 R。且因为 V 中列向量两两正交,有:

即:

另一方面,我们在上一小节中已经证明了协方差矩阵的特征值为数据在对应的特征向量方向上的方差值,由式 (3) 可知:L表示一个缩放矩阵 S

整体的线性变换矩阵为:

由 (6)(7)(8)(9) ,协方差矩阵 Σnew 可以继续分解为:

这样我们就得到了线性变换矩阵 T 与新数据协方差矩阵 Σnew 之间的关系!

最后,让我们以一个直观的例子结束这部分的介绍。

"一个例子" 上图中数据集是由白数据通过转换矩阵 T 得来的,其中:

根据上面的推导,对应的协方差矩阵 Σ 可以由转换矩阵 T 得到:

对应的最大、最小特征向量分别为:

在图中标出两个模为对应特征值的特征向量(为了显示效果,将两个向量同时向反方向进行了延长处理)

6. 总结

本篇推文中,我们从几何意义的角度重新认识了协方差矩阵,推导了数据“形状”与协方差矩阵特征向量、特征值之间的关系。

同时,我们得到了线性变换矩阵 T 与新数据协方差矩阵 Σnew 之间的关系: Σnew 与白数据的线性变换矩阵 T 有直接的关系。

具体地说,此线性变换矩阵 T 完全由 Σnew 的特征向量和特征值确定:单位化的特征向量组成旋转矩阵,而特征值对应于缩放矩阵中元素的平方。

7. 参考资料