范文一:主成分分析法
load hald;
indredients=[7 26 6 60;1 29 15 52;11 56 8 20; 11 31 8 47;7 52 6 33];
[pc,score,latent,tsquare]=princomp(ingredients); pc
score
latent
pc =
0.0678 0.6460 -0.5673 0.5062 0.6785 0.0200 0.5440 0.4933 -0.0290 -0.7553 -0.4036 0.5156 -0.7309 0.1085 0.4684 0.4844
score =
-36.8218 6.8709 4.5909 0.3967 -29.6073 -4.6109 2.2476 -0.3958 12.9818 4.2049 -0.9022 -1.1261 -23.7147 6.6341 -1.8547 -0.3786 0.5532 4.4617 6.0874 0.1424 10.8125 3.6466 -0.9130
-0.1350
32.5882 -8.9798 1.6063 0.0818 -22.6064 -10.7259 -3.2365 0.3243 9.2626 -8.9854 0.0169 -0.5437 3.2840 14.1573 -7.0465 0.3405 -9.2200 -12.3861 -3.4283 0.4352 25.5849 2.7817 0.3867 0.4468 26.9032 2.9310 2.4455 0.4116
latent =
517.7969
67.4964
12.4054
0.2372
范文二:R主成分分析法
考虑进出口总额Y与3个自变量国内生产总值X1,存储量X2,总消费量X3(单位为10亿法郎)之间的关系。现收集了11年度的数据,做经典回归分析和主成分回归分析。
运用R做经典回归分析和主成分回归分析 输入数据并作线性回归
> conomy=data.frame(
+ X1=c(149.3,161.1,171.5,175.5,180.8,190.7,202.1,212.4,226.1,231.9,239.0), + X2=c(4.2,4.1,3.1,3.1,1.1,2.2,2.1,5.6,5.0,5.1,0.7),
+ X3=c(108.1,114.8,123.2,126.9,132.1,137.7,146.0,154.1,162.3,164.3,167.6), + Y=c(15.9,16.4,19.0,19.1,18.8,20.4,22.7,26.5,28.1,27.6,26.3)) > lm.sol=lm(Y~X1+X2+X3,data=conomy) > summary(lm.sol)
输出的结果为
Call:
lm(formula = Y ~ X1 + X2 + X3, data = conomy)
Residuals:
Min 1Q Median 3Q Max -0.5297 -0.3875 0.0567 0.2251 0.7834
Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) -10.12620 1.21802 -8.314 7.12e-05 ***
X1 -0.05061 0.07082 -0.715 0.497974 X2 0.58693 0.09481 6.191 0.000449 *** X3 0.28573 0.10303 2.773 0.027558 * ---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4897 on 7 degrees of freedom
Multiple R-squared: 0.9919, Adjusted R-squared: 0.9884 F-statistic: 284.7 on 3 and 7 DF, p-value: 1.124e-07
从计算结果可以看出,按三个变量得到的回归方程为 Y=—10.12799-0.05140X1+0.58695X2+28685X3
但是该方程并不合理,Y是进口量,X1是国内总产值,而对应系数的符号为负,也就是说,国内总产值越高其进口量就越少,这与实际情况不符,究其原因,三个变量之间存在多重共线的问题。
为克服多重共线问题的影响,对变量做主成分分析 > conomy.pr=princomp(~X1+X2+X3,data=conomy,cor=T) > summary(conomy.pr,loadings=TRUE)
输出结果为:
Importance of components:
Comp.1 Comp.2 Comp.3 Standard deviation 1.4139225 0.9990828 0.0515442188 Proportion of Variance 0.6663922 0.3327221 0.0008856022 Cumulative Proportion 0.6663922 0.9991144 1.0000000000
Loadings:
Comp.1 Comp.2 Comp.3 X1 -0.706 0.707 X2 -0.999 X3 -0.707 -0.707
前两个主成分已经达到99%的贡献率,第一主成分是关于国内生产总值和总消费,因此称第一主成分为产销因子,第2主成分只与存储量有关,称为存储因子。又comp3的方差为
0.0515442188^2=0.002690889,近似为0,所以变量存在多重共线的问题。
下面做主成分回归,计算样本的主成分的预测值,并将第一主成分的预测值和第2主成分的预测值放在数据框conomy中,然后对主成分做回归分析。 > pre=predict(conomy.pr)
> conomy$Z1=pre[,1];conomy$Z2=pre[,2] > lm.sol=lm(Y~Z1+Z2,data=conomy) > summary(lm.sol) 运行结果为 Call:
lm(formula = Y ~ Z1 + Z2, data = conomy)
Residuals:
Min 1Q Median 3Q Max -0.89272 -0.26026 0.08328 0.35628 0.66748
Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 21.8909 0.1654 132.368 1.19e-14 *** Z1 -2.9891 0.1170 -25.556 5.89e-09 *** Z2 -0.8295 0.1655 -5.011 0.00104 ** ---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.5485 on 8 degrees of freedom Multiple R-squared: 0.9883, Adjusted R-squared: 0.9854 F-statistic: 339.1 on 2 and 8 DF, p-value: 1.847e-08
回归系数和回归方程均通过检验,而且效果显著,回归方程为 Y=21.8909 -2.9891Z1*-0.8295Z2*
但这只是得到响应变量与与主成分的关系,我们更希望得到响应变量与原变量的关系
> beta=coef(lm.sol);A=loadings(conomy.pr) > X.bar=conomy.pr$center;X.sd=conomy.pr$scale
> coef=(beta[2]*A[,1]+beta[3]*A[,2])/X.sd > beta0=beta[1]-sum(X.bar*coef) > c(beta0,coef)
(Intercept) X1 X2 X3 -9.12564405 0.07275176 0.60944425 0.10626595
即回归方程为
Y=-9.12564405+0.07275176 X1+0.60944425 X2+0.10626595 X3 此时相应的X1 X2 X3的系数均为正数,比原回归方程更合理。
范文三:PCA-主成分分析法
EEW 系统估算地震参数的 PCA 法
野田俊太 , 研究员 山本俊六 , 博士 . 工程师 , 实验室掌门人 佐藤新二 , 高级研究员 地震防灾研究室 /实验室,防灾科技司 /部门
摘要 :新干线 EEW 系统从单站数据快速估计地震参数。这个系统中估计震中位臵,是利用主成分分析法和 B-Δ方法。本文提出了新的方法以提高估算震央的性能 --通过引入可变的时间窗口,替代传统的固定时窗。 结果发现, 新的方法对提高 ? back-azimuth (后方位角) ? 估算的精度和速度分别达到了 28%和 0.25秒。 与以往的方法相比,震中距的估计速度提高了约 1.32秒,相对于传统方法。
关键词:EEW ,早期地震侦测,单台方法, P 波,主要成分分析, B-Δ方法
1. 介绍 .
EEW 系统是一个减少地震灾害的非常有效的方法。 一般情况下, EEW 系统能够提供地震参数, 例如, 在第一个台站侦测到 P 波之后的几秒钟之内,它可以算出震级和震中位臵。中村丰开发的 UrEDAS ,是为 了在地震期间,安全地终止新干线的运行 ---这在大约 20年以前就投入了实际应用。
RTRI 为新干线开发了一个新的 EEW 系统,而且 UrEDAS 被 2004-2005年的新系统所取代。当新系统 中的台站侦测到地震时,震中位臵的确定通过组合以下两种方法(称为 单站方法 ) :
1. 主成分分析(PCA ) (图 1) :通过对初始 P 波的微弱运动的第一主成分的?后 -方位角?的估算;
2.B-Δ法(图 2) :通过拟合函数 Bt*exp(-At )的系数 B 来估算震中距Δ。
在此过程中,第一次震中位臵的确定,此时震级的估算来自于?已被确定?的震中距和?通过使用预 定的地震动烈度的衰减关系?得出的振幅。在单站方法中 , 改进精确度和速度 , 去估算震中位臵是必要的 , 为 了使 new EEW系统更加稳定、有效。
图 1. 使用 PCA 法对“后 -方位角”的估算 . 图 2. 使用 B-Δ法对震中距的估算 .
这项研究首先提出了一个新的方法, 通过引入一个可变的时间窗口来代替传统的固定时窗, 去提高 PCA 的精度和速度。第二,另一个新方法被提议去提高 B-Δ法的速度,通过 使用系数 A 和 B 的时间序列特征 。 2. 改善后方位角的估算 .
本章陈述了一个新的方法,都提高了通过 PCA 法估算的后方位角的精度和速度。
首先, 通过使用 CPA , 估算后方位角的传统程序被描述; 它的精确性被评估, 通过分析 “强震动台站” 、 K-net 观测到的波形纪录,它由国家地球科学研究所和 NIED 运营。
传统方法估算后方位角通过以下步骤:
1.P 波到时的估算通过使用 STA/LTA(短时平均对长时平均)法。
2. 位移数据的计算来自于观测到的加速度记录,通过使用 IIR (无限持久脉冲响应)滤波器。
非常重要的
, 平均长度将比传统的时间窗口的长度更短。
当使用一个可变的时间窗口时 , 误差引起的频率在图 8中显示出来。因此 , 比较图 8(可变时间窗口 ) 与图 5(固定的 0.6秒时间窗口) ,它显示的是, 误差小于 30度的数据的增长 , 然而对于其他的数据是减少的。 误差的均方根 (RMS)为 49.0度,当可变的时间窗口被应用时。这些结果表明, 使用可变时间窗口提高的精 度约为 28%, 或者涉及到与 ?使用固定的 0.6秒时间窗 (在固定时间窗内达到的最高级别的精度) ?相比, 通过常规方法得到的结果是 14%。 其他关于 提高范围 (延伸 ) 的 详细分析,通过?使用参考文件 [8]中被讨论 的可变时间窗口? 。
图 7. 可变时间窗长度的频率分布 图 8. 可变时间窗的情况下,后方位角的估算误差的频率 据认为新推荐的后方位角估算方法对于 EEW 系统是相当有效的。这种方法提供了以下竞争优势 : 1. 除了改善后方位角估算的精度之外 , 与传统的方法相比 , 它也缩短了估算时间,因为可变时窗的平均长度 短于传统时窗 (长度 ) 。
2. 它是如此简单,即使进行实时处理, CPU 负载几乎没有增加。
3. 只有少量代码行被需要,除了常规方法源代码,为了计算?初始 P 波的波长周期?的前一半的长度。 3. 震中距估算的改善 .
本章提出了一个建议:一个新的方法来提高震中距估算的速度,通过 B-Δ法 [6](图 2) 获得。
传统的方法估计震中距通过实施下列步骤:
1. P波到达时间被估算,是通过 STA / LTA方法中,与后方位角估算相同的方式。
2. 通过使用 IIR 滤波器, 为 10-20Hz ,观察到的加速度数据,被执行带通滤波。
3. 震中距 (Δ) , 是从 ?初始 P 波的带通加速度数据? 的 UD 分量包络线的 ?拟合方程 y=Bt*exp(-At)(1) 确定的系数 B ?估算的。
系数 B 与震中距具有良好的相关性, (log 10Δ=-0.498 log10B+1.965) 。参数 B 几乎独立于震级,因为参数 A 与包络线的总体形状有关—主要受到震级的影响;而系数 B 与 P 波的初始部分的斜率有关,主要受震中 距的影响 , 分别地 。 我们称这个常规方法为 ? 2秒 B-Δ法? , 因为第三步中时间窗口的长度, 通常设臵为 2秒。 log 10ΔJMA 和 log 10Δ2s 之间的差异的均方根(RMS ) ,大约是 0.32(ΔJMA 和Δ2s 分别代表着由 JMA 发布的震 中距、其估算通过? 2秒 B-Δ法? ) 。本章中的数据集和前面的章节中的是一样的。
研究的焦点在?为了提高 B-Δ法的速度?的参数 A 和 B 的时间序列上。计算参数 A 和 B ,是在第三 步中以? 0.1至 2.0秒之间?的时间窗为长度, 0.1秒为增量。
因此, A 和 B 的时间序列具有两个重要特性:
1.A 和 B 是高度相关的。
2.A 值随着时间的推移 , 收敛到 各自的 定值。
图 9显示第一特征。 A 和 Log 10B 的时间序列之间的相关系数的平均值和标准偏差,分别为 0.85和 0.15。图 10说明了第二个特点。
图 10(a)和 (b)分别显示 A 的时间序列和它的微分值。 虽然在分析中使用的所有数据被绘制在这些图里,
A 清楚地收敛在 2
秒以内。这些特性反映了?拟合方程 (1)对于初始的 P 波包络线?的 适宜性 。 图 9. 系数 A 和 B 的时间序列的例子 , 显示了 A 、 B 的相关性。
a. 记录结果源于编码为 AKT020的台站 (震中距为 8km), 地震事件 (M5.1),发生于 1999年 2月 26日 14:18,在秋田县。
b. 记录结果源于编码为 WKY006的台站 (震中距为 180km), 地震事件 (M6.9),发生于 2004年 9月 5日 19:07,在纪伊半岛东南。
c. 记录结果源于编码为 NIG013的台站 (震中距为 53km), 地震事件 (M6.8),发生于 1999年 2月 26日 17:56,在中部新泻县。
d. 记录结果源于编码为 MYG006的台站 (震中距为 50km), 地震事件 (M7.2),发生于 2008年 6月 14日 8:43,在岩手—仙台县。
图 10. 系数 A (a )和它的微分值(b )的收敛特性的时间序列。注意 (b)图的纵坐标为 A 值的微分。 这些特点会产生一个可能的新的方法,来提高 B-Δ法的速度。也就是?当 A 的收敛性被确定下来(我们称 这个时间点为 ? (临界 ) 收敛时间? ) 时, 从参数 B ?能估算震中距, 与 2秒 B-Δ法有相同的精确度。 因此 , 该 提案是 -?通过使用两个参数,监控系数 A 的收敛?的一个方法,进行实时处理(图 11) 。第一个参数是的 Tad ,这是 A 的微分值的阈值; 另一个参数是 Dd ,这是在? A 小于 Tad 的‘整段时期内’ ?的微分值的阈 值。 这种新方法被命名为?可变 B-Δ法? 。
图
11.
“可变 B-Δ法”的原理说明图 .
图 12显示了 Dd 和? log 10Δ2s 及 log 10Δvar 之差的均方根 (RMS)?之间的关系 (ΔVAR 意味着通过使用?可 变 B-Δ法?估计的震中距 ) ;图 13显示了 Dd 和平均收敛时间之间的关系。这些图像说明较小的 Tad 和较 长的 Dd , 是分别地(表示收敛的判断更加严重) ,越短的收敛时间和越小的? log 10Δ2 s和 log 10Δvar 之间的 差异的?均方根(RMS )是,分别是一样好。换句话说,这一结果表明,一个交换,在收敛时间和可变 B-Δ法的精度,与 2秒 B-Δ法对照。
图 12. Dd和“ RMS(log10Δ2 s— log10Δvar) ”之间的关系 . 图 13. Dd和平均收敛时间之间的关系 . B-Δ法的速度能从 ‘可变 B-Δ法’, 通过适当地选择参数 (Tad和 Dd) 去提高, 同时保持 ? 2秒 B-Δ法? 的准确性。例如,当 Tad=30(1/s2) , Dd=0.3(s)被选择,收敛时间的均值是约 0.68秒, log 10ΔJMA 和 log 10Δvar 之间的差异的均方根 (RMS)约为 0.30。这个均方根值几乎等于如上所述的,使用? 2秒的 B-Δ法?获 得的值。因此结论是 :对 EEW 系统,可变 B-Δ法是相当有效的,因为速度可以提高约 66%, 同时保持震中 距的估算精度,通过与 2秒 B-Δ法相比较。
4. 结论 .
这项研究提出了两种方法 , 来提高震中位臵的估算性能。在后估计方位角里,准确性和快速性分别提高 了 28%和 23%, 通过引入可变长度的时间窗。对于震中距离估计,速度提高了 66%--同时保持相同水平的 精度,通过使用的系数 A 和 B 的时间序列特征 .
对 EEW 系统,应用这项研究中提出的方法,期望产生更强大、更快捷的地震预警。
致谢
我们想要表达我们衷心的感谢, 向国家地球科学和防灾研究所, 对允许我们使用 K-NET 观察到的波形记录。 参考文献
……
范文四:主成分分析法
主成分分析也称主分量分析,旨在利用降维的思想,把多指标转化为少数几个综
合指标。在实证问题研究中,为了全面、系统地分析问题,我们必须考虑众多影
响因素。这些涉及的因素一般称为指标,在多元统计分析中也称为变量。因为每
个变量都在不同程度上反映了所研究问题的某些信息,并且指标之间彼此有一定
的相关性,因而所得的统计数据反映的信息在一定程度上有重叠。在用统计方法
研究多变量问题时,变量太 多会增加计算量和增加分析问题的复杂性,人们希
望在进行定量分析的过程中,涉及的变量较少,得到的信息量较多。
主成分分析法简介-principal component analysis(PCA) 主成分分析法是
一种数学变换的方法, 它把给定的一组相关变量通过线性变换转成另一组不
相关的变量,这些新的变量按照方差依次递减的顺序排列。在数学变换中保
持变量的总方差不变,使第一变量具有最大的方差,称为第一主成分,第二
变量的方差次大,并且和第一变量不相关,称为第二主成分。依次类推,I
个变量就有I个主成分。
其中Li为p维正交化向量(Li*Li=1),Zi之间互不相关且按照方差由
大到小排列,则称Zi为X的第I个主成分。设X的协方差矩阵为Σ,则Σ
必为半正定对称矩阵,求特征值λi(按从大到小排序)及其特征向量,可
以证明,λi所对应的正交化特征向量,即为第I个主成分Zi所对应的系数
向量Li,而Zi的方差贡献率定义为λi/Σλj,通常要求提取的主成分的数
量k满足Σλk/Σλj>0.85。
编辑本段主成分分析的主要目的
是希望用较少的变量去解释原来资料中的大部分变异,将我们手中许多
相关性很高的变量转化成彼此相互独立或不相关的变量。通常是选出比原始
变量个数少,能解释大部分资料中的变异的几个新变量,即所谓主成分,并
用以解释资料的综合性指标。由此可见,主成分分析实际上是一种降维方法。
编辑本段分析步骤
数据标准化;
求相关系数矩阵;
一系列正交变换,使非对角线上的数置0,加到主对角上;
得特征根xi(即相应那个主成分引起变异的方差),并按照从大到小的顺
序把特征根排列;
求各个特征根对应的特征向量;
用下式计算每个特征根的贡献率Vi;
Vi=xi/(x1+x2+........)
根据特征根及其特征向量解释主成分物理意义。
编辑本段主成分分析法在社会调查中的应用
在社会调查中,对于同一个变量,研究者往往用多个不同的问题来测量
一个人的意见。这些不同的问题构成了所谓的测度项,它们代表一个变量的
不同方面。主成分分析法被用来对这些变量进行降维处理,使它们“浓缩”
为一个变量,称为因子。
在用主成分分析法进行因子求解时,我们最多可以得到与测度项个数一
样多的因子。如果保留所有的因子,就起不到降维的目的了。但是我们知道
因子的大小排列,我们可以对它们进行舍取。那么多小的因子需要舍弃呢?
在一般的行为研究中,我们常常用到的判断方法有两个:特征根大于1法与
碎石坡法。
因为因子中的信息可以用特征根li来表示,所以我们有特征根大于1这
个规则。如果一个因子的特征根大于1就保留,否则抛弃。这个规则,虽然
简单易用,却只是一个经验法则(rule of thumb),没有明确的统计检验。不
幸的是,统计检验的方法在实际中并不比这个经验法则更有效(Gorsuch,
1983)。所以这个经验法则至今仍是最常用的法则。作为一个经验法则,它不
总是正确的。它会高估或者低估实际的因子个数。它的适用范围是20-40个
的测度项,每个理论因子对应3-5个测度项,并且样本量是大的 ( 3100)。
碎石坡法是一种看图方法。如果我们以因子的次序为X轴、以特征根大
小为Y轴,我们可以把特征根随因子的变化画在一个坐标上,因子特征根呈
下降趋势。这个趋势线的头部快速下降,而尾部则变得平坦。从尾部开始逆
向对尾部画一条回归线,远高于回归线的点代表主要的因子,回归线两旁的
点代表次要因子。但是碎石坡法往往高估因子的个数。这种方法相对于第一
种方法更不可靠,所以在实际研究中一般不用。
抛弃小因子、保留大因子之后,降维的目的就达到了。
编辑本段因子旋转
在对社会调查数据进行分析时,除了把相关的问题综合成因子并保留大
的因子,研究者往往还需要对因子与测度项之间的关系进行检验,以确保每
一个主要的因子(主成分)对应于一组意义相关的测度项。为了更清楚的展
现因子与测度项之间的关系,研究者需要进行因子旋转。常见的旋转方法是
VARIMAX旋转。旋转之后,如果一个测度项与对应的因子的相关度很高(>0.5)
就被认为是可以接受的。如果一个测度项与一个不对应的因子的相关度过高
(>0.4),则是不可接受的,这样的测度项可能需要修改或淘汰。
用主成分分析法得到因子,并用因子旋转分析测度项与因子关系的过程
往往被称为探索性因子分析。
在探索性因子分析被接受之后,研究者可以对这些因子之间的关系进行
进一步测试,比用如结构方程分析来做假设检验。
编辑本段主成分分析应用中的问题
1问题的提出主成分分析是一种降维的方法,便于分析问题,在诸多领域
中都有广泛的应用。但有些教科书与论文使用主成分分析时,出现了一些错误
与不足,不能解决实际问题。如一些多元统计分析的教材中,用协方差矩阵的
主成分分析出现了如下错误与不足:①没有明确和判断该数据降维的条件是
否成立。②主成分系数的平方和不为1。③没有明确和判断所用数据是否适
合作单独的主成分分析。④选取的主成分对原始变量没有代表性。以下从相
关性等理论与结果上依次解决上述问题,并给出相应建议。2数据在行为与心
理研究中,常常要求分析某种身份的人的行为特征,如本例中的小学生的日常
行为特征,从而根据这些特征引导小学生向更积极的行为态度发展。这里用文
献[1]的数据见表1,其来自某课题组的调查结果。课题组对北方某小学480
名5~6年级学生的日常行为进行调查,共调查了15项指标如下:S1~对老师
提问的反应、S2~对班级事务的关心、S3~自习课上的表现、S4~对家庭作
业的态度、S5~关心同学的程度、S6~对待劳动的态度、S7~学习上的特殊
兴趣、S8~对待体育锻炼的态度、S9~在娱乐上的偏好、S10~解决问题的思
考方式、S11~对未来的打算
题 名】 主成分分析方法(举例)
【出 处】 /
【日 期】 /
【作 者】 /
【关键词】 /
【正 文】
3. 主成分分析方法应用实例
1) 实例1: 流域系统的主成分分析(张超,1984)
表3.5.1(点击显示该表)给出了某流域系统57个流域盆地的9项变量指标。其中,x1代表流域盆地总高度(m),x2代表流域盆地山口的海拔高度(m),x3代表流域盆地周长(m),x4代表河道总长度(m),x5代表河道总数,x6代表平均分叉率,x7代表河谷最大坡度(度),x8代表河源数, x9代表流域盆地面积(km2)。
注:表中数据详见书本87和88页。
(1) 分析过程:
① 将表3.5.1中的原始数据作标准化处理,然后将它们代入相关系数公式计算,得到相关系数矩阵(表3.5.2)。
② 由相关系数矩阵计算特征值,以及各个主成分的贡献率与累计贡献率(见表3.5.3)。由表3.5.3可知,第一,第二,第三主成分的累计贡献率已高达86.5%,故只需求出第一、第二、第三主成分z1,z2,z3即可。
z3上的载荷
(表3.5.4)。
(2) 结果分析:
▲ 第一主成分z1与x1,x3,x4,x5,x8,x9有较大的正相关,可以看作是流域盆地规模的代表;
▲ 第二主成分z2与x2有较大的正相关,与x7有较大的负相关,分可以看作是流域侵蚀状况的代表;
▲ 第三主成分z3与x6有较大的正相关,可以看作是河系形态的代表;
▲ 根据主成分载荷,该流域系统的9项要素可以被归纳为三类,即流域盆地的规模,流域侵蚀状况和流域河系形态。如果选取其中相关系数绝对值最大者作为代表,则流域面积、流域盆地出口的海拔高度和分叉率可作为这三类要素的代表。
范文五:主成分分析法
一、概述
在处理信息时, 当两个变量之间有一定相关关系时, 可以解释为这两个变量反映此课题 的信息有一定的重叠, 例如, 高校科研状况评价中的立项课题数与项目经费、 经费支出等之 间会存在较高的相关性; 学生综合评价研究中的专业基础课成绩与专业课成绩、 获奖学金次 数等之间也会存在较高的相关性。 而变量之间信息的高度重叠和高度相关会给统计方法的应 用带来许多障碍。
为了解决这些问题, 最简单和最直接的解决方案是削减变量的个数, 但这必然又会导致 信息丢失和信息不完整等问题的产生。 为此, 人们希望探索一种更为有效的解决方法, 它既 能大大减少参与数据建模的变量个数, 同时也不会造成信息的大量丢失。 主成分分析正式这 样一种能够有效降低变量维数,并已得到广泛应用的分析方法。
主成分分析以最少的信息丢失为前提, 将众多的原有变量综合成较少几个综合指标, 通 常综合指标(主成分)有以下几个特点:
↓主成分个数远远少于原有变量的个数 原有变量综合成少数几个因子之后, 因子将可以替代原有变量参与数据建模, 这将大大 减少分析过程中的计算工作量。
↓主成分能够反映原有变量的绝大部分信息 因子并不是原有变量的简单取舍, 而是原有变量重组后的结果, 因此不会造成原有变量 信息的大量丢失,并能够代表原有变量的绝大部分信息。
↓主成分之间应该互不相关
通过主成分分析得出的新的综合指标 (主成分) 之间互不相关, 因子参与数据建模能够 有效地解决变量信息重叠、多重共线性等给分析应用带来的诸多问题。
↓主成分具有命名解释性
总之,主成分分析法是研究如何以最少的信息丢失将众多原有变量浓缩成少数几个因 子,如何使因子具有一定的命名解释性的多元统计分析方法。
二、基本原理
主成分分析是数学上对数据降维的一种方法。其基本思想是设法将原来众多的具有一 定相关性的指标 X1, X2,…, XP (比如 p 个指标) ,重新组合成一组较少个数的互不相关的 综合指标 Fm 来代替原来指标。那么综合指标应该如何去提取,使其既能最大程度的反映原 变量 Xp 所代表的信息,又能保证新指标之间保持相互无关(信息不重叠) 。
设 F1表 示 原 变 量 的 第 一 个 线 性 组 合 所 形 成 的 主 成 分 指 标 , 即
11112121... p p
F a X a X a X =+++, 由数学知识可知,每一个主成分所提取的信息量可用其方
差来度量,其方差 Var(F1)越大,表示 F1包含的信息越多。常常希望第一主成分 F1所含的 信息量最大,因此在所有的线性组合中选取的 F1应该是 X1, X2,…, XP 的所有线性组合中 方差最大的,故称 F1为第一主成分。如果第一主成分不足以代表原来 p 个指标的信息,再 考虑选取第二个主成分指标 F2,为有效地反映原信息, F1已有的信息就不需要再出现在 F2中,即 F2与 F1要保持独立、不相关,用数学语言表达就是其协方差 Cov(F1, F2)=0,所以 F2是与 F1不相关的 X1, X2, …, XP 的所有线性组合中方差最大的, 故称 F2为第二主成分, 依此类推构造出的 F1、 F2、……、 Fm 为原变量指标 X1、 X2…… XP 第一、第二、……、第 m 个主成分。
11111221221122221122... ... ...... ... p p p p
m m m mp p
F a X a X a X F a X a X a X F a X a X a X =+++??=+++??
??=+++? 根据以上分析得知:
(1) Fi与 Fj 互不相关,即 Cov(Fi, Fj) = 0,并有 Var(Fi)=ai’Σai ,其中Σ为 X 的
协方差阵
(2)F1是 X1, X2,…, Xp 的一切线性组合(系数满足上述要求)中方差最大的 , …… , 即 Fm 是与 F1, F2, ……, Fm -1都不相关的 X1, X2, …, XP 的所有线性组合中方差最大者。
F1, F2,…, Fm (m ≤ p )为构造的新变量指标,即原变量指标的第一、第二、……、 第 m 个主成分。
由以上分析可见,主成分分析法的主要任务有两点:
(1)确定各主成分 Fi (i=1, 2,…, m )关于原变量 Xj (j=1, 2 ,…, p)的表达式, 即系数 ij a ( i=1, 2,…, m ; j=1, 2 ,…, p ) 。从数学上可以证明,原变量协方差矩阵 的特征根是主成分的方差, 所以前 m 个较大特征根就代表前 m 个较大的主成分方差值; 原变 量协方差矩阵前 m 个较大的特征值 i λ(这样选取才能保证主成分的方差依次最大)所对应 的特征向量就是相应主成分 Fi 表达式的系数 i a , 为了加以限制, 系数 i a 启用的是 i λ对应的 单位化的特征向量,即有 ' ai ai = 1。
(2)计算主成分载荷,主成分载荷是反映主成分 Fi 与原变量 Xj 之间的相互关联程度:
(, ) (, 1,2,
, ; 1,2, , ) k i ki P Z x i p k m ===
三、主成分分析法的计算步骤 主成分分析的具体步骤如下: (1)计算协方差矩阵
计算样品数据的协方差矩阵:Σ=(sij )p ?p ,其中
1
1()() 1n
ij ki i kj j k s x x n ==---∑ i, j=1, 2,…, p (2)求出Σ的特征值 i λ及相应的正交化单位特征向量 i a
Σ的前 m 个较大的特征值 λ1≥λ2≥… λm>0,就是前 m 个主成分对应的方差, i λ对应的单 位特征向量 i a 就是主成分 Fi 的关于原变量的系数,则原变量的第 i 个主成分 Fi 为:
Fi =' i a X
主成分的方差(信息)贡献率用来反映信息量的大小, i α为:
1
/m
i i i i αλλ==∑
(3)选择主成分
最终要选择几个主成分, 即 F1,F2, …… ,Fm 中 m 的确定是通过方差 (信息) 累计贡献率 G(m)来确定
1
1
() /p
m i k i k G m λλ===∑∑
当累积贡献率大于 85%时, 就认为能足够反映原来变量的信息了, 对应的 m 就是抽取的 前 m 个主成分。
(4)计算主成分载荷
主成分载荷是反映主成分 Fi 与原变量 Xj 之间的相互关联程度,原来变量 Xj (j=1, 2 , …, p) 在诸主成分 Fi (i=1, 2, …, m ) 上的荷载 lij( i=1, 2, …, m ; j=1, 2 , …, p ) 。 :
(, ) (1,2,
, ; 1,2, , ) i j ij l Z X i m j p ===
在 SPSS 软件中主成分分析后的分析结果中, “成分矩阵”反应的就是主成分载荷矩阵。
(5)计算主成分得分
计算样品在 m 个主成分上的得分:
1122... i i i pi p F a X a X a X =+++ i = 1, 2,…, m
实际应用时,指标的量纲往往不同,所以在主成分计算之前应先消除量纲的影响。消 除数据的量纲有很多方法,常用方法是将原始数据标准化,即做如下数据变换:
*1,2,..., ; 1,2,..., ij j
ij j
x x i n j p s -=
==
其中:1
1n j ij i x n ==∑, 2
211() 1n j ij j i s x n ==--∑ 根据数学公式知道,①任何随机变量对其作标准化变换后,其协方差与其相关系数是 一回事, 即标准化后的变量协方差矩阵就是其相关系数矩阵。 ②另一方面, 根据协方差的公 式可以推得标准化后的协方差就是原变量的相关系数, 亦即, 标准化后的变量的协方差矩阵 就是原变量的相关系数矩阵。也就是说,在标准化前后变量的相关系数矩阵不变化。
根据以上论述,为消除量纲的影响,将变量标准化后再计算其协方差矩阵,就是直接 计算原变量的相关系数矩阵,所以主成分分析的实际常用计算步骤是: ☆计算相关系数矩阵
☆求出相关系数矩阵的特征值 i λ及相应的正交化单位特征向量 i a ☆选择主成分 ☆计算主成分得分
总结:原指标相关系数矩阵相应的特征值 λi 为主成分方差的贡献,方差的贡献率为
1
/p
i i i i αλλ==∑, i α越大,说明相应的主成分反映综合信息的能力越强,可根据 λi 的大小
来提取主成分。每一个主成分的组合系数(原变量在该主成分上的载荷) i a 就是相应特征 值 λi 所对应的单位特征向量。 二、主成分分析的计算步骤 1、计算相关系数矩阵
r ij (i , j =1, 2,…, p )为原变量 x i 与 x j 的相关系数, r ij =r ji ,其计算公式为
2
解特征方程 ,常用雅可比法(Jacobi )求出特征值,并使其按大小顺序排
列 ;
分别求出对应于特征值 的特征向量 ,要求 =1 其中
表示向量 的第 j 个分量。 3、计算主成分贡献率及累计贡献率
贡献率:
累计贡献率:
一般取累计贡献率达 85%-95%的特征值, 所对应的第 1、第 2、…、第 m
(m ≤ p )个主成分。 4、计算主成分载荷
?
?
??
???
????
??
?=pp p p p p r r r r r r r r r R 2
122221
11211∑∑∑===----=
n
k n
k j kj
i ki
n
k j kj i ki
ij x
x
x x
r 11
2
2
1
) ()
()
)((0=-R I λ021≥≥≥≥p λλλ i λ) , , 2, 1(p i e i L =i
e 1
1
2
=∑=p
j ij
e
ij e i e )
, , 2, 1(1
p i p
k k
i
L =∑=λ
λ)
, , 2, 1(11
p i p
k k
i
k k
L =∑∑==λ
λm λλλ, , , 21L )
, , 2, 1, () , (p j i e x z p l ij i j i ij L ===λ
5、各主成分得分
三、主成分分析法在 SPSS 中的操作 1、指标数据选取、收集与录入(表 1)
2、 Analyze → Data Reduction → Factor Analysis,弹出 Factor Analysis 对话框:
3、把指标数据选入 Variables 框, Descriptives: Correlation Matrix 框组中选中 Coefficients, 然后点击 Continue, 返回 Factor Analysis 对话框,单击 OK 。
注意:SPSS 在调用 Factor Analyze 过程进行分析时 , SPSS 会自动对原始数据进行标 准化处理 , 所以在得到计算结果后的变量都是指经过标准化处理后的变量 , 但 SPSS 并不直
?
??
??
???????=nm n n m m z z z z z z z z z Z 2
1
22221
11211
接给出标准化后的数据 , 如需要得到标准化数据 , 则需调用 Descriptives 过程进行计算。
从表 3 可知 GDP 与工业增加值 , 第三产业增加值、 固定资产投资、 基本建设投资、 社会 消费品零售总额、 地方财政收入这几个指标存在着极其显著的关系 , 与海关出口总额存在着 显著关系。可见许多变量之间直接的相关性比较强 , 证明他们存在信息上的重叠。
主成分个数提取原则为主成分对应的特征值大于 1的前 m 个主成分。 特征值在某种程度上 可以被看成是表示主成分影响力度大小的指标 , 如果特征值小于 1, 说明该主成分的解释力 度还不如直接引入一个原变量的平均解释力度大 , 因此一般可以用特征值大于 1作为纳入标 准。 通过表 4( 方差分解主成分提取分析 ) 可知 , 提取 2个主成分 , 即 m=2, 从表 5( 初始因子 载荷矩阵 ) 可知 GDP 、工业增加值、第三产业增加值、固定资产投资、基本建设投资、社会 消费品零售总额、 海关出口总额、 地方财政收入在第一主成分上有较高载荷 , 说明第一主成 分基本反映了这些指标的信息 ; 人均 GDP 和农业增加值指标在第二主成分上有较高载荷 ,
说明第二主成分基本反映了人均 GDP 和农业增加值两个指标的信息。 所以提取两个主成分是 可以基本反映全部指标的信息 , 所以决定用两个新变量来代替原来的十个变量。 但这两个新 变量的表达还不能从输出窗口中直接得到 , 因为“ Component Matrix ”是指初始因子载荷矩 阵 , 每一个载荷量表示主成分与对应变量的相关系数。
用表 5( 主成分载荷矩阵 ) 中的数据除以主成分相对应的特征值开平方根便得到两个主 成分中每个指标所对应的系数。 将初始因子载荷矩阵中的两列数据输入 ( 可用复制粘贴的方 法 ) 到数据编辑窗口 ( 为变量 B1、 B2) , 然后利用“ Transform → Compute Variable” , 在 Compute Variable对话框中输入“ A1=B1/SQR(7.22)” [注 : 第二主成分 SQR 后的括号中填 1.235, 即可得到特征向量 A 1(见表 6) 。同理 , 可得到特征向量 A 2。将得到的特征向量与标准 化后的数据相乘 , 然后就可以得出主成分表达式 [注 : 因本例只是为了说明如何在 SPSS 进 行主成分分析 , 故在此不对提取的主成分进行命名 , 有兴趣的读者可自行命名。
标准化:通过 Analyze → Descriptive Statistics → Descriptives 对话框来实现 : 弹出 Descriptives 对话框后 , 把 X 1~X 10选入 Variables 框 , 在 Save standardized values as
variables 前的方框打上钩 , 点击“ OK ” , 经标准化的数据会自动填入数据窗口中 , 并以 Z 开头命名。
以每个主成分所对应的特征值占所提取主成分总的特征值之和的比例作为权重计算主 成分综合模型 , 即用第一主成分 F1 中每个指标所对应的系数乘上第一主成分 F1 所对应的 贡献率再除以所提取两个主成分的两个贡献率之和 , 然后加上第二主成分 F2 中每个指标所 对应的系数乘上第二主成分 F2 所对应的贡献率再除以所提取两个主成分的两个贡献率之和 , 即可得到综合得分模型
: