Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:姚和平 (中山大学)
邮箱:yaohp@mail2.sysu.edu.cn
目录
特征值和特征向量是线性代数里十分重要的概念,在经济学、统计学、物理、化学等学科中有着广泛的应用。这篇文章将为读者细致地讲解这两个概念,并介绍其重要应用—PCA (Principal Component Analysis)。作为引入,我们先来看一个现实案例。
1990 年,环境保护学家与木材行业在砍伐太平洋西北部原始森林,是否会造成斑点猫头鹰灭绝的问题上争论不休。数学家们也因此开始了对斑点猫头鹰种群的动力学研究,R.Lamberson (2005) 及其同事建立了以每年种群数量为区间的演化模型,时间为
根据统计数据,新的幼年猫头鹰在
我们用数学公式来表达这一模型:
将上面式子改写成矩阵形式,就将该问题转换成一个形式为
上式表达除了可以用来表示猫头鹰群落的演化过程,还可以表示类似城市人口流动、传染病 SEIR 模型等,在日常生活中应用十分广泛。但是
学过线性代数的读者都知道,矩阵乘法
设
因此我们可以看到,对于一个矩阵
聪明的读者可能会想到,通过定义条件
把
因此我们可以通过解特征方程
求
因式分解后得到的上式也被称为特征方程,因此
特征值和特征向量还有一个非常有趣的性质:
这个过程也叫矩阵的对角化,此时矩阵
我们回到斑点猫头鹰的演化问题,特征值与特征向量是解决动力系统的关键。你可能会问,这种将矩阵乘法变为标量乘法的方法对于任意一个特定方阵来说,只有少数几个特征向量有效,比如在斑点猫头鹰的问题里,种群量
以斑点猫头鹰为例,这里
其中
这样我们就把一个复杂的矩阵乘法问题转换为一个简单的标量乘法问题,但这样会对运算产生什么样的影响呢?具体到斑点猫头鹰问题,我们可以计算出转移矩阵
由于所有的特征值都小于
对于求一个
事实上,Matlab 的算法就是先计算出特征值
QR 算法的核心思路是由矩阵
根据上面两条定理,我们可以很容易看出,
QR 分解是指将矩阵 qrd()
函数可以完成这一过程。
// --QR分解--
mata //进入 Mata 环境
A = ( 1, 2, 3 \ 4, 5, 6 \ 7, 8, 9 )
qrd(A, Q, R)
st_matrix("Q",Q) //将 Meta 里的矩阵转换到 Stata 环境中
st_matrix("R",R) //同上
end //退出 Mata 环境
mat list Q
mat list R
当矩阵
随着新矩阵
到目前为止,我们已经讲解清楚了特征值和特征向量,那么接下来让我们来看一下主成分分析 (Principal Component Analysis) ,它是特征值和特征向量的一个主要应用方向。
在大数据时代,我们面对的数据往往更多的特征。过多的变量会使程序变得更慢,也会让我们更难以发现数据之间的关系。幸运的是,在真实的世界里,我们通常可以在不损失太多信息的前提下,把一些彼此关联十分紧密的变量进行合并和去掉一些无关变量。
例如,下面左图中的数据点看似在三维空间中但实际上都落在一个二维平面上,这是因为描述数据点的某一个变量与其他两个变量完全线性相关,因此只需要调整一下坐标系就可以去掉一个变量。右图中数字周围的空白部分不提供信息,同样可以去掉。通过对原数据进行一些调整来减少变量个数的方法就叫做降维,而 PCA 是迄今为止使用最广泛降维算法之一。
接下来,我们将介绍 PCA 的基本原理和 Stata 实操的内容。
从上面的介绍中,我们可以看出降维算法主要有两个功能:去相关和去冗余。通俗来说,去相关是指尽可能减少变量间的相关性,去冗余是指去除那些方差很小的无关变量。PCA 算法可以完美地做到上述的功能,让我们来看一下它是怎么实现的。
假设现有数据集
协方差矩阵的求解公式如下:
为达到上述两个目标,我们希望对原数据集进行一个简单变换,使变换后矩阵的协方差矩阵除对角线外元素都等于 0,对角线上元素表示原来变量的方差。很自然我们会想到上文提到的对角化,因为对角矩阵恰好满足上述性质,因此将对原数据集的协方差矩阵进行对角化,得到其特征值和特征向量。
这时我们会发现矩阵
因此,用原数据集
我们将
接下来 我们将进入降维的步骤,去除那些相对不重要的主成分。当去除一些主成分 (坐标轴) 时,我们相当于将数据从高维空间投影到一个低维空间。下图是一个简单的例子,现在将数据从二维平面投影到一条直线上,我们显然会选择投影到
因此在 PCA 算法中,我们只会保留那些投影后保留了原数据最大方差的主成分。而上式告诉我们,特征值越大,主成分方差越大。于是我们按照特征值从大到小的顺序将对应的特征向量从左到右排列,只保留变换矩阵
换一种理解方式,现在
因此 PCA 算法的核心逻辑和特征值的定义式
值得注意的是,虽然我们上面用协方差矩阵进行推导,但在实际中我们常常用相关系数矩阵进行代替,两者作用相似,但前者比起后者更容易受到量纲的干扰。另一方面,并不是所有的数据都适合使用 PCA 算法。PCA 算法对数据变量之间的相关关系有限制,在使用前应先对数据进行检验。
接下来,我们使用书中的一个例子来进行演示 (Jackson,2005)。该数据集记录了一百名成年男性左右耳的听力情况,测量结果是左耳和右耳在四个不同频率上的最小可辨识强度,比如变量 lft1000 指的是左耳在 1000 赫兹时的情况。
. use "https://www.stata-press.com/data/r17/audiometric", clear // 载入数据
. sum lft* rght* // 观察统计性质
. correlate lft* rght* // 相关系数矩阵
. pca lft* rght* // 对所有变量进行 pca 降维
Principal components/correlation Number of obs = 100
Number of comp. = 8
Trace = 8
Rotation: (unrotated = principal) Rho = 1.0000
--------------------------------------------------------------------------
Component | Eigenvalue Difference Proportion Cumulative
-------------+------------------------------------------------------------
Comp1 | 3.92901 2.31068 0.4911 0.4911
Comp2 | 1.61832 .642997 0.2023 0.6934
Comp3 | .975325 .508543 0.1219 0.8153
Comp4 | .466782 .126692 0.0583 0.8737
Comp5 | .34009 .0241988 0.0425 0.9162
Comp6 | .315891 .11578 0.0395 0.9557
Comp7 | .200111 .0456375 0.0250 0.9807
Comp8 | .154474 . 0.0193 1.0000
--------------------------------------------------------------------------
上表是 Stata 给出的解释方差表,第二列给出了每个主成分对应的特征值。该特征值的大小表示了该主成分解释数据集方差 (信息量) 的大小,所有特征值加起来就是原数据集总的方差 (信息量) 大小。该数据集总方差为 8,第一个主成分方差约为 3.93,解释了总方差的 49%(3.93/8)。同理,第二个主成分解释了总方差的 20%(1.6/8)。
每个主成分解释的百分比对应于解释方差表的第四列。第五列第 n 行数据表示前 n 个主成分累计解释方差的占比。解释方差表可以帮助我们判断提取主成分的个数,目前主要有 3 种方法可以帮助大家判断提取主成分的数量,分别是:
. screeplot // 画出陡坡图
. pca lft* rght* // 对所有变量进行 pca 降维
Principal components (eigenvectors)
------------------------------------------------------------------------------------------------------------
Variable | Comp1 Comp2 Comp3 Comp4 Comp5 Comp6 Comp7 Comp8 | Unexplained
-------------+--------------------------------------------------------------------------------+-------------
lft500 | 0.4011 -0.3170 0.1582 -0.3278 0.0231 0.4459 0.3293 -0.5463 | 0
lft1000 | 0.4210 -0.2255 -0.0520 -0.4816 -0.3792 -0.0675 -0.0331 0.6227 | 0
lft2000 | 0.3664 0.2386 -0.4703 -0.2824 0.4392 -0.0638 -0.5255 -0.1863 | 0
lft4000 | 0.2809 0.4742 0.4295 -0.1611 0.3503 -0.4169 0.4269 0.0839 | 0
rght500 | 0.3433 -0.3860 0.2593 0.4876 0.4975 0.1948 -0.1594 0.3425 | 0
rght1000 | 0.4114 -0.2318 -0.0289 0.3723 -0.3513 -0.6136 -0.0837 -0.3614 | 0
rght2000 | 0.3115 0.3171 -0.5629 0.3914 -0.1108 0.2650 0.4778 0.1466 | 0
rght4000 | 0.2542 0.5135 0.4262 0.1591 -0.3960 0.3660 -0.4139 -0.0508 | 0
------------------------------------------------------------------------------------------------------------
上表被称为旋转成分矩阵,其实就是我们上面说到的特征向量矩阵
我们前面提到过,进行 PCA 之前要进行检验,KMO 检验是一个不错的选择,我们来看一下怎么解读:
. estat kmo // 进行 KMO 检验
Kaiser-Meyer-Olkin measure of sampling adequacy
-----------------------
Variable | kmo
-------------+---------
lft500 | 0.7701
lft1000 | 0.7767
lft2000 | 0.7242
lft4000 | 0.6449
rght500 | 0.7562
rght1000 | 0.8168
rght2000 | 0.6673
rght4000 | 0.6214
-------------+---------
Overall | 0.7328
-----------------------
KMO 检验主要用于主成分提取的数据情况。一般来说,KMO 检验系数分布在 0 到 1 之间。如果系数值大于 0.6,则认为样本符合数据结构合理的要求,但大于 0.8 会更好。对于单个变量的分析结果,如果系数大于 0.5,则认为单个变量满足要求;如果系数大于0.8,则认为单个变量结果很好。在本研究中,任一变量的 KMO 检验结果均大于 0.6,总体检验结果大于 0.7,结果不算太好,但仍可以接受。Kaiser (1974) 给出了一个通用的判断 KMO 检验值的法则:
在了解完上面所有信息后,我们可以生成降维后的数据了。我们选择提取前三个主成分,通过原数据生成三个主成分 f_1、f_2、f_3。
. pca lft* rght*, components(3) // 选择前三个主成分降维
. predict f1-f3 // 生成新数据
至此我们已经完成了 PCA 的全部过程。在进行实证研究或机器学习中,经过 PCA 降维后的数据往往会比原数据在模型上拟合精度更高,也更有利于在繁杂的数据中快速找到规律,这其实就是一个数据净化的过程。需要注意的是,PCA 是需要提前对数据进行中心化的,程序里的函数已经自动执行这一步骤,各位读者如果选择自行处理的话要记得加上数据中心化的步骤。
Note:产生如下推文列表的 Stata 命令为:
lianxh 主成分, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh