Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者: 周翔宇 ( 中山大学 )
邮箱: zhouxy98@mail2.sysu.edu.cn
目录
本文基于 Vincent Spruyt 推文 A geometric interpretation of the covariance matrix 的内容,结合相关文献介绍了协方差矩阵的几何意义。我们会通过数据的线性变换来从几何意义的角度直观理解协方差矩阵。
具体来说,我们将首先回顾如何描述一组一维数据,随后介绍如何用协方差矩阵描述二维乃至多维数据。过程中适当的地方会用Stata作图举例。为了尽量兼顾数学上的严谨性,涉及到的若干证明我会附在文中供大家参考。
协方差矩阵 ( Covariance Matrix ) 可以很好地刻画随机变量之间的相关性,在金融工程、统计和机器学习等领域都有非常广泛的应用。
协方差矩阵的特征值分解在计量经济学和具体的实证研究中也应用广泛,其中主成分分析(Principal Component Analysis)极具代表性。实证研究中,变量数过多会给后续的处理造成极大阻碍。有时我们拿到的数据集包含成百上千个变量,但其中绝大多数变量常常并没有对数据集的信息做出太多贡献。面对这种情况,PCA 可以通过正交变换“合成”少数几个线性不相关的新变量,使得它们能够描述原数据中的绝大部分信息。这样我们就可以舍弃“低效率”的原数据集,从而更高效地解决问题。
简单来讲,PCA 的步骤如下:
了解了协方差矩阵特征值分解的巨大威力后,我们从回顾方差开始本文的第二部分。
对于一组给定的数据,该怎么描述它的“形状”?这个问题对于一维数据并不难回答:只要我们知道了样本均值和方差,就可以对它的分布有一个大体的认知。下面我们从方差开始。
方差是一组一维数据离散程度的度量。若样本量为
以服从一维正态分布的数据为例:
Stata代码如下:
//Fig01:不同标准差的正态分布密度函数图
twoway function y=normalden(x), range(-5 5) ///
|| function y=normalden(x,0.8), ///
range(-5 5) lpattern(dash) ///
|| function y=normalden(x,1.2), ///
range(-5 5) 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 一般越小。
但是,如果沿用一维数据的分析思路,计算水平方向的方差
生成上图利用了对白数据进行线性变换的方法。相关细节会在第 4 节展开
Stata代码如下:
//Fig02:生成一组二维数据并绘图
clear
set obs 1000 //设置观察值N=1000
set seed 2022 //设置随机数种子,令所生成数据可重现
gen x = rnormal() //生成服从标准正态分布的一组x
set obs 1000 //设置观察值N=1000
set seed 1204 //设置seed,令所生成数据可重现
gen y = rnormal() //生成服从标准正态分布的一组y
set matsize 5000 //设置矩阵尺寸以便存入数据
mkmat x y, matrix(data) //将x、y作为列向量存入矩阵data
matrix 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/4
matrix newdata1T = newdata1' //取处理过的newdata1的转置,命名为newdata1T
svmat newdata1T, names(col) //将newdata1T中的列提出,作为新的变量r1,r2
rename r1 x1 //将提取出来的变量r1重命名为x1
rename r2 y1 //将提取出来的变量r2重命名为y1
twoway scatter y1 x1, xlabel(-16 -8 0 8 16) msize(tiny) //作图
书接上回,当遇到方差无法很好地描述的高维数据时,协方差矩阵就派上用场了。
在统计学与概率论中,协方差矩阵是一个矩阵,其
假设
且
那么,协方差矩阵的第
而协方差矩阵为:
即:
又因为:
对于经过标准化的二维数据集:
易证明其协方差矩阵
这意味着它必定是一个对称矩阵,非对角线上的是协方差,对角线上的元素是相应随机变量与自己的协方差,也就是其方差。
Stata代码如下 ( 数据线性变换工作与上图代码类似。为保持简洁,这里只展示作图部分 ) :
*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的格式呈现在同一张图中
四幅图所对应数据的协方差矩阵分别为:
从上图可以直观地看到,协方差矩阵对角线上的元素描述了数据沿着坐标轴方向的离散情况,非对角线上的元素则描述了数据在对角线方向上的分散情况。
现在我们已经知道,协方差矩阵定义了数据在不同方向上的离散情况,但是协方差矩阵本身似乎还不够直观。以上图[1]为例,数据明显朝着左上方和右下方延伸,我们能否通过它的协方差矩阵找到某个指向数据最大分散方向的向量,且其模(长度)等于数据在该方向上的“方差”呢?也就是如何从原始数据
要实现上面的想法,可以把数据降维处理。以上图[1]为例,问题可以转化为:找到某个单位向量
我们将二维数据
其中
此时
使新数据
将(1)式代入上式可得:
现在,问题转化成了在
对
神奇的事情发生了,协方差矩阵
现在我们知道:协方差矩阵的最大特征向量总是指向数据“最大方差”的方向,且它的模等于相应的最大特征值,第二大特征向量(对于
上图中数据集的协方差矩阵为:
对应的最大、最小特征向量分别为:
在图中标出两个模为对应特征值的特征向量,可以直观地感受到“指向最大方差方向”和“指向最小方差方向”的含义。
注意到
Stata代码如下:
*Fig04: 在图中标出两个特征向量
set matsize 5000 //设置矩阵尺寸以便存入数据
mkmat x y, matrix(data) //将x、y作为列向量存入矩阵data
matrix 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/4
matrix newdata1T = newdata1' //取处理过的newdata1的转置,命名为newdata1T
svmat newdata1T, names(col) //将newdata1T中的列提出,作为新的变量r1,r2
rename r1 x1 //将提取出来的变量r1重命名为x1
rename 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(vtiny)) ///
(line v_max a, lwidth(thick) lcolor(orange)) ///
(line v_min b , lwidth(thick) lcolor(green))
接下来让我们进入较为技术化的部分,探究对白数据 (White data) 进行的线性变换矩阵与新数据的协方差矩阵之间的关系。
首先介绍 ( 数据 ) 矩阵的线性变换:用Stata生成一组二维白数据
Stata代码如下:
//Fig05:生成一组二维白数据
clear
set obs 1000 //设置观察值N=1000
set seed 2022 //设置随机数种子,令所生成数据可重现
gen x = rnormal() //生成服从标准正态分布的一组x
set obs 1000 //设置观察值N=1000
set seed 1204 //设置随机数种子,令所生成数据可重现
gen y = rnormal() //生成服从标准正态分布的一组y
twoway scatter y x, msize(tiny)
可以直观地看到,
若
线性变换
考察线性变换对协方差矩阵的影响:白数据
新数据
即放缩矩阵
例如:
此时,
如图,用Stata实现该线性变换,可以看到,数据在
Stata代码如下:
//Fig06:缩放效果图
clear
set obs 1000 //设置观察值N=1000
set seed 2022 //设置seed,令所生成数据可重现
gen x = rnormal() //生成服从标准正态分布的一组x
set obs 1000 //设置观察值N=1000
set seed 1204 //设置seed,令所生成数据可重现
gen y = rnormal() //生成服从标准正态分布的一组y
set matsize 5000 //设置矩阵尺寸以便存入数据
mkmat x y, matrix(data) //将x、y作为列向量存入矩阵data
matrix whitedata = data' //转置data得到矩阵whitedata以便后续操作
matrix S = (4,0\0,1) //输入尺度矩阵
matrix T = S
matrix newdata1 = T*whitedata //进行线性变换
matrix newdata1T = newdata1' //取处理过的newdata1的转置,命名为newdata1T
svmat newdata1T, names(col) //将newdata1T中的列提出,作为新的变量r1,r2
drop x y
rename r1 x //将提取出来的变量r1重命名为x
rename r2 y //将提取出来的变量r2重命名为y
twoway scatter y x, xlabel(-16 -8 0 8 16) ///
ylabel(-10 -5 0 5 10) msize(tiny) //画出散点图
若
线性变换
例如:
现在对上面通过拉长得到的
用Stata实现该线性变换,新数据
Stata代码如下:
//Fig07
clear
set obs 1000 //设置观察值N=1000
set seed 2022 //设置seed,令所生成数据可重现
gen x = rnormal() //生成服从标准正态分布的一组x
set obs 1000 //设置观察值N=1000
set seed 1204 //设置seed,令所生成数据可重现
gen y = rnormal() //生成服从标准正态分布的一组y
*白数据旋转缩放后绘图
set matsize 5000 //设置矩阵尺寸以便存入数据
mkmat x y, matrix(data) //将x、y作为列向量存入矩阵data
matrix 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/3
matrix newdata1T = newdata1' //取处理过的newdata1的转置,命名为newdata1T
svmat newdata1T, names(col) //将newdata1T中的列提出,作为新的变量r1,r2
drop x y
rename r1 x //将提取出来的变量r1重命名为x
rename r2 y //将提取出来的变量r2重命名为y
twoway scatter y x, xlabel(-16 -8 0 8 16) msize(tiny)
对于放缩矩阵,我们得到了式
对于由白数据
我们取它的单位化后的特征向量
其中
保留
与
注意到我们在上一小节中已经证明了协方差矩阵的最大特征向量指向方差最大的方向,且不同的特征向量两两正交。因为这里的
即:
另一方面,我们在上一小节中已经证明了协方差矩阵的特征值为数据在对应的特征向量方向上的方差值,由式
整体的线性变换矩阵为:
由
这样我们就得到了线性变换矩阵
最后,让我们以一个直观的例子结束这部分的介绍。
上图中数据集是由白数据通过转换矩阵
根据上面的推导,对应的协方差矩阵:
对应的最大、最小特征向量分别为:
在图中标出两个模为对应特征值的特征向量(为了显示效果,将两个向量同时向反方向进行了延长处理)
本篇推文中,我们从几何意义的角度重新认识了协方差矩阵,推导了数据“形状”与协方差矩阵特征向量、特征值之间的关系。
同时,我们得到了线性变换矩阵
具体地说,此线性变换矩阵
Note:产生如下推文列表的 Stata 命令为:
lianxh 方差 主成分
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh