第5章主成分分析与经验正交分解
5.1 主分量分析的数学模型
当存在若干个随机变量时,寻求它们的少量线性组合(即主成分),用以解释这些随机
变量,是很必要的。首先我们看一个例子。
例5.1 为了调查学生的身材状况,可以测量他们的身高(
)、体重(
)、胸围(
)和坐高(
)。可是用这4个指标表达学生身材状况不方便。但若用
=3.6356
+3.3242
+2.4770
+2.1650
表示学生身体魁梧程度;用
=-3.9739
+1.3582
+3.7323
-1.5729
表示学生胖瘦程度。则这两个指标(
,
)很好概括了4个指标(
-
)。
例中,学生不同,身高(
)、体重(
)、胸围(
)和坐高(
)不同;(
,
,
,
)是4维随机向量;
,
是他们的2个线性组合,
,
能很好表示
,
,
,
的特性。类似的问题在许多地方出现:可观测的随机变量很多,需要选出所有所有随机变量的少数线性组合,使之尽可能刻划全部随机变量的特性,选出的线性组合就是诸多随机变量的主成分,又称为主分量。寻求随机向量主成分,并加以解释,称为主成分分析,又称为主分量分析。主成分分析在许多学科中都有应用,细节可参看张尧廷(1991)、Richard(2003),主成分分析在气象等科学中称为PCA方法,见吴洪宝(2005)。
主成分分析的数学模型是:对于随机向量X,想选一些常数向量
,用
尽可能多反映随机向量X的主要信息,也即
尽量大。但是
的模可以无限增大,从而使
无限变大,这是我们不希望的;于是限定
模的大小,而改变
各分量的比例,使
最大;通常取
的模为1最方便。
定义5.1 设随机向量
二阶矩存在,若常数向量
,在条件
=1下
使
最大,则称
是X的第一主成分或第一主分量。
由定义可见,
尽可能多地反映原来p个随机变量变化的信息。但是一个主成分往往不能完全反映随机向量特色,必须建立其它主成分,它们也应当最能反映随机向量变化,而且他们应当与第一主成分不相关(不包含
的信息)。
定义5.2 若常数向量c=
在条件
=l,
下,使
最大,
则称
是 X的第二主成分;若常数向量c=
在条件
=l,
,
下,使
最大,则称
是 X的第三主成分;…。
当随机向量方差已知时,定理5.1给出主成分的计算公式。
定理5.1 设随机向量
方差存在为
。
特征值从大到小为
,
对应的彼此正交单位特征向量为
。则X的第j个主成分
为
与X的内积,即
(5.1)
且
证明:任取p维单位向量c,必有
。于是
,而在条件
下,当
,
即
时,
最大,所以X的第一主成分是
与X的内积
。由条件
,可得
,于是
,从而
;
所以在条件
=1、
下,当
时,
最大,所以X的第2个主成分为
与X的内积
。对第三,第四……主成分同样可证。
由证明过程可见:
。它称为第i个主成分的方差贡献,表示第i个主成分变化大小,从而反映第i个主成分提供的信息的大小。
例5.2 设
,且
则
=3.87939,
=[0.293128,-0.84403,-0.449099]
=1.6527,
=[0.449099,-0.293128,0.84403]
=0.467911,
=[0.84403,0.449099,-0.293128]
所以第一主成分就是
=0.293128
-0.84403
-0.449099
;
第二主成分就是
=0.449099
-0.293128
+0.84403
;
第三主成分就是
=0.84403
+0.449099
-0.293128
。
它们的方差贡献分别是
;
;
。
定义5.3
称为主成分
的方差贡献率;
称为前k个主成分的累计方差贡献率;
与X第k个分量的相关系数
称为因子负荷量。
当某个主成分的方差贡献率很小时,认为它提供的信息很少,可以略去此主成分。通常取q,使前q个主成分的累计方差贡献率达到70%-80%,然后只考虑前q个主分量,用它们解释随机向量X的特性,其余主成分认为是观测误差等随机因素造成的。
在实际问题中,X的每一分量可取不同单位,单位取小时(例如长度单位取毫米,甚至微米)该分量的方差会变大,从而在主成分中变得突出;而单位选取不应影响主成分。为了避免量纲对主成分的影响。常常将随机变量都标化,即令
,它就是无量纲量,令
再求X*的主成分,即
标准
excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载
化后的主成分。将
代入,可求随机向量X的主成分。容易证明
定理5.2 设随机向量X的相关阵为
,
特征值为
,
对应的彼此正交单位特征向量为
,则标准化后X的第j个主成分是
。
因此,标准化后的主成分称为由相关阵决定的主成分。直接由随机向量的协方差阵算出的主成分称为由协差阵决定的主成分。
同样一组随机变量,用它们的协差阵和相关阵求出的主成分是不一样的。这是因为优化的
准则
租赁准则应用指南下载租赁准则应用指南下载租赁准则应用指南下载租赁准则应用指南下载租赁准则应用指南下载
(目标函数)不同:前者
要求
对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗
=
最大,而后者要求
==
最大,其中
。
例5.3 (协差阵和相关阵决定的主成分不同)设随机变量
;其协方差阵是
,特征值和特征向量是
,
。因而由协方差阵决定的主成分是:
,
。
但随机变量
标准化后得到
;其中
。
X*的协差阵即X的相关阵是
,其特征值和特征向量是
,
从而由相关阵决定的主成分是:
。
由于主成分由方差决定,可以略去常数,因而由相关阵得到的主成分可写为:
,
可见由协方差阵与相关阵决定的主成分不同。
5.2 样本主成分及其计算
5.2.1 样本主成分
实际问题中随机向量的协差阵、相关阵都是未知的,只能得到样品
。这时总用样本协差阵与样本相关阵代替协差阵、相关阵求主成分。
定义5.4 样本协差阵与样本相关阵的特征向量,计算主成分。所得的主成分称为样本主成分。
这样求主成分是有道理的:若总体
,
的特征值和正交单位特征向量是
和
;
是
的极大似然估计,即
。
的特征值为
,
相应正交单位特征向量为
,则可证
定理5.3 若X服从正态分布,则
是
的极大似然估计;
是
的极大似然估计。
因此,若X服从正态分布,应当用第j个样本主成分
作为总体主成分
的估计值。从样本协差阵或样本相关阵出发,做主成分分析,所得样本主成分通常简称为主成分。
通常取
为样本协差阵(
的无偏估计),由
或R算出的样本相关阵是相同的,所产生(相关差阵决定)的主成分当然相同。而R与
有相同的特征向量,R的特征值是
特征值的n/(n-1)倍。因而由R与
所产生的(协方差阵决定的)主成分相同。
若X不一定服从正态分布,这时仍可由样本协差阵R或相关阵
出发,计算主成分。
同上节指出的一样:样本相关阵和样本协差阵决定的主成分是不同的。
5.2.2 SAS软件计算样本主成分
样本主成分的计算量很大,通常用软件计算,以下介绍用SAS软件计算的基本方法。
SAS调用PRINCOMP过程(即主成分过程)作主成分分析。PROC PRINCOMP 过程对输入资料文件执行主成分分析。其输入资料文件可以是原始数据,也可以是一个相关系数矩阵,或是协方差阵。输出资料则包括特征根、特征向量及标准化的主成分值。
主成分分析是一个多变量统计程序,可用来鉴定多个数值变量之间的关系。主成分分析除了用来概述变量之间的关系外,还可用来削减回归或集群分析中变量的数目。它的主要目的是求出一组变量的线性组合(即主成分),这些线性组合就是原变量矩阵的特征向量。每个向量的内积就是该向量对原变量群能解释的方差百分比。这些特征向量之间应该是彼此线性独立的。
PROC PRINCOMP语法
PROC PRINCOMP DATA= SAS-data-set /*输入资料文件名称*/
OUT= SAS-data-set /*输出资料文件名称*/
OUTSTAT= SAS-data-set /*输出资料文件名称*/
NOINT
COVARIANCE(COV)
N= n
STANDARD(STD)
PREFIX= name
NOPRINT
SINGULAR= value
VARDEF= DF|N|WEIGHT|WDF; 或N,或WGT,或WDF)
VAR variable-list; /*指明那些数值变量作主成分分析*/
PARTIAL variable-list;
FREQ variable;
WEIGHT variable;
BY variable-list;
调用PRLNCOMP过程时常用两个语句:即PROC PRINCOMP ,VAR。
(1) PROC PRINCOMP语句。
一般形式是 PROC PRINCOMP;其功能是调用PRINCOMP过程。加选项cov指示电脑用协差阵计算样本主成分,不加选项cov则电脑用相关阵计算主成分;加选项out=文件名,指示电脑将每个观测的主成分得分存入一个数据集,即“文件名”所表示的数据集,加选项n=k指示电脑只计算k个主成分,不加选项n=k则电脑计算全部p个主成分。例如proc princomp data=wang1 out=wang2 n=3;指示电脑对数据集wang1中数据做主成分分析,求3个主成分,并将各次观测的主成分得分存入数据集wang2。
(2) VAR语句
其功能是规定要分析的变量。例如var x1-x3 u1 v2;表示将变量x1,x2,x3,u1,v作为随机向量进行主成分分析。
计算主成分固然重要,解释主成分的意义更重要。下面我们介绍用SAS作主成分分析的实例,并对于算出的主成分加以解释,希望学者对练习题中的主成分也试作解释。
例5.4 北京1951~1976年冬季的气温资料如表5-1,第一列为年度,第二列为该年12月的月平均温度。第三、四列为次年1、2月的月平均温度。试做主成分分析。
表 5-1 北京1951~1976年冬季月平均气温
year
x1
x2
x3
1951
1.0
-2.7
-4.3
1952
-5.3
-5.9
-3.5
1953
-2.0
-3.4
-0.8
1954
-5.7
-4.7
-1.1
1955
-0.9
-3.8
-3.1
1956
-5.7
-5.3
-5.9
1957
-2.1
-5.0
-1.6
1958
0.6
-4.3
-0.2
1959
-1.7
-5.7
2.0
1960
-3.6
-3.6
1.3
1961
-3.0
-3.1
-0.8
1962
0.1
-3.9
-1.1
1963
-2.6
-3.0
-5.2
1964
-1.4
-4.9
-1.7
1965
-3.9
-5.7
-2.5
1966
-4.7
-4.8
-3.3
1967
-6.0
-5.6
-4.9
1968
-1.7
-6.4
-5.1
1969
-3.4
-5.6
-2.0
1970
-3.1
-4.2
-2.9
1971
-3.8
-4.9
-3.9
1972
-2.0
-4.1
-2.4
1973
-1.7
-4.2
-2.0
1974
-3.6
-3.3
-2.0
1975
-2.7
-3.7
0.1
1976
-2.4
-7.6
-2.2
解:因为所有变量单位相同,可用协方差阵求主成分。以变量year Dec Jan Feb分别表示年度、12月、1月、2月的温度。采用下列程序
data temperat; /*建立数据集temperat*/
input year Dec Jan Feb; /*建立变量year、Dec、Jan和Feb*/
cards; /*以下为数据体*/
1951 1.0 -2.7 -4.3
1952 -5.3 -5.9 -3.5
1953 -2.0 -3.4 -0.8
1954 -5.7 -4.7 -1.1
1955 -0.9 -3.8 -3.1
1956 -5.7 -5.3 -5.9
1957 -2.1 -5.0 -1.6
1958 0.6 -4.3 -0.2
1959 -1.7 -5.7 2.0
1960 -3.6 -3.6 1.3
1961 -3.0 -3.1 -0.8
1962 0.1 -3.9 -1.1
1963 -2.6 -3.0 -5.2
1964 -1.4 -4.9 -1.7
1965 -3.9 -5.7 -2.5
1966 -4.7 -4.8 -3.3
1967 -6.0 -5.6 -4.9
1968 -1.7 -6.4 -5.1
1969 -3.4 -5.6 -2.0
1970 -3.1 -4.2 -2.9
1971 -3.8 -4.9 -3.9
1972 -2.0 -4.1 -2.4
1973 -1.7 -4.2 -2.0
1974 -3.6 -3.3 -2.0
1975 -2.7 -3.7 0.1
1976 -2.4 -7.6 -2.2
; /*空语句,结束数据体*/
proc princomp cov; /* 用协差阵做主成分分析*/
var Dec Jan Feb; /* 对变量Dec Jan Feb 作主成分分析*/
run;
执行上述程序,得到得许多表,主要的是:基本统计量(Simple Statistic);协方差矩阵(Covariance Matrix);样本协差阵的特阵值表(Eigenvalues of the Covariance Matrix)、方差贡献、方差贡献率及累计方差贡献率;样本协差阵的特征向量表(即主成分的系数表,Eigenvectors)。这些表及分析如下
Eigenvalues
Eigenvalue Difference Proportion Cumulative
PRIN1 4.79742 2.06927 0.552919 0.55292
PRIN2 2.72815 1.57720 0.314429 0.86735
PRIN3 1.15095 . 0.132652 1.00000
上表是样本协差阵的特征值表(表头为Eigenvalues),其中PRIN1、PRIN2、PRIN3表示3个主成分,上表第2列给出样本协差阵的特征值,第4列给出方差贡献,第5列给出方差贡献累计百分比。由于前两个特阵值方差贡献累计百分比等于0.867354,它大于0.7,所以只需取两个主成分。
Eigenvectors
PRIN1 PRIN2 PRIN3
DEC 0.643587 0.709882 -.286116
JAN 0.213039 0.192899 0.957812
FEB 0.735126 -.677390 -.027085
上表是特征向量表(表头为Eigenvectors)上表给出所考察变量样本协差阵的特征向量(0.643587,0.213039,0.735126)’、(0.709882,0.192899,-0.677390)’和(-0.286116,0.957812,-0.027085)’。因此第一、二、三主成分分别是
=0.643587Dec+0.213039Jan+0.735126Feb,
=0.709882Dec+0.192899Jan-0.677390Feb,
=-0.286116Dec+0.957812Jan-0.027085Feb
由于第一主成分中Dec,Feb系数是较大正数,Jan系数是较小正数,说明第一主成分主要表示冬季气温偏高的程度,由于1月分的系数变化较小,冬季气温偏高主要由12月,2月温度的偏高形成。第二主成分Dec系数与Feb系数反号较大,反映第二主成分主要表示12月与2月温度距平的反差,即12月温度距平减去2月温度距平所得值的反差。
例5.5 美国各州犯罪率情况如表5-2。试以murder(谋杀),rape(强奸),robbery(抢劫),assult(斗殴),burglary(夜盗),larceny(偷窃),auto(汽车犯罪)为7元随机向量,做主成分分析。
表 5-2 美国各州犯罪率(十万人中犯罪人数)
murder
rape
robbery
assult
burglary
larceny
auto
Albama
14.2
25.2
96.8
278.3
1135.5
1881.9
280.7
Alaska
10.8
51.6
96.8
284.0
1331.7
3369.8
753.3
Arirona
9.5
34.2
138.2
312.3
2346.1
4467.4
439.5
Arkansas
8.8
34.2
138.2
312.3
2346.1
4467.4
439.5
Califonia
11.5
49.4
287.0
358.0
2139.4
3499.8
663.5
Colorado
6.3
42.0
170.7
292.9
1935.2
3903.2
477.1
Conecticat
4.2
16.8
129.5
131.8
1346.0
2620.7
593.2
Delaware
6.0
24.9
157.0
194.2
1682.6
3678.4
467.0
Florida
10.2
39.6
187.9
449.1
1859.9
3840.5
351.4
Geogia
11.7
31.1
140.5
256.5
1351.1
2170.2
297.9
Hawaii
7.2
25.5
128.0
64.1
1911.5
3920.4
489.4
Idaho
5.5
19.4
39.6
172.5
1050.8
2599.6
237.6
Illinois
9.9
21.8
211.3
209.0
1085.0
2828.5
528.6
Indiana
7.4
26.5
123.2
153.5
1086.2
2498.7
377.4
Iowa
2.3
10.6
41.2
89.8
812.5
2685.1
219.9
Kansas
6.6
22.0
100.7
180.5
1270.4
2739.3
244.3
Kentaky
10.1
19.1
81.1
123.3
872.2
1662.1
245.4
Loisana
15.5
30.9
142.9
335.5
1165.5
2469.9
337.7
Maine
2.4
13.5
38.7
170.0
1253.1
2350.7
246.9
Maryland
8.0
34.8
292.1
358.9
1400.0
3177.7
428.5
Masschusetts
3.1
20.8
169.1
231.6
1532.2
2311.3
1140.1
Michigan
9.3
38.9
261.9
274.6
1522.7
3159.0
545.5
Minnesota
2.7
19.5
85.9
85.8
1134.7
2559.3
343.1
Mississippi
14.3
19.6
65.7
189.1
915.6
1239.9
144.4
Missouri
9.6
28.3
189.0
233.5
1318.3
2424.2
378.4
Montana
5.4
16.7
39.2
156.8
804.9
2773.2
309.3
Nebraska
3.9
18.1
64.7
112.7
760.0
2316.1
249.1
Nevada
15.8
49.1
323.1
355.0
2453.1
4212.6
559.2
Mew Hampashare
3.2
10.7
23.2
76.0
1041.7
2343.9
293.4
New Jersey
5.6
21.0
180.4
185.1
1435.8
2774.5
511.5
New Maxico
8.8
39.1
109.6
343.4
1418.7
3008.6
259.5
New York
10.7
29.4
472.6
319.1
1728.0
2782.0
745.8
North Carolina
10.6
17.0
61.3
318.3
1154.1
2037.8
192.1
North Dakoda
100.9
9.0
13.3
43.8
446.1
1843.0
144.7
Ohio
7.8
27.3
190.5
181.1
1216.0
2696.8
400.4
Oklahoma
8.6
29.2
73.8
205.0
1288.2
2228.1
326.8
Oregan
4.9
39.9
124.1
286.9
1636.4
3506.1
388.9
Pennsyvania
5.6
19.0
130.3
128.0
877.5
1624.1
333.2
Rhode Island
3.6
10.5
86.5
201.0
1849.5
2844.1
791.4
South Carolina
11.9
33.0
105.9
485.3
1613.6
2342.4
245.1
South Dakoda
2.0
13.5
17.9
155.7
570.5
1704.4
147.5
Tennessee
10.1
29.7
145.8
203.9
1259.7
1776.5
314.0
Texas
13.3
33.8
152.4
208.2
1603.1
2988.7
397.6
Utah
3.5
20.3
68.8
147.3
1171.6
3004.6
334.5
Vermont
1.4
15.9
30.8
101.2
1348.2
2201.0
265.2
Virginia
9.0
23.3
92.1
165.7
986.2
2521.2
226.7
Wasinton
4.3
39.6
106.2
224.8
1605.6
3386.9
360.3
West Viginia
6.0
13.2
42.2
90.9
597.4
1341.7
163.3
Wiskonsin
2.8
12.9
52.2
63.7
846.9
2614.2
220.7
Wyoming
5.4
21.9
39.7
173.9
811.6
2772.2
282.0
解:评估美国各州犯罪率时,用7种犯罪率为7维随机向量,以50个州的统计数据为50次观测。考虑不同犯罪的犯罪率差异很大,用相关阵计算主成分。采用程序
data crime; /*建立数据集crime*/
input state $ 1-15 murder rape robbery assult burglary larceny auto;
/*建立变量state murder rape robbery assult burglary larceny auto。state $ 1-15表示前15列存州名。murder rape robbery assult burglary larceny auto 表7种罪的犯罪率*/
cards; /*以下为数据体*/
Albama 14.2 25.2 96.8 278.3 1135.5 1881.9 280.7
Alaska 10.8 51.6 96.8 284.0 1331.7 3369.8 753.3
Arirona 9.5 34.2 138.2 312.3 2346.1 4467.4 439.5
Arkansas 8.8 34.2 138.2 312.3 2346.1 4467.4 439.5
Califonia 11.5 49.4 287.0 358.0 2139.4 3499.8 663.5
Colorado 6.3 42.0 170.7 292.9 1935.2 3903.2 477.1
Conecticat 4.2 16.8 129.5 131.8 1346.0 2620.7 593.2
Delaware 6.0 24.9 157.0 194.2 1682.6 3678.4 467.0
Florida 10.2 39.6 187.9 449.1 1859.9 3840.5 351.4
Geogia 11.7 31.1 140.5 256.5 1351.1 2170.2 297.9
Hawaii 7.2 25.5 128.0 64.1 1911.5 3920.4 489.4
Idaho 5.5 19.4 39.6 172.5 1050.8 2599.6 237.6
Illinois 9.9 21.8 211.3 209.0 1085.0 2828.5 528.6
Indiana 7.4 26.5 123.2 153.5 1086.2 2498.7 377.4
Iowa 2.3 10.6 41.2 89.8 812.5 2685.1 219.9
Kansas 6.6 22.0 100.7 180.5 1270.4 2739.3 244.3
Kentaky 10.1 19.1 81.1 123.3 872.2 1662.1 245.4
Loisana 15.5 30.9 142.9 335.5 1165.5 2469.9 337.7
Maine 2.4 13.5 38.7 170.0 1253.1 2350.7 246.9
Maryland 8.0 34.8 292.1 358.9 1400.0 3177.7 428.5
Masschusetts 3.1 20.8 169.1 231.6 1532.2 2311.3 1140.1
Michigan 9.3 38.9 261.9 274.6 1522.7 3159.0 545.5
Minnesota 2.7 19.5 85.9 85.8 1134.7 2559.3 343.1
Mississippi 14.3 19.6 65.7 189.1 915.6 1239.9 144.4
Missouri 9.6 28.3 189.0 233.5 1318.3 2424.2 378.4
Montana 5.4 16.7 39.2 156.8 804.9 2773.2 309.3
Nebraska 3.9 18.1 64.7 112.7 760.0 2316.1 249.1
Nevada 15.8 49.1 323.1 355.0 2453.1 4212.6 559.2
Mew Hampashare 3.2 10.7 23.2 76.0 1041.7 2343.9 293.4
New Jersey 5.6 21.0 180.4 185.1 1435.8 2774.5 511.5
New Maxico 8.8 39.1 109.6 343.4 1418.7 3008.6 259.5
New York 10.7 29.4 472.6 319.1 1728.0 2782.0 745.8
North Carolina 10.6 17.0 61.3 318.3 1154.1 2037.8 192.1
North Dakoda 100.9 9.0 13.3 43.8 446.1 1843.0 144.7
Ohio 7.8 27.3 190.5 181.1 1216.0 2696.8 400.4
Oklahoma 8.6 29.2 73.8 205.0 1288.2 2228.1 326.8
Oregan 4.9 39.9 124.1 286.9 1636.4 3506.1 388.9
Pennsyvania 5.6 19.0 130.3 128.0 877.5 1624.1 333.2
Rhode Island 3.6 10.5 86.5 201.0 1849.5 2844.1 791.4
South Carolina 11.9 33.0 105.9 485.3 1613.6 2342.4 245.1
South Dakoda 2.0 13.5 17.9 155.7 570.5 1704.4 147.5
Tennessee 10.1 29.7 145.8 203.9 1259.7 1776.5 314.0
Texas 13.3 33.8 152.4 208.2 1603.1 2988.7 397.6
Utah 3.5 20.3 68.8 147.3 1171.6 3004.6 334.5
Vermont 1.4 15.9 30.8 101.2 1348.2 2201.0 265.2
Virginia 9.0 23.3 92.1 165.7 986.2 2521.2 226.7
Wasinton 4.3 39.6 106.2 224.8 1605.6 3386.9 360.3
West Viginia 6.0 13.2 42.2 90.9 597.4 1341.7 163.3
Wiskonsin 2.8 12.9 52.2 63.7 846.9 2614.2 220.7
Wyoming 5.4 21.9 39.7 173.9 811.6 2772.2 282.0
;
proc princomp out=crimprin; /*调用PRINCOMP过程,用相关阵做主成分分析*/
var murder rape robbery assult burglary larceny auto; /*对这7个变量做分析*/
run;
执行以上程序,电脑按相关阵做主成分分析;输出主要数表有:样本相关阵的特征值(表头为Eigenvalues of the Correlation Matrix)表,方差贡献、方差贡献率及累计方差贡献率;样本相关阵的特征向量(表头为Eigenvectors)。表及解释如下
Eigenvalues of the Correlation Matrix
Eigenvalue Difference Proportion Cumulative
1 3.81730007 2.78454963 0.5453 0.5453
2 1.03275044 0.22145080 0.1475 0.6929
3 0.81129963 0.14770303 0.1159 0.8088
4 0.66359660 0.35782066 0.0948 0.9036
5 0.30577594 0.06348335 0.0437 0.9472
6 0.24229259 0.11530785 0.0346 0.9819
7 0.12698474 0.0181 1.0000
Eigenvectors
Prin1 Prin2 Prin3 Prin4 Prin5 Prin6 Prin7
murder -.094836 0.893895 0.335604 0.264209 0.087862 0.037372 -.020129
rape 0.433768 0.218170 -.298382 -.102754 -.033667 -.772201 -.259286
robbery 0.398823 0.091935 0.367321 -.422729 -.696268 0.173693 0.062497
assult 0.39223 0.2585 -.37199 -.431946 0.445511 0.353752 0.361585
burglary 0.463531 -.067937 -.044742 0.305199 0.096792 0.445645 -.690946
larceny 0.402967 -.071041 -.14078 0.678773 -.216768 -.005623 0.55226
auto 0.335705 -.261558 0.709373 -.021882 0.501519 -.219929 0.123736
由特征值表(表头为 Eigenvalues of the Correlation Matrix),第5列可见,前3个特征值所占比例之和为0.80,只要取3个主成分就够了。由特征向量表(表头为Eigenvectors),从第2列起,每列是1个特征向量。第1个特征向量各个分量值大体相同,近似于
=0.38;所以第1主成分表示各州犯罪程度的严重性。第2个特征向量各分量对应murder,rape, assult,分量值为负的,对应burglary,larceny,auto分量是正的,murder,rape, assult暴力程度重, burglary,larceny,auto暴力程度轻,因此第二主成分反映暴力程度的轻重,第二主成分的值越大,暴力成分越轻。第三主成分的特性不明显,不考虑。第一、第二主成分分别为:
y1=0.303311murder+0.432675rape+0.391443robbery+0.401331assult
+0.4434023burglary+0.361074larceny+0.29296226auto;
y2=-0.6634076murder-0.167388rape+0.019456robbery-0.335621assult
+0.237752burglary+0.391665 arceny+0.496972 auto
许多统计资料简化成样本协差阵,或样本相关阵;这时仍可用SAS的princomp过程计算,只是在data步输入数据时要用“_type_=COV”语句说明。
例5.6 测量雄龟甲的长、宽、厚,并求其自然对数,得到变量
;所得24只龟数据的协方差阵如下表,试作主成分分析。
表5.3 龟甲数据的协方差阵
由于观测资料已被处理为协方差阵,而协方差阵是对称的,只需要输入下三角阵即可,协差阵乘以常数不改变特征向量和累积方差贡献率,所以0.001不用输入。我们采用如下程序
data turtle(type=cov); /*建立数据集turtle*/
_type_='cov'; /*数据集为协方差阵类型*/
input name $ x1-x3; /*建立变量name x1 x2 x3 */
cards; /*以下是数据体*/
x1 11.072 . .
x2 8.019 6.417 .
x3 8.160 6.005 6.773
; /*空语句,结束数据体*/
proc princomp COV; /*用协方差阵计算3个主成分*/
var x1-x3; /*对变量x1 x2 x3求主成分*/
run;
执行后电脑按相关阵做主成分分析;输出主要数表有:协方差阵的特征值(表头为Eigenvalues),特征向量表(表头为Eigenvectors)。解释如下
Eigenvalues
Eigenvalue Difference Proportion Cumulative
PRIN1 23.3035 22.7048 0.960493 0.96049
PRIN2 0.5987 0.2389 0.024676 0.98517
PRIN3 0.3598 . 0.014831 1.00000
上表是特征值表,由表第2列可见,特征值分别是23.303496、0.5986906、0.3598188;由上表第5列可见,第1特征值占总变差的96%,所以只需1个主成分,就能解释全部变化。
Eigenvectors
PRIN1 PRIN2 PRIN3
X1 0.683103 -.158344 -.712950
X2 0.510212 -.595012 0.621002
X3 0.522546 0.787964 0.325666
上是特征向量表,由表可见,第1主成分的系数0.683103、0.510212、0.522546相差不多,所以第1主成分表示龟甲的尺寸的自然对数和,即龟甲体积的自然对数。
5.3 主成分得分
实际问题中常需要知道主成分的值,例如例3中需要知道哪个州犯罪程度严重,哪个州犯罪程度较轻,这就需要计算每个州第一主成分的值;需要知道哪个州暴力犯罪严重,哪个州暴力犯罪较轻,这就需要计算每个州第二主成分的值。同时由于经验正交分解的需要和计算等原因,我们也往往要计算主成分得分。
将各变量值代入主成分的表达式,就能计算主成分的值。例如例2中北京气温的第一主成分是
prin1=0.643587Dec +0.213039Jan+0.735126Feb,
而1951年Dec、Jan、Feb的值分别是1.0、-2.7、-4.3;所以1951年第一主分量的值就是prin1=0.643587*1.0+0.213039*(-2.7)+0.735126*(-4.3)。
定义5.5 当用样本协方差阵求主成分时,求各观测值距平(观测值减去其平均值),再代入主成分的公式,所得称为(协方差阵生成的)主成分得分。
例如例2中第一主成分是0.643587*Dec+0.213039*Jan+0.735126*Feb; Dec,Jan,Feb的样本均值分别是-2.74,-4.59,-2.27;1951年Dec,Jan,Feb的值分别是1.0,-2.7,-4.3;所以1951年(协方差阵生成的)的第一主成分得分就是
0.643587*(1.0+2.74)+0.213039*(-2.7+4.59)+0.735126*(-4.3+2.27)=1.32。
定义5.6 当用样本相关阵阵求主成分时,将各观测标准化(观测值减去其平均值,除以样本标准差)再代入主成分的公式,所得称为(相关阵生成的)主成分得分。
例如例2用相关阵计算时,第一主成分是0.6388*Dec+0.5734*Jan+0.5129*Feb。而1951年标准化的Dec,Jan,Feb的值分别是2.013,1.613,-1.034;于是1951年的(相关阵生成的)第一主成分得分就是
6388*2.013+0.5734*1.613 +0.5129*(-1.034)=1.681
由主成分得分的值很容易算出主成分的值,但由于主成分得分与主成分的值差一常数,因而在比较各次观测主成分的值时,只需比较主成分得分的值即可。
SAS-PRINCOMP过程作主成分分析时,能计算主成分得分,在PROC PRINCOMP语句中加选项OUT=文件名,主成分得分的值即存在该文件中。
例5.4(续) 北京1951~1976年冬季的气温资料,求(协方差阵生成的)各年主成分得分。
解 :采用下列程序
data temperat;
input year Dec Jan Feb;
cards;
1951 1.0 -2.7 -4.3
1952 -5.3 -5.9 -3.5
1953 -2.0 -3.4 -0.8
1954 -5.7 -4.7 -1.1
1955 -0.9 -3.8 -3.1
1956 -5.7 -5.3 -5.9
1957 -2.1 -5.0 -1.6
1958 0.6 -4.3 -0.2
1959 -1.7 -5.7 2.0
1960 -3.6 -3.6 1.3
1961 -3.0 -3.1 -0.8
1962 0.1 -3.9 -1.1
1963 -2.6 -3.0 -5.2
1964 -1.4 -4.9 -1.7
1965 -3.9 -5.7 -2.5
1966 -4.7 -4.8 -3.3
1967 -6.0 -5.6 -4.9
1968 -1.7 -6.4 -5.1
1969 -3.4 -5.6 -2.0
1970 -3.1 -4.2 -2.9
1971 -3.8 -4.9 -3.9
1972 -2.0 -4.1 -2.4
1973 -1.7 -4.2 -2.0
1974 -3.6 -3.3 -2.0
1975 -2.7 -3.7 0.1
1976 -2.4 -7.6 -2.2
;
proc princomp cov out=prin; /*各次观测的主成分值存入数据集prin。*/
var Dec Jan Feb; /* 对变量Dec Jan Feb 作主成分分析*/
proc print data=prin; /* 打印数据集prin所存各次观测的的主成分得分*/
run;
proc sort data=prin; /*将主成分得分按第一主成分得分排序*/
by prin1;
proc print; /* 打印数据集排序后的各次观测的主成分得分*/
run;
proc sort data=prin; /*将主成分得分按第二主成分得分排序*/
by prin2;
proc print; /* 打印数据集排序后的各次观测的主成分得分*/
run;
执行上述程序,与例5.4相比,增加的SAS输出是下表,其中PRIN1、PRIN2、PRIN3分别表示第1、2、3主成分得分。
表5-4 北京冬季气温主成分得分
OBS YEAR DEC JAN FEB PRIN1 PRIN2 PRIN3
1 1951 1.0 -2.7 -4.3 1.32159 4.39464 0.79664
2 1952 -5.3 -5.9 -3.5 -2.82663 -1.23681 -0.48750
3 1953 -2.0 -3.4 -0.8 1.81464 -0.24090 0.88972
4 1954 -5.7 -4.7 -1.1 -1.06412 -2.91502 0.71132
5 1955 -0.9 -3.8 -3.1 0.74659 2.02081 0.25417
6 1956 -5.7 -5.3 -5.9 -4.72054 0.22071 0.26664
7 1957 -2.1 -5.0 -1.6 0.82132 -0.07862 -0.59250
8 1958 0.6 -4.3 -0.2 3.73731 1.02475 -0.73246
9 1959 -1.7 -5.7 2.0 3.57608 -2.36830 -1.47492
10 1960 -3.6 -3.6 1.3 2.28606 -2.83781 1.09907
11 1961 -3.0 -3.1 -0.8 1.23497 -0.89291 1.46318
12 1962 0.1 -3.9 -1.1 2.83912 1.35662 -0.18190
13 1963 -2.6 -3.0 -5.2 -1.72085 2.39084 1.56369
14 1964 -1.4 -4.9 -1.7 1.21963 0.50533 -0.69429
15 1965 -3.9 -5.7 -2.5 -1.14787 -0.88178 -0.72358
16 1966 -4.7 -4.8 -3.3 -2.05911 -0.73417 0.38901
17 1967 -6.0 -5.6 -4.9 -4.24241 -0.72751 0.03805
18 1968 -1.7 -6.4 -5.1 -1.79244 2.30614 -1.95308
19 1969 -3.4 -5.6 -2.0 -0.43721 -0.84625 -0.78440
20 1970 -3.1 -4.2 -2.9 -0.60750 0.24643 0.49508
21 1971 -3.8 -4.9 -3.9 -1.94226 0.29187 0.05198
22 1972 -2.0 -4.1 -2.4 0.48932 0.70789 0.26259
23 1973 -1.7 -4.2 -2.0 0.95514 0.63061 0.07014
24 1974 -3.6 -3.3 -2.0 -0.07594 -0.54455 1.47579
25 1975 -2.7 -3.7 0.1 1.96183 -1.40534 0.77828
26 1976 -2.4 -7.6 -2.2 -0.36673 -0.38669 -2.98072
Obs year Dec Jan Feb Prin1 Prin2 Prin3
1 1956 -5.7 -5.3 -5.9 -4.72054 0.22071 0.26664
2 1967 -6.0 -5.6 -4.9 -4.24241 -0.72751 0.03805
3 1952 -5.3 -5.9 -3.5 -2.82663 -1.23681 -0.48750
4 1966 -4.7 -4.8 -3.3 -2.05911 -0.73417 0.38901
5 1971 -3.8 -4.9 -3.9 -1.94226 0.29187 0.05198
6 1968 -1.7 -6.4 -5.1 -1.79244 2.30614 -1.95308
7 1963 -2.6 -3.0 -5.2 -1.72085 2.39084 1.56369
8 1965 -3.9 -5.7 -2.5 -1.14787 -0.88178 -0.72358
9 1954 -5.7 -4.7 -1.1 -1.06412 -2.91502 0.71132
10 1970 -3.1 -4.2 -2.9 -0.60750 0.24643 0.49508
11 1969 -3.4 -5.6 -2.0 -0.43721 -0.84625 -0.78440
12 1976 -2.4 -7.6 -2.2 -0.36673 -0.38669 -2.98072
13 1974 -3.6 -3.3 -2.0 -0.07594 -0.54455 1.47579
14 1972 -2.0 -4.1 -2.4 0.48932 0.70789 0.26259
15 1955 -0.9 -3.8 -3.1 0.74659 2.02081 0.25417
16 1957 -2.1 -5.0 -1.6 0.82132 -0.07862 -0.59250
17 1973 -1.7 -4.2 -2.0 0.95514 0.63061 0.07014
18 1964 -1.4 -4.9 -1.7 1.21963 0.50533 -0.69429
19 1961 -3.0 -3.1 -0.8 1.23497 -0.89291 1.46318
20 1951 1.0 -2.7 -4.3 1.32159 4.39464 0.79664
21 1953 -2.0 -3.4 -0.8 1.81464 -0.24090 0.88972
22 1975 -2.7 -3.7 0.1 1.96183 -1.40534 0.77828
23 1960 -3.6 -3.6 1.3 2.28606 -2.83781 1.09907
24 1962 0.1 -3.9 -1.1 2.83912 1.35662 -0.18190
25 1959 -1.7 -5.7 2.0 3.57608 -2.36830 -1.47492
26 1958 0.6 -4.3 -0.2 3.73731 1.02475 -0.73246
Obs year Dec Jan Feb Prin1 Prin2 Prin3
1 1954 -5.7 -4.7 -1.1 -1.06412 -2.91502 0.71132
2 1960 -3.6 -3.6 1.3 2.28606 -2.83781 1.09907
3 1959 -1.7 -5.7 2.0 3.57608 -2.36830 -1.47492
4 1975 -2.7 -3.7 0.1 1.96183 -1.40534 0.77828
5 1952 -5.3 -5.9 -3.5 -2.82663 -1.23681 -0.48750
6 1961 -3.0 -3.1 -0.8 1.23497 -0.89291 1.46318
7 1965 -3.9 -5.7 -2.5 -1.14787 -0.88178 -0.72358
8 1969 -3.4 -5.6 -2.0 -0.43721 -0.84625 -0.78440
9 1966 -4.7 -4.8 -3.3 -2.05911 -0.73417 0.38901
10 1967 -6.0 -5.6 -4.9 -4.24241 -0.72751 0.03805
11 1974 -3.6 -3.3 -2.0 -0.07594 -0.54455 1.47579
12 1976 -2.4 -7.6 -2.2 -0.36673 -0.38669 -2.98072
13 1953 -2.0 -3.4 -0.8 1.81464 -0.24090 0.88972
14 1957 -2.1 -5.0 -1.6 0.82132 -0.07862 -0.59250
15 1956 -5.7 -5.3 -5.9 -4.72054 0.22071 0.26664
16 1970 -3.1 -4.2 -2.9 -0.60750 0.24643 0.49508
17 1971 -3.8 -4.9 -3.9 -1.94226 0.29187 0.05198
18 1964 -1.4 -4.9 -1.7 1.21963 0.50533 -0.69429
19 1973 -1.7 -4.2 -2.0 0.95514 0.63061 0.07014
20 1972 -2.0 -4.1 -2.4 0.48932 0.70789 0.26259
21 1958 0.6 -4.3 -0.2 3.73731 1.02475 -0.73246
22 1962 0.1 -3.9 -1.1 2.83912 1.35662 -0.18190
23 1955 -0.9 -3.8 -3.1 0.74659 2.02081 0.25417
24 1968 -1.7 -6.4 -5.1 -1.79244 2.30614 -1.95308
25 1963 -2.6 -3.0 -5.2 -1.72085 2.39084 1.56369
26 1951 1.0 -2.7 -4.3 1.32159 4.39464 0.79664
以上是26年观测资料与每年没排序和排序后的第1,2,3主成分得分,分别用变量prin1,prin2,prin3表示。
从主成分得分表可见1958,1959年第1主成分最强(冬季温度偏高),1956年第1主成分负方向最强(冬季温度偏低),1974年第1主成分绝对值最小(冬季温度最接近常年);1951年第2主成分最强(12月比2月温度高得多),1954年第2主成分负方向最强(12月温度比2月温度低得多),1957年第2主成分绝对值最小(12月温度与2月温度差接近历年平均值)。
例5.5(续) 对于例5-5美国各州犯罪率数据求主成分得分。为了比较各州犯罪轻重程度和比较各州暴力犯罪轻重程度,分别按第一、第二主成分得分排序后输出,采用下列程序:
data crime; /*建立数据集crime*/
input state $ 1-15 murder rape robbery assult burglary larceny auto;
/*建立变量state murder rape robbery assult burglary larceny auto。state $ 1-15表示前15列存州名。murder rape robbery assult burglary larceny auto 表7种罪的犯罪率*/
cards; /*以下为数据体*/
Albama 14.2 25.2 96.8 278.3 1135.5 1881.9 280.7
Alaska 10.8 51.6 96.8 284.0 1331.7 3369.8 753.3
Arirona 9.5 34.2 138.2 312.3 2346.1 4467.4 439.5
Arkansas 8.8 34.2 138.2 312.3 2346.1 4467.4 439.5
Califonia 11.5 49.4 287.0 358.0 2139.4 3499.8 663.5
Colorado 6.3 42.0 170.7 292.9 1935.2 3903.2 477.1
Conecticat 4.2 16.8 129.5 131.8 1346.0 2620.7 593.2
Delaware 6.0 24.9 157.0 194.2 1682.6 3678.4 467.0
Florida 10.2 39.6 187.9 449.1 1859.9 3840.5 351.4
Geogia 11.7 31.1 140.5 256.5 1351.1 2170.2 297.9
Hawaii 7.2 25.5 128.0 64.1 1911.5 3920.4 489.4
Idaho 5.5 19.4 39.6 172.5 1050.8 2599.6 237.6
Illinois 9.9 21.8 211.3 209.0 1085.0 2828.5 528.6
Indiana 7.4 26.5 123.2 153.5 1086.2 2498.7 377.4
Iowa 2.3 10.6 41.2 89.8 812.5 2685.1 219.9
Kansas 6.6 22.0 100.7 180.5 1270.4 2739.3 244.3
Kentaky 10.1 19.1 81.1 123.3 872.2 1662.1 245.4
Loisana 15.5 30.9 142.9 335.5 1165.5 2469.9 337.7
Maine 2.4 13.5 38.7 170.0 1253.1 2350.7 246.9
Maryland 8.0 34.8 292.1 358.9 1400.0 3177.7 428.5
Masschusetts 3.1 20.8 169.1 231.6 1532.2 2311.3 1140.1
Michigan 9.3 38.9 261.9 274.6 1522.7 3159.0 545.5
Minnesota 2.7 19.5 85.9 85.8 1134.7 2559.3 343.1
Mississippi 14.3 19.6 65.7 189.1 915.6 1239.9 144.4
Missouri 9.6 28.3 189.0 233.5 1318.3 2424.2 378.4
Montana 5.4 16.7 39.2 156.8 804.9 2773.2 309.3
Nebraska 3.9 18.1 64.7 112.7 760.0 2316.1 249.1
Nevada 15.8 49.1 323.1 355.0 2453.1 4212.6 559.2
Mew Hampashare 3.2 10.7 23.2 76.0 1041.7 2343.9 293.4
New Jersey 5.6 21.0 180.4 185.1 1435.8 2774.5 511.5
New Maxico 8.8 39.1 109.6 343.4 1418.7 3008.6 259.5
New York 10.7 29.4 472.6 319.1 1728.0 2782.0 745.8
North Carolina 10.6 17.0 61.3 318.3 1154.1 2037.8 192.1
North Dakoda 100.9 9.0 13.3 43.8 446.1 1843.0 144.7
Ohio 7.8 27.3 190.5 181.1 1216.0 2696.8 400.4
Oklahoma 8.6 29.2 73.8 205.0 1288.2 2228.1 326.8
Oregan 4.9 39.9 124.1 286.9 1636.4 3506.1 388.9
Pennsyvania 5.6 19.0 130.3 128.0 877.5 1624.1 333.2
Rhode Island 3.6 10.5 86.5 201.0 1849.5 2844.1 791.4
South Carolina 11.9 33.0 105.9 485.3 1613.6 2342.4 245.1
South Dakoda 2.0 13.5 17.9 155.7 570.5 1704.4 147.5
Tennessee 10.1 29.7 145.8 203.9 1259.7 1776.5 314.0
Texas 13.3 33.8 152.4 208.2 1603.1 2988.7 397.6
Utah 3.5 20.3 68.8 147.3 1171.6 3004.6 334.5
Vermont 1.4 15.9 30.8 101.2 1348.2 2201.0 265.2
Virginia 9.0 23.3 92.1 165.7 986.2 2521.2 226.7
Wasinton 4.3 39.6 106.2 224.8 1605.6 3386.9 360.3
West Viginia 6.0 13.2 42.2 90.9 597.4 1341.7 163.3
Wiskonsin 2.8 12.9 52.2 63.7 846.9 2614.2 220.7
Wyoming 5.4 21.9 39.7 173.9 811.6 2772.2 282.0
;
proc princomp out=crimprin n=2;
var murder rape robbery assult burglary larceny auto;
run;
proc sort data=crimprin out=sortprin1;
by prin1; /*将数据集crimprin第一主成分按照从小到大顺序,重新排序*/
run;
proc print data=sortprin1;
var state prin1; /*打印第一主成分和州名*/
run;
proc sort data=crimprin out=sortprin2;
by prin2; /*将数据集crimprin第二主成分按照从小到大顺序,从新排序*/
run;
proc print data=sortprin2;
var state prin2; /*打印第二主成分和州名*/
run;
执行此程序计算了2个主成分得分,存入数据集crimprin,分别记为 PRIN1、PRIN2。按主成分大小顺序相同将第1主成分得分从小到大排序后打印,得到的输出是
表 5-5美国各州犯罪第一主成分得分排序表(从小到大)
Obs state Prin1
1 North Dakoda -4.24229
2 West Viginia -3.18590
3 South Dakoda -2.86744
4 Mississippi -2.26684
……………………
47 Arkansas 2.83175
48 New York 3.18532
49 Califonia 3.94339
50 Nevada 4.56461
由上表可见North Dakoda州的犯罪程度最轻,Nevada州的犯罪程度最重。
将第2主成分得分排序后打印,得到的主要输出是
表5-6 美国各州犯罪第二主成分得分排序表(从小到大)
Obs state Prin2
1 Masschusetts -1.44958
2 Rhode Island -1.40935
3 Conecticat -1.00923
4 Mew Hampashare -0.97043
. . .
47 Nevada 0.90345
48 Loisana 0.93489
49 South Carolina 1.15880
50 North Dakoda 5.60596
由上表可见Masschusetts州的暴力犯罪趋势最轻,North Dakoda州的暴力犯罪趋势最重。
5.4 主成分聚类与主成分回归
5.4.1 变量聚类与样品分类
主成分分析可用于聚类:变量聚类与样品聚类。
变量聚类:由主成分系数的差异,可将变量聚类。例如例5.5中第2主成分中murder,rape, assult系数为负的, burglary,larceny, auto系数是正的。按系数正负可把7个变量分为两类: murder, rape, assult属于暴力程度严重的一类;burglary,larceny,auto属于暴力程度较轻的一类。按照这种方法,根据主成分系数的正负可以将变量聚类。
样品聚类:如果2个主成分能很好的概括随机向量的信息,计算每个样品的这两个主成分得分,把他们的散点图画出来,就能从图上将样品分类。
例5.5(续2) 按照第一、第二主成分得分,画出散点图
data crime; /*建立数据集crime*/
input state $ 1-15 murder rape robbery assult burglary larceny auto;
/*建立变量state murder rape robbery assult burglary larceny auto。state $ 1-15表示前15列存州名。murder rape robbery assult burglary larceny auto 表7种罪的犯罪率*/
cards; /*以下为数据体*/
Albama 14.2 25.2 96.8 278.3 1135.5 1881.9 280.7
Alaska 10.8 51.6 96.8 284.0 1331.7 3369.8 753.3
Arirona 9.5 34.2 138.2 312.3 2346.1 4467.4 439.5
Arkansas 8.8 34.2 138.2 312.3 2346.1 4467.4 439.5
Califonia 11.5 49.4 287.0 358.0 2139.4 3499.8 663.5
Colorado 6.3 42.0 170.7 292.9 1935.2 3903.2 477.1
Conecticat 4.2 16.8 129.5 131.8 1346.0 2620.7 593.2
Delaware 6.0 24.9 157.0 194.2 1682.6 3678.4 467.0
Florida 10.2 39.6 187.9 449.1 1859.9 3840.5 351.4
Geogia 11.7 31.1 140.5 256.5 1351.1 2170.2 297.9
Hawaii 7.2 25.5 128.0 64.1 1911.5 3920.4 489.4
Idaho 5.5 19.4 39.6 172.5 1050.8 2599.6 237.6
Illinois 9.9 21.8 211.3 209.0 1085.0 2828.5 528.6
Indiana 7.4 26.5 123.2 153.5 1086.2 2498.7 377.4
Iowa 2.3 10.6 41.2 89.8 812.5 2685.1 219.9
Kansas 6.6 22.0 100.7 180.5 1270.4 2739.3 244.3
Kentaky 10.1 19.1 81.1 123.3 872.2 1662.1 245.4
Loisana 15.5 30.9 142.9 335.5 1165.5 2469.9 337.7
Maine 2.4 13.5 38.7 170.0 1253.1 2350.7 246.9
Maryland 8.0 34.8 292.1 358.9 1400.0 3177.7 428.5
Masschusetts 3.1 20.8 169.1 231.6 1532.2 2311.3 1140.1
Michigan 9.3 38.9 261.9 274.6 1522.7 3159.0 545.5
Minnesota 2.7 19.5 85.9 85.8 1134.7 2559.3 343.1
Mississippi 14.3 19.6 65.7 189.1 915.6 1239.9 144.4
Missouri 9.6 28.3 189.0 233.5 1318.3 2424.2 378.4
Montana 5.4 16.7 39.2 156.8 804.9 2773.2 309.3
Nebraska 3.9 18.1 64.7 112.7 760.0 2316.1 249.1
Nevada 15.8 49.1 323.1 355.0 2453.1 4212.6 559.2
Mew Hampashare 3.2 10.7 23.2 76.0 1041.7 2343.9 293.4
New Jersey 5.6 21.0 180.4 185.1 1435.8 2774.5 511.5
New Maxico 8.8 39.1 109.6 343.4 1418.7 3008.6 259.5
New York 10.7 29.4 472.6 319.1 1728.0 2782.0 745.8
North Carolina 10.6 17.0 61.3 318.3 1154.1 2037.8 192.1
North Dakoda 100.9 9.0 13.3 43.8 446.1 1843.0 144.7
Ohio 7.8 27.3 190.5 181.1 1216.0 2696.8 400.4
Oklahoma 8.6 29.2 73.8 205.0 1288.2 2228.1 326.8
Oregan 4.9 39.9 124.1 286.9 1636.4 3506.1 388.9
Pennsyvania 5.6 19.0 130.3 128.0 877.5 1624.1 333.2
Rhode Island 3.6 10.5 86.5 201.0 1849.5 2844.1 791.4
South Carolina 11.9 33.0 105.9 485.3 1613.6 2342.4 245.1
South Dakoda 2.0 13.5 17.9 155.7 570.5 1704.4 147.5
Tennessee 10.1 29.7 145.8 203.9 1259.7 1776.5 314.0
Texas 13.3 33.8 152.4 208.2 1603.1 2988.7 397.6
Utah 3.5 20.3 68.8 147.3 1171.6 3004.6 334.5
Vermont 1.4 15.9 30.8 101.2 1348.2 2201.0 265.2
Virginia 9.0 23.3 92.1 165.7 986.2 2521.2 226.7
Wasinton 4.3 39.6 106.2 224.8 1605.6 3386.9 360.3
West Viginia 6.0 13.2 42.2 90.9 597.4 1341.7 163.3
Wiskonsin 2.8 12.9 52.2 63.7 846.9 2614.2 220.7
Wyoming 5.4 21.9 39.7 173.9 811.6 2772.2 282.0
;
proc princomp out=crimprin n=2;
var murder rape robbery assult burglary larceny auto;
run;
PROC PLOT data=crimprin;
PLOT PRIN2*PRIN1=STATE/VPOS=31;
TITLE2 ‘PLOT OF THE FIRST TWO PRINCIPAL COMPONENTS’;
RUN;
例5.7 (气温分析) 本例的输入资料文件(TEMPERAT)是美国六十四个城市一月与七月的平均日温。
DATA TEMPERAT;
TITLE2 'MEAN TEMPERATURE IN JANUARY AND JULY FOR SELECTED CITIES';
INPUT CITY $1-15 JANUARY JULY;
CARDS;
MOBILE 51.2 81.6
PHOENIX 51.2 91.2
LITTLE ROCK 39.5 81.4
SACRAMENTO 45.1 75.2
DENVER 29.9 73.0
HARTFORD 24.8 72.7
WILMINGTON 32.0 75.8
WASHINGTON DC 35.6 78.7
JACKSONVILLE 54.6 81.0
MIAMI 67.2 82.3
ATLANTA 42.4 78.0
BOISE 29.0 74.5
CHICAGO 22.9 71.9
PEORIA 23.8 75.1
INDIANAPOLIS 27.9 75.0
DES MOINES 19.4 75.1
WICHITA 31.3 80.7
LOUISVILLE 33.3 76.9
NEW ORLEANS 52.9 81.9
PORTLAND, MAINE 21.5 68.0
BALTIMORE 33.4 76.6
BOSTON 29.2 73.3
DETROIT 25.5 73.3
SAULT STE MARIE 14.2 63.8
DULUTH 8.5 65.6
MINNEAPOLIS 12.2 71.9
JACKSON 47.1 81.7
KANSAS CITY 27.8 78.8
ST LOUIS 31.3 78.6
GREAT FALLS 20.5 69.3
OMAHA 22.6 77.2
RENO 31.9 69.3
CONCORD 20.6 69.7
ATLANTIC CITY 32.7 75.1
ALBUQUERQUE 35.2 78.7
ALBANY 21.5 72.0
BUFFALO 23.7 70.1
NEW YORK 32.2 76.6
CHARLOTTE 42.1 78.5
RALEIGH 40.5 77.5
BISMARCK 8.2 70.8
CINCINNATI 31.1 75.6
CLEVELAND 26.9 71.4
COLUMBUS 28.4 73.6
OKLAHOMA CITY 36.8 81.5
PORTLAND, OREG 38.1 67.1
PHILADELPHIA 32.3 76.8
PITTSBURGH 28.1 71.9
PROVIDENCE 28.4 72.1
COLUMBIA 45.4 81.2
SIOUX FALLS 14.2 73.3
MEMPHIS 40.5 79.6
NASHVILLE 38.3 79.6
DALLAS 44.8 84.8
EL PASO 43.6 82.3
HOUSTON 52.1 83.3
SALT LAKE CITY 28.0 76.7
BURLINGTON 16.8 69.8
NORFOLK 40.5 78.3
RICHMOND 37.5 77.9
SPOKANE 25.4 69.7
CHARLESTON, WV 34.5 75.0
MILWAUKEE 19.4 69.9
CHEYENNE 26.6 69.1
;
PROC PLOT;
PLOT JULY*JANUARY=CITY/VPOS=36;
PROC PRINCOMP COV OUT=PRIN;
VAR JULY JANUARY;
PROC PLOT;
PLOT PRIN2*PRIN1=CITY/VPOS=26;
TITLE3 'PLOT OF PRINCIPAL COMPONENTS';
Run;
例5.8 美国大学生篮球队排名
data bballm;
label csn = 'Community Sports News (Chapel Hill NC)'
dursun = 'Durham Sun'
durher = 'Durham Morning Herald'
waspost = 'Washington Post'
usatoda = 'USA Today'
spormag = 'Sport Magazine'
insport = 'Inside Sports'
upi = 'United Press International'
ap = 'Associated Press'
sporill = 'Sports Illustrated'
;
title1 'Pre-Season 1985 College Basketball Rankings';
input school $13. csn dursun durher waspost usatoda
spormag insport upi ap sporill;
format csn -- sporill 5.1;
cards;
Louisville 1 8 1 9 8 9 6 10 9 9
Georgia Tech 2 2 4 3 1 1 1 2 1 1
Kansas 3 4 5 1 5 11 8 4 5 7
Michigan 4 5 9 4 2 5 3 1 3 2
Duke 5 6 7 5 4 10 4 5 6 5
UNC 6 1 2 2 3 4 2 3 2 3
Syracuse 7 10 6 11 6 6 5 6 4 10
Notre Dame 8 14 15 13 11 20 18 13 12 .
Kentucky 9 15 16 14 14 19 11 12 11 13
LSU 10 9 13 . 13 15 16 9 14 8
DePaul 11 . 21 15 20 . 19 . . 19
Georgetown 12 7 8 6 9 2 9 8 8 4
Navy 13 20 23 10 18 13 15 . 20 .
Illinois 14 3 3 7 7 3 10 7 7 6
Iowa 15 16 . . 23 . . 14 . 20
Arkansas 16 . . . 25 . . . . 16
Memphis State 17 . 11 . 16 8 20 . 15 12
Washington 18 . . . . . . 17 . .
UAB 19 13 10 . 12 17 . 16 16 15
UNLV 20 18 18 19 22 . 14 18 18 .
NC State 21 17 14 16 15 . 12 15 17 18
Maryland 22 . . . 19 . . . 19 14
Pittsburg 23 . . . . . . . . .
Oklahoma 24 19 17 17 17 12 17 . 13 17
Indiana 25 12 20 18 21 . . . . .
Virginia 26 . 22 . . 18 . . . .
Old Dominion 27 . . . . . . . . .
Auburn 28 11 12 8 10 7 7 11 10 11
St. Johns 29 . . . . 14 . . . .
UCLA 30 . . . . . . 19 . .
St. Joseph's . . 19 . . . . . . .
Tennessee . . 24 . . 16 . . . .
Montana . . . 20 . . . . . .
Houston . . . . 24 . . . . .
Virginia Tech . . . . . . 13 . . .
;
/* Proc MEANS is used to output a dataset containing the maximum */
/* value of each of the newspaper and magazine rankings. The */
/* output data set, maxrank, is then used to set the missing */
/* values to the highest rank plus thirty-six, divided by two */
/* (that is, the mean of the missing ranks). This ad hoc method of */
/* replacing missing values is based more on intuition than on */
/* rigorous statistical theory. Observations are weighted by the */
/* number of nonmissing values. */
proc means data=bballm;
output out=maxrank
max=mcsn mdurs mdurh mwas musa mspom mins mupi map mspoi;
data bball;
set bballm;
if _n_=1 then set maxrank;
array services{10} csn -- sporill;
array maxranks{10} mcsn -- mspoi;
keep school csn -- sporill weight;
weight=0;
do i=1 to 10;
if services{i}=. then services{i}=(maxranks{i}+36) / 2;
else weight=weight+1;
end;
/* Use the PRINCOMP procedure to transform the observed ranks. */
/* Use n=1 because the data should be related to a single */
/* underlying variable. Sort the data and print the */
/* resulting component. */
proc princomp data=bball n=1 out=pcbball standard;
var csn --sporill;
weight weight;
proc sort data=pcbball;
by prin1;
proc print;
var school prin1;
title2 'College Teams as Ordered by PRINCOMP';
run;
例5.9 55个地区或国家的赛跑纪录如表5-7,试作主成分分析,并将55个国家或地区按赛跑成绩分类。
表 5-7 55个地区或国家的赛跑纪录
序号
国家或地区
100m
(秒)
200m
(秒)
400m
(秒)
800m
(分)
1500m
(分)
5000m
(分)
10000m(分)
马拉松(分)
1
argentin
10.39
20.81
46.84
1.81
3.70
14.04
29.36
137.72
2
australi
10.31
20.06
44.84
1.74
3.57
13.28
27.66
128.30
3
austria
10.44
20.81
46.82
1.79
3.60
13.26
27.72
135.90
4
belgium
10.34
20.68
45.04
1.73
3.60
13.22
27.45
129.95
5
bermuda
10.28
20.58
45.91
1.80
3.75
14.68
30.55
146.62
6
brazil
10.22
20.43
45.21
1.73
3.66
13.62
28.62
133.13
7
burma
10.64
21.52
48.30
1.80
3.85
14.45
30.28
139.95
8
canada
10.17
20.22
45.68
1.76
3.63
13.55
28.09
130.15
9
chile
10.34
20.80
46.20
1.79
3.71
13.61
29.30
134.03
10
china
10.51
21.04
47.30
1.81
3.73
13.90
29.13
133.53
11
columbia
10.43
21.05
46.10
1.82
3.74
13.49
27.88
131.35
12
cookis
12.18
23.20
52.94
2.02
4.24
16.70
35.38
164.70
13
costa
10.94
21.90
48.66
1.87
3.84
14.03
28.81
136.58
14
czech
10.35
20.65
45.64
1.76
3.58
13.42
28.19
134.32
15
denmark
10.56
20.52
45.89
1.78
3.61
13.50
28.11
130.78
16
domrep
10.14
20.65
46.80
1.82
3.82
14.91
31.45
154.12
17
finland
10.43
20.69
45.49
1.74
3.61
13.27
27.52
130.87
18
france
10.11
20.38
45.28
1.73
3.57
13.34
27.97
132.30
19
gdr
10.12
20.33
44.87
1.73
3.56
13.17
27.42
129.92
20
frg
10.16
20.37
44.50
1.73
3.53
13.21
27.61
132.23
21
gbni
10.11
20.21
44.93
1.70
3.51
13.01
27.51
129.13
22
greece
10.22
20.71
46.56
1.78
3.64
14.59
28.45
134.60
23
guatemal
10.98
21.82
48.40
1.89
3.80
14.16
30.11
139.33
24
hungary
10.26
20.62
46.02
1.77
3.62
13.49
28.44
132.58
25
india
10.60
21.42
45.73
1.76
3.73
13.77
28.81
131.98
26
indonesi
10.59
21.49
47.80
1.84
3.92
14.73
30.79
148.83
27
ireland
10.61
20.96
46.30
1.79
3.56
13.32
27.81
132.35
28
israel
10.71
21.00
47.80
1.77
3.72
13.66
28.93
137.55
29
italy
10.01
19.72
45.26
1.73
3.60
13.23
27.52
131.08
30
japan
10.34
20.81
45.86
1.79
3.64
13.41
27.72
128.63
31
kenya
10.46
20.66
44.92
1.73
3.55
13.10
27.38
129.75
32
korea
10.34
20.89
46.90
1.79
3.77
13.96
29.23
136.25
33
dprkorea
10.91
21.94
47.30
1.85
3.77
14.13
29.67
130.87
34
luxembou
10.35
20.77
47.40
1.82
3.67
13.64
29.08
141.27
35
malaysia
10.40
20.92
46.30
1.82
3.80
14.64
31.01
154.10
36
mauritiu
11.19
22.45
47.70
1.88
3.83
15.06
31.77
152.23
37
mexico
10.42
21.30
46.10
1.80
3.65
13.46
27.95
129.20
38
netherla
10.52
20.95
45.10
1.74
3.62
13.36
27.61
129.02
39
nz
10.51
20.88
46.10
1.74
3.54
13.21
27.70
128.98
40
norway
10.55
21.16
46.71
1.76
3.62
13.34
27.69
131.48
41
png
10.96
21.78
47.90
1.90
4.01
14.72
31.36
148.22
42
philippi
10.78
21.64
46.24
1.81
3.83
14.74
30.64
145.27
43
poland
10.16
20.24
45.36
1.76
3.60
13.29
27.89
131.58
44
portugal
10.53
21.17
46.70
1.79
3.62
13.13
27.38
128.65
45
rumania
10.41
20.98
45.87
1.76
3.64
13.25
27.67
132.50
46
singapor
10.38
21.28
47.40
1.88
3.89
15.11
31.32
157.77
47
spain
10.42
20.77
45.98
1.76
3.55
13.31
27.73
131.57
48
sweden
10.25
20.61
45.63
1.77
3.61
13.29
27.94
130.63
49
switzerl
10.37
20.46
45.78
1.78
3.55
13.22
27.91
131.20
50
taipei
10.59
21.29
46.80
1.79
3.77
14.07
30.07
139.27
51
thailand
10.39
21.09
47.91
1.83
3.84
15.23
32.56
149.90
52
turkey
10.71
21.43
47.60
1.79
3.67
13.56
28.58
131.50
53
usa
9.93
19.75
43.86
1.73
3.53
13.20
27.43
128.22
54
ussr
10.07
20.00
44.60
1.75
3.59
13.20
27.53
130.55
55
wsamoa
10.82
21.86
49.00
2.02
4.24
16.28
34.71
161.83
可用下列SAS程序作主成分分析,并将第1,2主成分画散点图
data runrecod;
input country $ x1-x8;
cards;
argentin 10.39 20.81 46.84 1.81 3.70 14.04 29.36 137.72
australi 10.31 20.06 44.84 1.74 3.57 13.28 27.66 128.30
austria 10.44 20.81 46.82 1.79 3.60 13.26 27.72 135.90
belgium 10.34 20.68 45.04 1.73 3.60 13.22 27.45 129.95
bermuda 10.28 20.58 45.91 1.80 3.75 14.68 30.55 146.62
brazil 10.22 20.43 45.21 1.73 3.66 13.62 28.62 133.13
burma 10.64 21.52 48.30 1.80 3.85 14.45 30.28 139.95
canada 10.17 20.22 45.68 1.76 3.63 13.55 28.09 130.15
chile 10.34 20.80 46.20 1.79 3.71 13.61 29.30 134.03
china 10.51 21.04 47.30 1.81 3.73 13.90 29.13 133.53
columbia 10.43 21.05 46.10 1.82 3.74 13.49 27.88 131.35
cookis 12.18 23.20 52.94 2.02 4.24 16.70 35.38 164.70
costa 10.94 21.90 48.66 1.87 3.84 14.03 28.81 136.58
czech 10.35 20.65 45.64 1.76 3.58 13.42 28.19 134.32
denmark 10.56 20.52 45.89 1.78 3.61 13.50 28.11 130.78
domrep 10.14 20.65 46.80 1.82 3.82 14.91 31.45 154.12
finland 10.43 20.69 45.49 1.74 3.61 13.27 27.52 130.87
france 10.11 20.38 45.28 1.73 3.57 13.34 27.97 132.30
gdr 10.12 20.33 44.87 1.73 3.56 13.17 27.42 129.92
frg 10.16 20.37 44.50 1.73 3.53 13.21 27.61 132.23
gbni 10.11 20.21 44.93 1.70 3.51 13.01 27.51 129.13
greece 10.22 20.71 46.56 1.78 3.64 14.59 28.45 134.60
guatemal 10.98 21.82 48.40 1.89 3.80 14.16 30.11 139.33
hungary 10.26 20.62 46.02 1.77 3.62 13.49 28.44 132.58
india 10.60 21.42 45.73 1.76 3.73 13.77 28.81 131.98
indonesi 10.59 21.49 47.80 1.84 3.92 14.73 30.79 148.83
ireland 10.61 20.96 46.30 1.79 3.56 13.32 27.81 132.35
israel 10.71 21.00 47.80 1.77 3.72 13.66 28.93 137.55
italy 10.01 19.72 45.26 1.73 3.60 13.23 27.52 131.08
japan 10.34 20.81 45.86 1.79 3.64 13.41 27.72 128.63
kenya 10.46 20.66 44.92 1.73 3.55 13.10 27.38 129.75
korea 10.34 20.89 46.90 1.79 3.77 13.96 29.23 136.25
dprkorea 10.91 21.94 47.30 1.85 3.77 14.13 29.67 130.87
luxembou 10.35 20.77 47.40 1.82 3.67 13.64 29.08 141.27
malaysia 10.40 20.92 46.30 1.82 3.80 14.64 31.01 154.10
mauritiu 11.19 22.45 47.70 1.88 3.83 15.06 31.77 152.23
mexico 10.42 21.30 46.10 1.80 3.65 13.46 27.95 129.20
netherla 10.52 20.95 45.10 1.74 3.62 13.36 27.61 129.02
nz 10.51 20.88 46.10 1.74 3.54 13.21 27.70 128.98
norway 10.55 21.16 46.71 1.76 3.62 13.34 27.69 131.48
png 10.96 21.78 47.90 1.90 4.01 14.72 31.36 148.22
philippi 10.78 21.64 46.24 1.81 3.83 14.74 30.64 145.27
poland 10.16 20.24 45.36 1.76 3.60 13.29 27.89 131.58
portugal 10.53 21.17 46.70 1.79 3.62 13.13 27.38 128.65
rumania 10.41 20.98 45.87 1.76 3.64 13.25 27.67 132.50
singapor 10.38 21.28 47.40 1.88 3.89 15.11 31.32 157.77
spain 10.42 20.77 45.98 1.76 3.55 13.31 27.73 131.57
sweden 10.25 20.61 45.63 1.77 3.61 13.29 27.94 130.63
switzerl 10.37 20.46 45.78 1.78 3.55 13.22 27.91 131.20
taipei 10.59 21.29 46.80 1.79 3.77 14.07 30.07 139.27
thailand 10.39 21.09 47.91 1.83 3.84 15.23 32.56 149.90
turkey 10.71 21.43 47.60 1.79 3.67 13.56 28.58 131.50
usa 9.93 19.75 43.86 1.73 3.53 13.20 27.43 128.22
ussr 10.07 20.00 44.60 1.75 3.59 13.20 27.53 130.55
wsamoa 10.82 21.86 49.00 2.02 4.24 16.28 34.71 161.83
proc princomp out=w;
var x1-x8;
run;
由执行程序后得到的表可作如下分析。以下是特征值表
Eigenvalues
Eigenvalue Difference Proportion Cumulative
PRIN1 6.62215 5.74453 0.827768 0.82777
PRIN2 0.87762 0.71830 0.109702 0.93747
PRIN3 0.15932 0.03527 0.019915 0.95739
PRIN4 0.12405 0.04417 0.015506 0.97289
PRIN5 0.07988 0.01192 0.009985 0.98288
PRIN6 0.06797 0.02155 0.008496 0.99137
PRIN7 0.04642 0.02382 0.005802 0.99717
PRIN8 0.02260 . 0.002825 1.00000
两个主成分的累计方差贡献率为0.9375,可见只要两个主成分就够了。以下是由特征向量表
Eigenvectors
PRIN1 PRIN2 PRIN3 PRIN4 PRIN5 PRIN6 PRIN7 PRIN8
X1 0.317556 0.566878 0.332262 0.127628 0.262555 -.593704 0.136241 0.105542
X2 0.336979 0.461626 0.360657 -.259116 -.153957 0.656137 -.112640 -.096054
X3 0.355645 0.248273-.560467 0.652341 -.218323 0.156625 -.002854 -.000127
X4 0.368684 0.012430 -.532482 -.479999 0.540053 -.014692 -.238016 -.038165
X5 0.372810 -.139797 -.153443 -.404510 -.487715 -.157843 0.610011 0.139291
X6 0.364374 -.312030 0.189764 0.029588 -.253979 -.141299 -.591299 0.546697
X7 0.366773 -.306860 0.181752 0.080069 -.133176 -.219017 -.176871 -.796795
X8 0.341926 -.438963 0.263209 0.299512 0.497928 0.315285 0.398822 0.158164
第1主成分个分量系数是大体相同的正数,第1主成分反映各国或地区跑步实力,它越小速度快,实力越强。第2主成分在100米,200米,400米的系数为正,1500米,5000米,10000米和马拉松系数为负,且距离越长,负的越多,反映短距离跑速度与长距离跑速度差异,即长距离跑相对短距离跑的实力,它越小短距离相对长距离实力越强。
为了将第1主成分作为横轴,第2主成分作为立轴画出散点图,可用下列程序
proc gplot data=w;
plot prin2*prin1=country;
run;
得图5.1
图5.1 径塞
记录
混凝土 养护记录下载土方回填监理旁站记录免费下载集备记录下载集备记录下载集备记录下载
第1,2主成分散点图
从图5.1可见这些国家或地区可分为4类:第1类只有库克群岛,第1,2主成分都最大,说明跑步速度慢,长距离跑相对短距离跑的实力大(短距离跑速度特别慢);第2类包括百慕大,多米尼加,马来西亚,新加坡,泰国第1主成分较大,第2主成分为负值,说明总体来说,跑步速度不快,短距离跑相对长距离跑的实力强,即短距离跑得快而长距离跑得慢;第3类包括中国台北,朝鲜,哥斯达尼加,巴布亚新几内亚,危地马拉,毛里求斯,缅甸,菲律宾,印度尼西亚,第1主成分较大,说明总体来说,跑步速度不快,第2主成分为正值,说明短距离跑相对长距离跑的实力小(长距离跑速度较快);其他国家或地区为一类;第1主成分较小,说明跑步速度总的较快,第2主成分为正值,说明短距离跑得快而长距离跑得慢。
5.4.2 主成分回归
在回归分析中,常遇到自变量存在多重共线性问题,即自变量的观测值存在线性相关,或近似线性相关。这时设计矩阵满足
,第4章中已指出,用公式
估计参数会造成较大方差。选取自变量的主成分,它们是彼此正交的。用少量主成分作回归,再将主成分化为原始变量,这样得到的回归方程就不存在较大方差了。与逐步回归相比,这样做的好处是,所有变量都被计算,而逐步回归略去一些自变量,这些自变量不被计算。主成分回归原理见例5.10。
例5.10(主成分回归的原理) 设自变量x1,x2,x3和因变量y的观测值如表5-8,以x1,x2,x3为自变量y为因变量作回归。由于x3的值近似为x1与x2 的和,存在多重共线性,直接以x1,x2,x3为自变量不好,作主成分回归。
首先将x1,x2,x3零均值化,得到x1*,x2*,x3*
x1*=(x1-0.22)/0.49193,
x2*=(x2-0.24)/0.53666,
x3*=(x3-0.46)/0.50794
用相关阵计算主成分(计算所得表略去),主成分得分见表5-9,因为前一、两个主成分的累计百分比分别为0.576和0.992,取前两个主成分为自变量(略去1个主成分)。由表5-9数据作以y为因变量,Prin1和Prin2为自变量的回归方程
y=1.48+1.40128Prin1+0.31623Prin2
有主成分分析算得结果,Prin1,Prin2可用x1*,x2*,x3*表示:
prin1=0.461761x1*+0.461761x2*+0.757333x3*
prin2=0.707107x1*-0.707107x2*
通过变量替换,变成以x1,x2,x3为自变量的回归方程
y=-0.0598224+1.7699x1+0.789043x2+2.08929x3
这就是主成分回归所得方程。
表5-8 变量观测值
X1
x2
x3
y
1.1
0
1
4
0
1.2
1
3
0
0
0.3
0.3
0
0
0
0.1
0
0
0
0
表5-9 变量主成分得分(用协差阵计算)
Y
Prin1
Prin2
Prin3
4.0
0.55132
0.90525
0.01335
3.0
0.88824
-0.68739
0.01118
0.3
-0.33232
-0.04123
-0.13951
0.1
-0.55362
-0.08832
0.05749
0.0
-0.55362
-0.08832
0.05749
上述步骤很麻烦,可用SAS-REG过程实现:在proc reg语句中,加选项outest=文件名(存储主成分回归的输出数据集)在model语句自变量后加“/”号,再加选项pcomit=k(k为略去的主成分个数) 。执行后看最后一张表最后一行即可。
例5.11 1978-1993我国民航客运量y(万人),国民收入x1(亿元),消费额x2(亿元),铁路客运量x3(万人),民航航线里程数x4(万公里),来华旅游人数x5(万人)如表5-10。试建立以民航客运量为因变量的回归方程。
表5-10 我国民航数据表
year
y
x1
x2
x3
x4
x5
1978
231
3010
1888
81491
14.89
180.92
1979
298
3350
2195
86389
16.00
420.39
1980
343
3688
2531
92204
19.53
570.25
1981
401
3941
2799
95300
21.82
776.71
1982
445
4258
3054
99922
23.27
792.43
1983
391
4736
3358
106044
22.91
947.70
1984
554
5652
3905
110353
26.02
1285.22
1985
744
7020
4879
112110
27.72
1783.30
1986
997
7859
5552
108579
32.43
2281.95
1987
1310
9313
6386
112429
38.91
2690.23
1988
1442
11738
8038
122645
37.38
3169.48
1989
1283
13176
9005
113807
47.19
2450.14
1990
1660
14384
9663
95712
50.68
1746.20
1991
2178
16557
10969
95081
55.91
3335.65
1992
2886
20223
12985
99693
83.66
3311.50
1993
3383
24882
15949
105458
96.08
4152.70
解 直接建立回归方程发现,回归方程消费额x2的系数是负的,这与实际情况不符:消费额越多,应当乘飞机人数越多;而回归方程x2的系数是负的,导出错误结论:消费额越多,则乘飞机人数越少。通过检验发现,x2的系数是负的原因是:自变量出现多重共线性。本例中自变量有5个,主成分也有5个。通过主成分分析,前3个主成分累计百分比为99.61,可以选用前3个主成分,即去掉2个主成分,于是可用如下SAS程序计算
data airpline;
do year= 1978 to 1993;
input y x1-x5;
output;
end;
cards;
231 3010 1888 81491 14.89 180.92
298 3350 2195 86389 16.00 420.39
343 3688 2531 92204 19.53 570.25
401 3941 2799 95300 21.82 776.71
445 4258 3054 99922 23.27 792.43
391 4736 3358 106044 22.91 947.70
554 5652 3905 110353 26.02 1285.22
744 7020 4879 112110 27.72 1783.30
997 7859 5552 108579 32.43 2281.95
1310 9313 6386 112429 38.91 2690.23
1442 11738 8038 122645 37.38 3169.48
1283 13176 9005 113807 47.19 2450.14
1660 14384 9663 95712 50.68 1746.20
2178 16557 10969 95081 55.91 3335.65
2886 20223 12985 99693 83.66 3311.50
3383 24882 15949 105458 96.08 4152.70
;
proc reg;
model y=x1-x5/selection=stepwise;
model y=x1-x5/selection=stepwise;
run;
proc reg outest=w; /*主成分回归计算结果存数据集wu*/
model y=x1-x5/pcomit=2; /*实行主成分分析,去掉2个主成分*/
proc print data=w;
run;
执行后输出的最后一个表为
_ _ I
_ D _ P N
M _ E R C _ T
O T P I O R E
D Y V D M M R
O E P A G I S C
B L E R E T E E X X X X X
S _ _ _ _ _ _ P 1 2 3 4 5 Y
1 MODEL1 PARMS Y . . 75.9488 622.327 0.13433 -0.15715 -0.009738 18.4435 0.29278 -1
2 MODEL1 IPC Y . 2 96.1468 702.573 0.03878 0.05967 -0.010648 9.9057 0.21851 -1
这表倒数第2行和最后一行分别是:不减少主成分和减少两个主成分所得参数估计。由这表最后一行可见主成分方法得到的回归方程就是
y=702.573+0.03878x1+0.05967x2-.010648x3+9.9057 x4+0.21851x5
其中消费额x2的系数就是正的了。
5.5 经验正交函数分解
经验正交函数分解也称为经验正交分解或自然正交分解,其算法类似主成分分析,但含义不同,分析方法不同。主成分分析是对随机向量作分析,而经验正交函数分解是对确定性变量进行分析。由于经验正交函数分解在气象上的用途较大,在气象上称为EOF方法,见吴洪宝(2005)。
对气象要素场的分解,目前已有多种方法,如谐波分析、球函数分解、Chebixief分解。与它们相比,经验正交函数分解更有其优越性。(1)没有固定的函数形式,因而不需许多数学假设,更易符合实际。(2)能在不规则分布站点使用。(3)既可以在空间不同点进行分解,也可在同一站点对不同时间点分解,也可对同一站点不同要素做分解。
为了便于叙述经验正交分解,我们举不同站点,不同时间气温的例5.12和例5.13作为例子;容易看出气温可换为任何其他指标,站点这一属性变量也可换为别的属性变量,如例5.14中换为害虫的品种。
例5.12 设有p个站点,在i时刻j个站点温度观测值为
,温度矩阵为
,令
,
即
是各站点观测值减去对时间的平均值(距平),
形成矩阵
。
设样本协差阵为R,则
,若
为R的第j个彼此正交的单位特征向量,则
,
是第j个主成分得分所成向量。用分块矩阵写成
由特征向量正交性,通过分块矩阵乘法易证
,
于是,由分块矩阵乘法
。称
为第j个空间函数,
彼此正交。每一个
是一个空间正交
经验函数。和主成分分析一样,通常只取少数前k个
和
,
,即F可用少数空间正交经验
函数描述。
我们也可以考虑时间正交经验函数,令
, 设
形成矩阵
。
设
,若
为H的第j个彼此正交的单位特征向量,则
,
是第j个主成分得分所成向量。用分块矩阵写成
由特征向量正交性,通过分块矩阵乘法易证
,
于是,由分块矩阵乘法
。称
为第j个时间函数,
彼此正交。每一个
是
一个正交经验函数。和主成分分析一样,通常只取较小的数k,只考虑前k个
和
,
,
即G可用少数时间正交经验函数描述。
上述分解称为正交经验函数分解,也称为正交经验分解。
容易看出,正交经验函数分解与主成分分析的计算相同:若有q行p列数据矩阵F,将矩阵q个行变量作为p维随机变量的q次观测值,用协方差阵做主成分分析。各个主成分的系数构成p个正交经验函数,p个主成分得分就是正交经验函数的系数。
正交经验(函数)分解可以用SAS-PRINCOMP过程完成,例5.13说明其过程。
例5.13 设某三气象站(a,b,c)五次气温观测如表4-11,试做经验正交函数分解。
表5-11三气象站五次气温观测值
a
b
c
10
12
14
9
11
12
9
12
13
10
12
14
8
11
13
用空间经验函数分解,可采用下列程序
data temperat;
input a b c;
cards;
10 12 14
9 11 12
9 12 13
10 12 14
8 11 13
;
proc princomp cov out=w1;
var a b c;
proc print data=w1;
run;
执行该程序后,得到的主要输出为
Simple Statistics
A B C
Mean 9.200000000 11.60000000 13.20000000
StD 0.836660027 0.54772256 0.83666003
上表指出三地气温平均值分别是9.2,11.6,13.2;因而距平F应是观测值减去此三数,即{
}为
0.8
0.4
0.8
-0.2
-0.6
-1.2
-0.2
0.4
-0.2
0.8
0.4
0.8
-1.2
-0.6
-0.2
Eigenvectors
PRIN1 PRIN2 PRIN3
A 0.642542 0.707107 -.295194
B 0.417468 0.000000 0.908692
C 0.642542 -.707107 -.295194
上表给出三个特征向量,它们是彼此正交的经验函数(空间函数)。第一特征向量
中各分量相差不多,反映三地气温偏高的特征;第二特征向量
中第1,3个分量绝对值相同,符号相反,反映a地气温偏高,c地气温偏低的趋势。第三特征值很小,第三特征向量不考虑
OBS A B C PRIN1 PRIN2 PRIN3
1 10 12 14 1.19505 0.00000 -0.10883
2 9 11 12 -1.15004 0.70711 -0.13194
3 9 12 13 -0.09003 -0.00000 0.48155
4 10 12 14 1.19505 0.00000 -0.10883
5 8 11 13 -1.15004 -0.70711 -0.13194
上表给出
,
说明第一特征在第1,4年较强(1、4年气温偏高),第2,5年反方向较强(2、5年气温偏低),第三年很弱(气温接近平均值);第二特征第二年较强(a地气温距平比c地高),第五年反方向较强(c地气温比a地高)。
若想得到时间正交经验分解,则须把5年的观测值作为5个随机变量,a、b、c,3地代表不同观测次数,就可得另一分解。这时
采用程序
data temperat;
input year1-year5;
cards;
10 9 9 10 8
12 11 12 12 11
14 12 13 14 13
;
proc princomp cov out=w2;
var year1-year5;
proc print data=w2;
run;
执行程序后得到的主要结果是
Eigenvectors
PRIN1 PRIN2 PRIN3 PRIN4 PRIN5
YEAR1 0.436094 -.465405 -.037197 0.303046 0.707107
YEAR2 0.334245 0.240358 0.911322 0.000000 0.000000
YEAR3 0.450442 0.713419 -.353370 0.404061 0.000000
YEAR4 0.436094 -.465405 -.037197 0.303046 -.707107
YEAR5 0.552292 0.007656 -.204583 -.808122 0.000000
这5个特征向量就是彼此正交的时间经验函数。由上表可见第1特征向量是
反映五年气温偏高特征;第2特征向量是
,反映第三年偏高,第1,4年偏低特征。
OBS YEAR1 YEAR2 YEAR3 YEAR4 YEAR5 PRIN1 PRIN2 PRIN3 PRIN4 PRIN5
1 10 9 9 10 8 -4.82526 -0.22404 0 0 0
2 12 11 12 12 11 0.59581 0.55828 0 0 0
3 14 12 13 14 13 4.22945 -0.33425 0 0 0
上表表明
,中 -4.82546和4.22945绝对值相对很大符号一正一负,0.59581绝对值相对偏小,说明第一特征在a,c地突出,a地5年温度偏低,c地5年温度偏高;
。0.55828绝对值相对偏大,这说明第二特征在b地突出,b地第三年偏高,第1,4年偏低。
例5.14 已知山楂园昆虫群落调查表如表5-12,其中x1-x16表示不同时刻昆虫的群落x1:桃蚜,x2:山楂木虱,x3:草履蚧,x4:山楂叶螨, x5:梨网蝽,x6:黑绒金龟子,x7:苹毛金龟子, x8:顶梢叶蛾,x9:苹小卷叶蛾,x10:金纹细蛾,x11:舟形毛虫,x12:山楂粉蝶,x13:桃小食心虫,x14:梨小食心虫,x15:白小食心虫,x16:桑天蛾。T表示时段。做经验正交函数分解。
表5-12 山楂园昆虫群落调查表
t
x1
x2
x3
x4
x5
x6
x7
x8
x9
x10
x11
x12
x13
x14
x15
x16
1
0
10
4
0
0
9
0
0
0
0
0
0
0
0
9
0
2
24
18
7
0
0
18
7
0
0
0
0
13
0
0
0
0
3
329
53
13
0
8
274
182
7
1
0
0
22
0
0
5
0
4
675
86
2
11
43
419
619
12
1
0
0
34
0
0
7
0
5
266
123
28
7
47
64
253
31
14
4
0
46
0
0
23
0
6
38
205
35
16
94
13
47
64
17
9
0
31
0
16
32
0
7
2
180
17
34
125
4
0
23
44
35
0
7
0
72
13
0
8
0
71
6
53
207
0
0
11
17
60
0
0
0
115
11
1
9
0
23
0
89
391
0
0
7
8
125
0
0
0
9
143
4
10
0
10
0
74
647
0
0
13
11
153
10
0
37
26
19
6
11
0
0
0
93
561
0
0
1
4
65
126
0
72
289
8
15
12
0
0
0
64
174
0
0
0
6
70
284
0
346
21
3
5
13
0
0
0
23
13
0
0
0
3
213
133
0
295
93
23
0
14
0
0
0
8
0
0
0
0
1
145
13
0
82
41
15
0
15
0
0
0
3
0
0
0
0
0
43
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
解 设变量x1-x16表示各种昆虫群落,程序为
data insect;
input t x1-x16;
cards;
1 0 10 4 0 0 9 0 0 0 0 0 0 0 0 9 0
2 24 18 7 0 0 18 7 0 0 0 0 13 0 0 0 0
3 329 53 13 0 8 274 182 7 1 0 0 22 0 0 5 0
4 675 86 2 11 43 419 619 12 1 0 0 34 0 0 7 0
5 266 123 28 7 47 64 253 31 14 4 0 46 0 0 23 0
6 38 205 35 16 94 13 47 64 17 9 0 31 0 16 32 0
7 2 180 17 34 125 4 0 23 44 35 0 7 0 72 13 0
8 0 71 6 53 207 0 0 11 17 60 0 0 0 115 11 1
9 0 23 0 89 391 0 0 7 8 125 0 0 0 9 143 4
10 0 10 0 74 647 0 0 13 11 153 10 0 37 26 19 6
11 0 0 0 93 561 0 0 1 4 65 126 0 72 289 8 15
12 0 0 0 64 174 0 0 0 6 70 284 0 346 21 3 5
13 0 0 0 23 13 0 0 0 3 213 133 0 295 93 23 0
14 0 0 0 8 0 0 0 0 1 145 13 0 82 41 15 0
15 0 0 0 3 0 0 0 0 0 43 0 0 0 0 0 0
16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
;
proc princomp cov n=4 out=eo;
var x1-x16;
run;
proc print data=eo;
var prin1-prin4;
run;
执行后的主要输出为
Eigenvalues
Eigenvalue Difference Proportion Cumulative
PRIN1 84013.7 43593.6 0.538676 0.538676
PRIN2 40420.1 22524.4 0.259164 0.797841
PRIN3 17895.7 13220.9 0.114743 0.912584 PRIN4 4674.9 . 0.029974 0.942558
以上是协方差阵的前4个特阵值,累积比例0.9426表明只要取4个特阵向量-经验正交函数就够了。
Eigenvectors
PRIN1 PRIN2 PRIN3 PRIN4
X1 0.615750 0.268084 0.105590 -.024590
X2 0.075776 0.000319 -.192503 0.662614
X3 0.009459 -.008403 -.030454 0.098534
X4 -.065180 0.123261 0.042156 0.037380
X5 -.369040 0.888579 -.092851 -.056854
X6 0.379858 0.163405 0.080058 -.139139
X7 0.531713 0.249775 0.098918 0.059858
X8 0.009834 0.006608 -.046147 0.139279
X9 -.006596 0.005590 -.024813 0.099052
X10 -.120040 0.065880 0.218413 -.357796
X11 -.086355 0.034270 0.524309 0.237510
X12 0.039251 0.005635 -.018196 0.074936
X13 -.107747 -.044523 0.765976 0.047768
X14 -.108612 0.155797 0.098543 0.531365
X15 -.023971 0.046099 -.046134 -.138075
X16 -.006550 0.014258 0.005916 0.010357
以上是前4个特征向量,即4个主要正交经验函数。第1特征向量中x1、x7系数大,x5负系数大,表明桃蚜和苹毛金龟子高发而梨网蝽低发;第2特征向量中x5系数大,表明梨网蝽高发, 第3征向量中x11,x13系数大,表明舟形毛虫和桃小食心虫高发; 第4特征向量中x2,x14系数大,表明山楂木虱,梨小食心虫高发。
OBS PRIN1 PRIN2 PRIN3 PRIN4
1 -31.461 -188.305 -73.085 -29.684
2 -8.182 -179.016 -70.590 -23.294
3 369.967 -4.305 -8.985 -31.552
4 859.690 254.049 73.040 -13.554
5 279.934 -0.902 -44.508 51.884
6 -3.851 -75.585 -112.096 114.183
7 -79.377 -59.232 -104.294 112.322
8 -129.836 22.840 -79.173 44.987
9 -203.233 184.630 -87.230 -95.750
10 -304.774 407.810 -62.373 -98.047
11 -306.690 370.828 43.951 100.632
12 -176.478 -25.204 346.231 28.616
13 -121.268 -154.751 278.666 -18.080
14 -68.164 -175.742 33.614 -58.972
15 -40.817 -186.957 -61.825 -49.482
16 -35.460 -190.160 -71.343 -34.209
上表为主成分得分。Prin1得分中3,4,5时段得分数大,表明该三个时段桃蚜和苹毛金龟子高发高发;Prin2得分中4,9,10,11时段得分数大,表明上述时段梨网蝽高发(与实际情况略有差异,因有第1主成分影响);Prin3得分中第12,13时段得分大,表明这2时段舟形毛虫和桃小食心虫高发Prin4得分中6,7,10时短得分大,表明这三时段山楂木虱,梨小食心虫高发。
练习题
1. 1984洛杉矶奥运会发布国家女子径赛纪录如表5-13(100m,200m,400m赛成绩以秒为单位,其余一分钟为单位)。作如下统计分析:(1)求得各项纪录的样本相关阵。(2)写出前两个主成分,为了解释标准化变量的变化情况,这两个变量是否够了?(3)解释这两个主成分的含义。(4)按照第一主成分的大小将这25个国家排序,所排出的结果说明什麽?
表5-13 国家女子径赛纪录
Country
m100
m200
m400
m800
m1500
m3000
Marathon
Argentina
11.61
22.94
54.50
2.15
4.43
9.79
178.52
Australia
11.20
22.35
51.08
1.98
4.13
9.08
152.37
Austria
11.43
23.09
50.62
1.99
4.22
9.34
159.37
Belgium
11.41
23.04
52.00
2.00
4.14
8.88
157.85
Bermuda
11.46
23.05
53.30
2.16
4.58
9.81
169.98
Brazil
11.31
23.17
52.80
2.10
4.49
9.77
168.75
Burma
12.14
24.47
55.00
2.18
4.45
9.51
191.02
Canada
11.00
22.25
50.06
2.00
4.06
8.81
149.45
Chile
12.00
24.52
54.90
2.05
4.23
9.37
171.38
China
11.95
24.41
54.97
2.08
4.33
9.31
168.48
Columbia
11.60
24.00
53.26
2.11
4.35
9.46
165.42
Cook Isiand
12.90
27.10
60.40
2.30
4.84
11.10
233.22
Costa Rica
11.96
24.60
58.25
2.21
4.68
10.43
171.80
Czechoslovakia
11.09
21.97
47.99
1.89
4.14
8.92
158.85
Denmark
11.42
23.52
53.60
2.03
4.18
8.71
151.75
Dominicap
11.79
24.05
56.05
2.24
4.74
9.89
203.88
Finland
11.13
22.39
50.14
2.03
4.10
8.92
154.23
France
11.15
22.59
51.73
2.00
4.14
8.98
155.27
GDR
10.81
21.71
48.16
1.93
3.96
8.75
157.68
FRG
11.01
22.39
49.75
1.95
4.03
8.59
148.53
United Kingdom
11.00
22.13
50.46
1.98
4.03
8.62
149.72
Greece
11.79
24.08
54.93
2.07
4.35
9.87
182.20
Guatemala
11.84
24.54
56.09
2.28
4.86
10.54
215.08
Hungary
11.45
23.06
51.50
2.01
4.14
8.98
156.37
India
11.95
24.28
53.60
2.10
4.32
9.98
188.03
Indonesia
11.85
24.24
55.34
2.22
4.61
10.02
201.28
Ireland
11.43
23.51
53.24
2.05
4.11
8.89
149.38
Israel
11.45
23.57
54.90
2.10
4.25
9.37
160.48
Italy
11.29
23.00
52.01
1.96
3.98
8.63
151.82
Japan
11.73
24.00
53.73
2.09
4.35
9.20
150.50
Kenya
11.73
23.88
52.70
2.00
4.15
9.20
181.05
South Korea
11.96
24.49
55.70
2.15
4.42
9.62
164.65
North Korea
12.25
25.78
51.20
1.97
4.25
9.35
179.17
Luxembourg
12.03
24.96
56.10
2.07
4.38
9.64
174.68
Malaysia
12.23
24.21
55.09
2.19
4.69
10.46
182.17
Mauritius
11.76
25.08
58.10
2.27
4.79
10.90
261.13
Mexico
11.89
23.62
53.76
2.04
4.25
9.59
158.53
Netherlands
11.25
22.81
52.38
1.99
4.06
9.01
152.48
NewZealand
11.55
23.13
51.60
2.02
4.18
8.76
145.48
Norway
11.58
23.31
53.12
2.03
4.01
8.53
145.48
PNG
12.25
25.07
56.96
2.24
4.84
10.69
233.00
Philippines
11.76
23.54
54.60
2.19
4.60
10.16
200.37
Poland
11.13
22.21
49.29
1.95
3.99
8.97
160.82
Portugal
11.81
24.22
54.30
2.09
4.16
8.84
151.20
Rumania
11.44
23.46
51.20
1.92
3.96
8.53
165.45
Singapore
12.30
25.00
55.08
2.12
4.52
9.94
182.77
Spain
11.80
23.98
53.59
2.05
4.14
9.02
162.60
Sweden
11.16
22.82
51.79
2.02
4.12
8.84
154.48
Switzerland
11.45
23.31
53.11
2.02
4.07
8.77
153.42
Taipei
11.22
22.62
52.50
2.10
4.38
9.63
177.87
Thailand
11.75
24.46
55.80
2.20
4.72
10.28
168.45
Turkey
11.98
24.44
56.45
2.15
4.37
9.38
201.08
USA
10.79
21.83
50.62
1.96
3.95
8.50
142.72
USSR
11.06
22.19
49.19
1.89
3.87
8.45
151.22
Western Samoa
12.74
25.85
58.73
2.33
5.81
13.04
306.00
2. 1984洛杉矶奥运会发布国家男子径赛纪录如表5-14(100m,200m,400m赛成绩以秒为单位,其余以分钟为单位)。作如下统计分析:(1)求得各项纪录的样本相关阵。(2)写出前两个主成分,为了解释标准化变量的变化情况,这两个变量是否够了?(3)解释这两个主成分的含义。(4)按照第一主成分的大小将这25个国家排序,所排出的结果说明什麽?比较与女子成绩主成分的异同。
表5-14男子径赛纪录
Country
100m
200m
400m
800m
1500m
5Km
10Km
Marathon
Argentina
10.39
20.81
46.84
1.81
3.70
14.04
29.36
137.72
Australia
10.31
20.06
44.84
1.74
3.57
13.28
27.66
128.30
Austria
10.44
20.81
46.82
1.79
3.60
13.26
27.72
135.90
Belgium
10.34
20.68
45.04
1.73
3.60
13.22
27.45
129.95
Bermuda
10.28
20.58
45.91
1.80
3.75
14.68
30.55
146.62
Brazil
10.22
20.43
45.21
1.73
3.66
13.62
28.62
133.13
Burma
10.64
21.52
48.30
1.80
3.85
14.45
30.28
139.95
Canada
10.17
20.22
45.68
1.76
3.63
13.55
28.09
130.15
Chile
10.34
20.80
46.20
1.79
3.71
13.61
29.30
134.03
China
10.51
21.04
47.30
1.81
3.73
13.90
29.13
133.53
Columbia
10.43
21.05
46.10
1.82
3.74
13.49
27.88
131.35
Cook Islands
12.18
23.20
52.94
2.02
4.24
16.70
35.38
164.70
Costa Rica
10.94
21.90
48.66
1.87
3.84
14.03
28.81
136.58
Czechoslovakia
10.35
20.65
45.64
1.76
3.58
13.42
28.19
134.32
Denmark
10.56
20.52
45.89
1.78
3.61
13.50
28.11
130.78
Dominica
10.14
20.65
46.80
1.82
3.82
14.91
31.45
154.12
Finland
10.43
20.69
45.49
1.74
3.61
13.27
27.52
130.87
France
10.11
20.38
45.28
1.73
3.57
13.34
27.97
132.30
East Germany
10.12
20.33
44.87
1.73
3.56
13.17
27.42
129.92
West Germany
10.16
20.37
44.50
1.73
3.53
13.21
27.61
132.23
United Kingdom
10.11
20.21
44.93
1.70
3.51
13.01
27.51
129.13
Greece
10.22
20.71
46.56
1.78
3.64
14.59
28.45
134.60
Guatemala
10.98
21.82
48.40
1.89
3.80
14.16
30.11
139.33
Hungary
10.26
20.62
46.02
1.77
3.62
13.49
28.44
132.58
India
10.60
21.42
45.73
1.76
3.73
13.77
28.81
131.98
Indonesia
10.59
21.49
47.80
1.84
3.92
14.73
30.79
148.83
Ireland
10.61
20.96
46.30
1.79
3.56
13.32
27.81
132.35
Israel
10.71
21.00
47.80
1.77
3.72
13.66
28.93
137.55
Italy
10.01
19.72
45.26
1.73
3.60
13.23
27.52
131.08
Japan
10.34
20.81
45.86
1.79
3.64
13.41
27.72
128.63
Kenya
10.46
20.66
44.92
1.73
3.55
13.10
27.38
129.75
South Korea
10.34
20.89
46.90
1.79
3.77
13.96
29.23
136.25
North Korea
10.91
21.94
47.30
1.85
3.77
14.13
29.67
130.87
Luxembourg
10.35
20.77
47.40
1.82
3.67
13.64
29.08
141.27
Malaysia
10.40
20.92
46.30
1.82
3.80
14.64
31.01
154.10
Mauritius
11.19
22.45
47.70
1.88
3.83
15.06
31.77
152.23
Mexico
10.42
21.30
46.10
1.80
3.65
13.46
27.95
129.20
Netherlands
10.52
20.95
45.10
1.74
3.62
13.36
27.61
129.02
New Zealand
10.51
20.88
46.10
1.74
3.54
13.21
27.70
128.98
Norway
10.55
21.16
46.71
1.76
3.62
13.34
27.69
131.48
P.N.G.
10.96
21.78
47.90
1.90
4.01
14.72
31.36
148.22
Philippines
10.78
21.64
46.24
1.81
3.83
14.74
30.64
145.27
Poland
10.16
20.24
45.36
1.76
3.60
13.29
27.89
131.58
Portugal
10.53
21.17
46.70
1.79
3.62
13.13
27.38
128.65
Rumania
10.41
20.98
45.87
1.76
3.64
13.25
27.67
132.50
Singapore
10.38
21.28
47.40
1.88
3.89
15.11
31.32
157.77
Spain
10.42
20.77
45.98
1.76
3.55
13.31
27.73
131.57
Sweden
10.25
20.61
45.63
1.77
3.61
13.29
27.94
130.63
Switzerland
10.37
20.46
45.78
1.78
3.55
13.22
27.91
131.20
Taibei
10.59
21.29
46.80
1.79
3.77
14.07
30.07
139.27
Thailand
10.39
21.09
47.91
1.83
3.84
15.23
32.56
149.90
Turkey
10.71
21.43
47.60
1.79
3.67
13.56
28.58
131.50
USA
9.93
19.75
43.86
1.73
3.53
13.20
27.43
128.22
USSR
10.07
20.00
44.60
1.75
3.59
13.20
27.53
130.55
Western Samoa
10.82
21.86
49.00
2.02
4.24
16.28
34.71
161.83
3.为了测量显影液的质量,将每卷胶卷通过红、绿、蓝滤色片;并在高、中、低三种密度下进行测量。得到9个指标。由108卷胶卷所得数据产生样本协差阵如表5-15,试作主成分分析
表5-15
高红
高绿
高蓝
中红
中绿
中蓝
低红
低绿
低蓝
177
179
95
96
53
32
-7
-4
-3
179
419
245
131
181
127
-2
1
4
95
245
302
60
109
142
4
4
11
96
131
60
158
102
42
4
3
2
53
181
109
102
137
96
4
5
6
32
127
142
42
96
128
2
2
8
-7
-2
4
4
4
2
34
31
33
-4
1
4
3
5
2
31
39
89
-3
4
11
2
6
8
33
89
48
例1 利用武汉市2005年五月份的每天的各监测站平均的SO2监测值与每天早上八点钟的风力、气温、三小时降水作主分量分析。
解:采用以下程序
data temperat;
input date emo_so2
dd_08
t_08
r_08;
cards;
20050501
22.43
1
23.00
18.00
20050502
45.57
2
20.40
0.00
20050503
70.14
1
22.30
0.00
20050504
47.14
2
22.90
7.00
20050505
42.29
1
23.00
0.10
20050506
34.57
1
17.80
0.00
20050507
46.14
1
20.00
0.00
20050508
27.86
1
21.60
0.00
20050509
45.57
1
21.90
0.00
20050510
53.57
1
24.20
0.00
20050511
43.29
2
24.10
0.00
20050512
54.71
2
25.20
0.00
20050513
49.29
0
23.80
0.70
20050514
26.57
2
20.70
0.00
20050515
20.71
2
20.00
0.00
20050516
26.14
1
22.00
13.00
20050517
18.29
2
20.30
18.00
20050518
30.29
2
20.20
0.00
20050519
33.00
2
20.40
0.20
20050520
46.14
0
20.90
0.00
20050521
24.14
1
19.20
5.00
20050522
26.14
1
21.10
0.00
20050523
42.71
0
23.30
0.30
20050524
45.00
0
22.30
0.00
20050525
48.29
2
24.10
0.00
20050526
52.86
2
23.80
0.00
20050527
52.29
0
24.30
0.00
20050528
33.43
1
25.30
0.00
20050529
44.14
2
24.30
0.00
20050530
54.86
2
27.60
0.00
20050531
30.71
1
25.90
0.00
;
proc princomp cov out=prin;
var
dd_08
t_08
r_08;
proc print data=prin;
run;
执行上述程序,输出的结果是:
Eigenvalues of the Covariance Matrix
Eigenvalue Difference Proportion Cumulative
1 25.6674168 20.9350596 0.8299 0.8299
2 4.7323573 4.2041313 0.1530 0.9829
3 0.5282259 0.0171 1.0000
Eigenvectors
Prin1 Prin2 Prin3
dd_08 0.010759 0.003334 0.999937
t_08 -.074731 0.997201 -.002521
r_08 0.997146 0.074699 -.010978
上表给出所考察变量样本协差阵的特征值25.6674168、4.7323573、0.5282259和特征向量(0.010759,-0.074731,0.997146)、(0.003334, 0.997201, 0.074699)和(0.999937,-0.002521,-0.010978)。因此第一、二、三主分量分别是
prin1=0.010759*dd_08-0.074731*t_08+0.997146*r_08,
prin2=0.003334*dd_08+0.997201*t_08+0.074699*r_08,
prin3=0.999937*dd_08-0.002521*t_08-0.010978*r_08
由于第1个主分量累计贡献为0.8299,前两个第2个主分量所占比例为0.9829,所以对三个主分量只需要考虑前两个。
例2
资料:江苏省沿江七市1954—2002年梅雨量。资料说明:NJ-南京,SZ-苏州,WX-无锡,CZ-常州,ZJ-镇江,YZ-扬州,NT-南通。所用资料省略了小数位,6234表示623.4毫米。
方法:主分量分析。
程序主体如下:
DATA MEIYU;
TITLE 'MEIYU IN SELECTED CITIES';
INPUT YEAR NJ SZ WX CZ ZJ YZ NT;
CARDS;
1954 6234 3958 2799 3306 6441 6817 3500
1955 2083 2085 3528 3090 1321 913 4912
1956 3716 2254 2541 3753 3513 3650 2967
1957 2074 4360 4036 3610 2240 2728 3681
1958 612 887 826 281 919 639 1737
1959 1054 503 684 1293 1484 1746 969
1960 1060 901 1305 543 1622 1529 931
1961 1549 1971 2056 2226 846 593 1017
1962 3321 1449 3387 2303 1702 1771 2341
1963 1460 1624 2318 1548 1961 1548 1470
1964 1246 1320 1738 1803 727 588 1203
1965 562 223 142 379 1474 988 894
1966 1327 1456 1903 1380 1134 796 1602
1967 1280 1192 1822 2278 1793 1267 2287
1968 2085 1174 1687 2212 2943 2706 1756
1969 6073 1551 3039 4000 5984 5251 2796
1970 2515 2878 3940 4409 2317 2406 5934
1971 2290 1403 1281 1358 1371 1411 1019
1972 4097 787 2038 2102 6298 5582 2148
1973 1685 2477 1528 1314 1699 1385 1347
1974 2689 2775 3089 3481 3449 3034 2718
1975 4833 3880 4194 4630 5022 3936 4951
1976 2655 2575 2533 2091 1409 916 2089
1977 285 1156 312 426 392 604 227
1978 241 119 116 174 238 242 415
1979 2299 952 2003 2371 1672 2649 2619
1980 4582 2974 3053 3414 6908 6310 3303
1981 668 1276 1552 1398 1567 1793 2513
1982 4477 2981 3428 3280 3257 3899 3557
1983 3152 2648 2289 3374 2993 2208 2520
1984 2665 2347 2301 2583 2394 1277 2569
1985 1135 1446 2578 1719 1840 1142 2682
1986 2707 3296 2914 2874 3210 2408 4870
1987 3845 4079 4203 4504 3154 3546 3905
1988 935 917 596 658 825 1199 511
1989 2708 2478 2098 2027 2385 2438 1403
1990 1760 507 1136 1387 2316 1821 1553
1991 8799 5145 6777 8339 8217 7704 7349
1992 1597 530 664 697 1215 1162 739
1993 1992 1855 2052 1926 1501 1523 1937
1994 446 1100 1347 750 334 331 887
1995 2262 3774 2958 2493 2223 1428 2919
1996 5516 2374 4780 4002 4767 3570 3060
1997 1170 1107 848 878 1298 2006 884
1998 3034 1091 1118 1848 2739 2173 1501
1999 407 744 581 537 331 294 440
2000 556 344 661 1219 1542 1525 1552
2001 411 3031 2939 927 227 224 2334
2002 1710 1962 3282 2223 687 382 1905
;
proc princomp cov outstat=prin out=temp1;
VAR NJ SZ WX CZ ZJ YZ NT;
proc print data=temp1;
run;
结果分析:
第一主分量解释了总方差的80.7%,第二主分量解释了12.7%,二者共解释了总方差的93.3%。故只需分析第一和第二主分量就可概括江苏沿江地区梅雨的时空演变情况。
从第一主分量可以看出,第一特点是量值均为正数,表明江苏沿江地区梅雨量异常具有同步性,第二特点是宁、镇、扬地区为高值区,苏、锡、常、通为低值区。结合第一主分量时间系数曲线(图1)中存在的年际振荡的趋势,江苏沿江地区梅雨似乎存在一个振幅周期约3~4年的振荡。图中黑粗线为三阶多项式回归趋势线。可以看出,上世纪50年代末至70年代初为梅雨量较少的时期,70年代中后期至90年代初梅雨量偏多,1993年以来,梅雨量偏少,1999年以后似乎有回升的趋势。由于资料只到2002年,但我们知道2003年为我省历史罕见的梅雨量异常多年,因此也印证了这种变化趋势。
从第二主分量可以看出,宁、镇、扬地区为负值,苏、锡、常、通地区正值,表明这两个地区梅雨季节旱涝异常呈反位相。即当宁、镇、扬地区偏旱时,苏、锡、常、通地区偏涝;当宁、镇、扬地区偏涝时,苏、锡、常、通地区偏旱。也就是东旱西涝或东涝西旱。结合第二主分量时间系数曲线(图2)可以看出上世纪80年代末以前,以负系数为主(即东旱西涝),而90年代至今以正系数为主(既东涝西旱)。
综上可将江苏省沿江地区梅雨降水分成四种型:旱型、涝型、东旱西涝型、东涝西旱型。
例3 全国14个城市的春夏秋冬四季的平均季雷暴日如下表,试用主分量分析方法分析。
Spring Summer Autumn Winter
Shanghai 6.9 19.7 5.4 0.3
Ganyu 5.1 29.7 3.4 0.1
Xuzhou 5.5 26.1 2.9 0.1
Sheyang 5.7 27.8 5.3 0.2
Qingjiang 6.0 29.5 4.1 0.4
Xuyi 6.3 33.0 4.1 0.2
Dongtai 5.4 29.5 4.8 0.4
Gaoyou 5.7 30.1 5.8 0.4
Qidong 6.4 23.4 5.8 0.2
Nantong 6.4 25.7 5.3 0.3
Nanjing 7.4 25.9 3.2 0.4
Wujin 7.7 26.6 4.2 0.3
Liyang 9.1 31.7 5.4 0.5
Wuxian 10.0 29.6 7.2 0.4
Sas程序如下:
data temperat;
input city Spring Summer Autumn Winter;
cards;
Shanghai 6.9 19.7 5.4 0.3
Ganyu 5.1 29.7 3.4 0.1
Xuzhou 5.5 26.1 2.9 0.1
Sheyang 5.7 27.8 5.3 0.2
Qingjiang 6.0 29.5 4.1 0.4
Xuyi 6.3 33.0 4.1 0.2
Dongtai 5.4 29.5 4.8 0.4
Gaoyou 5.7 30.1 5.8 0.4
Qidong 6.4 23.4 5.8 0.2
Nantong 6.4 25.7 5.3 0.3
Nanjing 7.4 25.9 3.2 0.4
Wujin 7.7 26.6 4.2 0.3
Liyang 9.1 31.7 5.4 0.5
Wuxian 10.0 29.6 7.2 0.4
;
proc princomp data = temperat cov outstat = prin;
var Spring Summer Autumn Winter;
proc print data = prin;
run;
运行结果如下:
Eigenvalues of the Covariance Matrix
Eigenvalue Difference Proportion Cumulative
1 12.1230930 9.4749228 0.7772 0.7772
2 2.6481702 1.8309751 0.1698 0.9470
3 0.8171951 0.8077413 0.0524 0.9994
4 0.0094538 0.0006 1.0000
Eigenvectors
Prin1 Prin2 Prin3 Prin4
Spring 0.032046 0.818308 -.572537 -.039334
Summer 0.999396 -.020388 0.027273 -.006921
Autumn -.010823 0.572720 0.819409 -.021051
Winter 0.007957 0.044148 -.005087 0.998980
结果分析:
上表给出协差阵的四个特征值分别为:
λ1=12.1230930
λ2=2.6481702
λ3=0.8171951
λ4=0.0094538
对应得特征向量分别为:
(0.032046,0.999396,-.010823,0.007957)
(0.818308,-.020388,0.572720,0.044148)
(-.572537,0.027273,0.819409,-.005087)
(-.039334,-.006921,-.021051,0.998980)
由此可得:
Prin1= 0.032046Spring+ 0.818308Summer-.572537Autumn -.039334Winter
Prin2= 0.999396Spring -.020388+ Summer+ 0.027273Autumn -.006921Winter
Prin3= -.010823Spring+ 0.572720 Summer+ 0.819409Autumn -.021051Winter
Prin4= 0.007957Spring+ 0.044148 Summer -.005087Autumn +0.998980Winter
Spring 和Summer两个主分量解释了总方差的94.7%,可见只要用Spring 和Summer两个主分量就能很好的概括这组数据。
例4
某矿石中5种成分的百分比含量如下,资料的特点是每个样品的5个变量取值之和为100%,而且,每个变量的取值全为正数。试寻找少数几个综合指标反映原来5个指标85%以上的信息。
编号 X1 X2 X3 X4 X5 编号 X1 X2 X3 X4 X5
1 48.8 31.7 3.8 6.4 9.3 14 41.4 12.9 23.4 15.8 6.5
2 48.2 23.8 9.0 9.2 9.8 15 26.2 17.5 15.8 28.3 12.2
3 37.0 9.1 34.2 9.5 10.2 16 32.3 7.3 40.9 12.9 6.6
4 50.9 23.8 7.2 10.1 8.0 17 43.2 44.3 1.0 7.8 3.7
5 44.2 38.3 2.9 7.7 6.9 18 49.5 32.3 3.2 8.7 6.3
6 52.3 26.2 4.2 12.5 4.8 19 42.3 15.8 20.4 8.3 13.2
7 44.6 33.0 4.6 12.2 5.6 20 44.6 11.5 23.8 11.6 8.5
8 34.6 5.2 42.9 9.6 7.7 21 45.8 16.6 16.8 12.0 8.8
9 41.2 11.7 26.7 9.6 10.8 22 49.9 25.0 6.8 10.9 7.4
10 42.6 46.6 0.7 5.6 4.5 23 48.6 34.0 2.5 9.4 5.5
11 49.9 19.5 11.4 9.5 9.7 24 45.5 16.6 17.6 9.6 10.7
12 45.2 37.3 2.7 5.5 9.3 25 45.9 24.9 9.7 9.8 9.7
13 32.7 8.5 38.9 8.0 11.9
程序: DATA a1;
%LET m=5;
%LET n=25;
ARRAY a(&n);
ARRAY y(&m);
ARRAY x(&n,&m);
DO i=1 TO &n;
a(i)=0;
DO j=1 TO &m;
INPUT x(i,j) @@;
x(i,j)=log(x(i,j));
END;
DO j=1 TO &m;
a(i)=a(i)+x(i,j);
END;
a(i)=a(i)/&m;
DO j=1 TO &m;
x(i,j)=x(i,j)-a(i);
END;
DO j=1 TO &m;
y(j)=x(i,j);
END;
OUTPUT;
END;
CARDS;
48.8 31.7 3.8 6.4 9.3
48.2 23.8 9.0 9.2 9.8
37.0 9.1 34.2 9.5 10.2
50.9 23.8 7.2 10.1 8.0
44.2 38.3 2.9 7.7 6.9
52.3 26.2 4.2 12.5 4.8
44.6 33.0 4.6 12.2 5.6
34.6 5.2 42.9 9.6 7.7
41.2 11.7 26.7 9.6 10.8
42.6 46.6 0.7 5.6 4.5
49.9 19.5 11.4 9.5 9.7
45.2 37.3 2.7 5.5 9.3
32.7 8.5 38.9 8.0 11.9
41.4 12.9 23.4 15.8 6.5
26.2 17.5 15.8 28.3 12.2
32.3 7.3 40.9 12.9 6.6
43.2 44.3 1.0 7.8 3.7
49.5 32.3 3.2 8.7 6.3
42.3 15.8 20.4 8.3 13.2
44.6 11.5 23.8 11.6 8.5
45.8 16.6 16.8 12.0 8.8
49.9 25.0 6.8 10.9 7.4
48.6 34.0 2.5 9.4 5.5
45.5 16.6 17.6 9.6 10.7
45.9 24.9 9.7 9.8 9.7
;
DATA a2;
SET a1;
PROC PRINCOMP covariance PREFIX=z;
VAR y1-y&m;
RUN;
运行结果:
Eigenvalues of the Covariance Matrix
Eigenvalue Difference Proportion Cumulative
1 1.59570484 1.49729856 0.9197 0.9197
2 0.09840629 0.06163619 0.0567 0.9764
3 0.03677010 0.03260761 0.0212 0.9976
4 0.00416249 0.00416249 0.0024 1.000
5 0.00000000 0.0000 1.0000
Eigenvectors
z1 z2 z3 z4 z5
y1 -.205444 -.149494 0.740766 -.432100 0.447214
y2 -.589186 -.074714 -.059929 0.666098 0.447214
y3 0.780682 -.057476 0.167619 0.398918 0.447214
y4 -.016375 0.822418 -.251544 -.245124 0.447214
y5 0.030322 -.540734 -.596912 -.387792 0.447214
结果分析:
Z1主要受Y2、Y3的影响,而且,它们的符号相反。即当这种矿石的样品中X2的含量高时,X3的含量就低,反之亦然。
例5对徐州,丰县和沛县三个地区从1961年到2000年逐年降水资料进行主分量分析。单位:mm
程序如下:
data rain;
input year c1 c2 c3;
cards;
1961
727.8
586.8
591.5
1962
1013.1
1285.4
1049.6
1963
1213.4
855.9
1367.4
1964
969.5
1103.3
1109.9
1965
1633.9
1158.6
1144.5
1966
611.7
462.7
570.7
1967
756.3
842.7
729.9
1968
638.9
486.7
550.0
1969
969.6
779.6
744.1
1970
793.6
792.9
721.7
1971
887.6
801.3
1178.9
1972
873.0
673.0
671.4
1973
808.5
899.7
865.9
1974
1021.3
848.2
897.9
1975
1065.1
639.2
1450.0
1976
784.8
1335.4
630.2
1977
695.3
747.4
742.8
1978
2123.2
1610.3
628.0
1979
946.4
922.3
959.2
1980
776.5
1260.0
646.9
1981
612.2
418.6
492.4
1982
1017.3
558.9
546.3
1983
950.3
563.7
556.5
1984
952.3
818.1
657.3
1985
966.6
1098.3
992.1
1986
681.6
431.2
630.7
1987
620.5
694.5
686.6
1988
500.6
352.0
425.9
1989
692.7
541.6
632.6
1990
1088.8
803.1
833.4
1991
781.9
664.4
649.9
1992
888.6
875.5
899.2
1993
836.4
800.2
745.8
1994
821.1
1533.9
609.9
1995
825.3
651.3
719.7
1996
965.8
673.0
659.7
1997
794.3
427.9
591.0
1998
1128.9
725.6
819.9
1999
641.9
462.4
571.6
2000
979.6
630.4
865.4
;
proc princomp data=rain cov outstat=prin;
var c1 c2 c3;
proc print data=prin;
例6 为预报入霉期y(6月1日为1),考虑6个气象因子x1,x2,x3,x4,x5,x6.从过去22年的气象资料查得下表,试找出主要的预报因子和回归方程.
序号
X1
X2
X3
X4
X5
X6
Y
1
31
7
16
5
4
265
23
2
30
5
4
7
4
262
23
3
33
10
0
0
0
258
3
4
25
4
6
0
6
262
20
5
26
6
12
5
7
260
26
6
27
9
19
4
9
266
27
7
27
7
19
4
5
259
19
8
31
13
4
2
2
257
6
9
31
8
1
0
2
266
16
10
28
14
0
0
4
265
22
11
25
16
18
4
7
268
24
12
30
12
5
2
4
262
30
13
24
5
22
9
8
264
28
14
28
3
19
2
4
262
24
15
30
0
0
0
0
264
24
16
27
2
14
4
8
259
30
17
26
10
7
3
9
262
17
18
30
11
1
0
2
260
9
19
28
6
7
0
5
260
20
20
29
9
22
1
5
259
16
21
32
13
0
0
1
263
9
22
20
7
12
0
5
251
16
采用以下程序:
data temperat;
input num x1 x2 x3 x4 x5 x6 y;
cards;
1 31 7 16 5 4 265 23
2 30 5 4 7 4 262 23
3 33 10 0 0 0 258 3
4 25 4 6 0 6 262 20
5 26 6 12 5 7 260 26
6 27 9 19 4 9 266 27
7 27 7 19 4 5 259 19
8 31 13 4 2 2 257 6
9 31 8 1 0 2 266 16
10 28 14 0 0 4 265 22
11 25 16 18 4 7 268 24
12 30 12 5 2 4 262 30
13 24 5 22 9 8 264 28
14 28 3 19 2 4 262 24
15 30 0 0 0 0 264 24
16 27 2 14 4 8 259 30
17 26 10 7 3 9 262 17
18 30 11 1 0 2 260 9
19 28 6 7 0 5 260 20
20 29 9 22 1 5 259 16
21 32 13 0 0 1 263 9
22 20 7 12 0 5 251 16
;
proc princomp data=temperat cov outstat=prin;
var x1 x2 x3 x4 x5 x6;
proc print data=prin;
run;
运行结果为:
Eigenvalues of the Covariance Matrix
Eigenvalue Difference Proportion Cumulative
1 73.6070177 55.4521548 0.6238 0.6238
2 18.1548629 4.8083404 0.1538 0.7776
3 13.3465225 5.7223905 0.1131 0.8907
4 7.6241320 3.7131947 0.0646 0.9553
5 3.9109373 2.5500806 0.0331 0.9885
6 1.3608567 0.0115 1.0000
Eigenvectors
Prin1 Prin2 Prin3 Prin4 Prin5 Prin6
x1 -.207877 0.180535 0.248947 0.724525 0.257144 0.520730
x2 -.122493 0.757868 -.632378 0.004947 0.085417 -.058390
x3 0.922915 0.098946 -.075883 0.332226 -.148182 -.018669
x4 0.190738 0.089625 0.218956 -.128328 0.888245 -.319683
x5 0.231309 0.065301 -.022736 -.556459 0.149949 0.780758
x6 0.011238 0.609062 0.695634 -.196343 -.305219 -.115330
由上面累计特征量(cumulative)列可看出:只要用第一、二、三个主分量即用:
prin1=-.207877x1+0.180535x2+0.248947x3+0.724525x4+0.257144x5+ 0.520730x6
prin2=-.122493x1+0.757868x2-.632378x3+0.004947x4+0.085417x5 -.058390x6
prin3=0.922915x1+0.098946x2-.075883x3+0.332226x4-.148182x5 -.018669x6
就能很好地概括这组数据了,这三个主分量解释了总方差的89.07℅。
例7 回归分析
一 资料选取及分析方法简介
资料选取的是1951---1980年(计29年)的气象资料。该回归方程的目的是研究长江中下游六月份降水量Y与六个预报因子X1---X6之间的相关关系。其中,X1:12月份北半球极涡强度。 X2:9月份40°N,90°W、100°W两点500hPa高度平均。 X3:10月份65°N,150°E-170°W五点500hPa高度平均。 X4:12月份180°E-160°W三点500hPa高度平均。X5:2月份80°N、150°W-130°W;70°N、150°W-130°W六点500hPa高度平均。X6:9月份月平均温度减20度。以上因子X1—X6中凡是高度平均均以减去5000位势米。(注:该资料中的数据已标准化)
分析方法是正交组合的多元回归。该方法是将M个因子正交组合成M个新的因子(即相互独立的主成份),然后利用正交组合因子建立多元回归方程。具体做法:用SAS的PROC PRINCOMP过程来分析组合因子,这里组合因子为主分量分析中的主成份,用主成份来建立回归的好处:一、由于组合因子彼此正交,建立回归方程时特别简单,回归系数可以独立计算。而且由于信息的重新组合分配,使预报因子减少成为可能。二、因为组合因子和预报对象的相关系数与原变量与预报对象的相关系数有很大的不同。三、组合因子能利用到各个原始变量的信息,这与逐步回归有很大不同。把PROC PRINCOMP得到的主成份与预报对象组成新的数据来用正交回归PROC ORTHOREG过程建立回归方程,同时根据组合因子与预报对象Y的相关系数来对预报因子进行取舍。
二 程序
(一)组合因子的产生-----PROC PRINCOMP过程
1 程序PROC PRINCOMP
DATA chang;
TITLE '长江中下游6月份降水量的Y与X1,X2,X3,X4,X5,X6的关系之组合因子分析';
INPUT x1 x2 x3 x4 x5 x6;
cards;
-0.362 -0.684 -0.672 1.305 0.529 -0.548
1.834 1.004 2.328 -1.685 -0.440 0.284
-2.069 0.160 1.256 -1.570 -2.217 -0.363
-1.094 1.848 2.113 0.385 -0.186 1.208
0.858 1.004 -0.244 1.880 -1.571 2.040
-0.118 0.582 -1.101 0.270 0.529 -0.271
0.126 -1.106 -0.030 -0.420 1.337 -1.842
1.346 0.582 -1.101 -0.075 -1.086 -0.641
-0.362 -1.106 1.256 -0.650 2.306 0.191
-0.606 1.004 -0.244 -0.535 -1.086 0.746
0.126 -1.106 0.185 0.615 1.337 0.284
0.370 -1.106 -0.030 0.615 0.852 0.284
-0.362 2.270 -0.672 -0.880 -1.247 0.653
0.614 0.160 -0.458 1.535 0.368 1.670
0.370 -1.106 -0.887 0.960 -0.117 -0.918
0.126 -1.106 -0.887 -0.190 -0.440 -1.295
2.078 -0.684 -1.530 0.615 -0.278 0.099
-0.118 -1.528 -0.244 1.075 -0.117 0.469
1.346 0.160 0.613 -1.685 1.660 1.116
0.614 0.160 0.185 0.500 -0.601 -0.271
-1.582 -0.262 -1.315 1.075 -0.117 -0.548
-0.362 0.160 -0.672 0.615 0.045 -1.288
-0.362 0.160 -0.458 -1.110 -0.601 -1.658
-1.094 -1.106 1.042 -1.570 0.206 -0.179
-1.826 -1.106 -1.315 -0.650 -0.278 2.040
1.346 -0.262 0.613 -1.340 -0.117 -1.011
0.126 0.582 1.256 0.615 1.014 -0.086
-0.362 1.426 0.185 0.155 0.206 0.653
-0.606 1.004 0.828 0.155 1.014 -0.918
;
PROC PRINCOMP data=chang outstat=prin;
var x1 x2 x3 x4 x5 x6;
proc print data=prin;
run;
2 PROC PRINCOMP过程产生的新的组合因子
-1.4925 -0.6176 -0.2976 0.3746 0.6870 0.1539 79.5
2.1892 1.8433 1.6778 -1.1641 -0.3773 0.6343 293.7
2.2939 0.7754 -2.4981 0.2298 -0.4534 1.0176 389.2
2.9734 -0.3686 0.1947 1.1060 1.0423 0.7352 298.9
1.1684 -2.8455 1.3777 0.0285 -0.0918 -0.7724 274.3
-0.4992 -0.5474 -0.1514 -.2140 0.4159 -1.0957 172.7
-1.7837 1.6799 -0.4490 -0.4052 0.4572 -.1380 55.5
0.0030 -0.8312 0.0536 -2.0415 -0.2565 -0.1761 221.8
-0.8553 2.1460 0.7510 1.6607 -0.0109 -0.2100 211.9
1.5047 -0.7266 -0.5684 0.0152 -0.4221 -0.4042 191.9
-1.3259 0.3784 0.7588 1.0034 0.0242 0.1837 156.5
-1.2378 0.0701 0.6934 0.6076 -0.2017 0.3430 86.3
2.2302 -1.0085 -.5056 -0.7373 -.1537 -1.2968 305.5
-0.2875 -1.6629 1.5455 0.8072 -0.0486 -0.0495 137.4
-1.6780 -0.5551 -0.3623 -0.5675 0.1180 0.5770 137.5
-1.2263 0.1116 -1.0241 -0.9073 -0.4390 0.3318 130.4
-1.4191 -1.1684 1.0846 -1.4798 -0.9161 0.0787 71.6
-1.1812 -0.7249 0.0285 0.7725 -0.4592 1.0031 139.4
0.3357 1.6312 1.9445 0.1894 -1.0951 -1.1240 197.2
0.1120 -0.3702 0.2667 -0.6690 0.2933 0.6062 218.0
-0.9717 -1.1936 -1.5904 0.6061 0.6363 -0.3143 182.1
-0.7718 -0.2807 -0.8398 -0.5374 0.9663 -0.2113 159.8
0.0324 0.7587 -1.5532 -1.2658 0.0750 -0.3341 155.3
0.2575 1.8429 -0.9976 0.8823 -0.8008 0.3803 268.4
-0.0028 -1.0067 -1.1570 1.9727 -2.0726 -0.6394 243.3
0.0548 1.5514 0.3334 -1.5011 -0.4695 0.3055 143.4
0.2139 0.6385 0.9757 0.5581 1.2699 0.1265 115.1
1.0972 -0.4269 0.3775 0.4217 0.6268 -0.7700 238.2
0.2726 0.9013 -0.0723 0.2844 1.6478 -0.5158 219.1
(二)组合因子建立回归方程----PROC ORTHOREG
data preci;
input z1-z6 y;
cards;
-1.4925 -0.6176 -0.2976 0.3746 0.6870 0.1539 79.5
2.1892 1.8433 1.6778 -1.1641 -0.3773 0.6343 293.7
2.2939 0.7754 -2.4981 0.2298 -0.4534 1.0176 389.2
2.9734 -0.3686 0.1947 1.1060 1.0423 0.7352 298.9
1.1684 -2.8455 1.3777 0.0285 -0.0918 -0.7724 274.3
-0.4992 -0.5474 -0.1514 -.2140 0.4159 -1.0957 172.7
-1.7837 1.6799 -0.4490 -0.4052 0.4572 -.1380 55.5
0.0030 -0.8312 0.0536 -2.0415 -0.2565 -0.1761 221.8
-0.8553 2.1460 0.7510 1.6607 -0.0109 -0.2100 211.9
1.5047 -0.7266 -0.5684 0.0152 -0.4221 -0.4042 191.9
-1.3259 0.3784 0.7588 1.0034 0.0242 0.1837 156.5
-1.2378 0.0701 0.6934 0.6076 -0.2017 0.3430 86.3
2.2302 -1.0085 -.5056 -0.7373 -.1537 -1.2968 305.5
-0.2875 -1.6629 1.5455 0.8072 -0.0486 -0.0495 137.4
-1.6780 -0.5551 -0.3623 -0.5675 0.1180 0.5770 137.5
-1.2263 0.1116 -1.0241 -0.9073 -0.4390 0.3318 130.4
-1.4191 -1.1684 1.0846 -1.4798 -0.9161 0.0787 71.6
-1.1812 -0.7249 0.0285 0.7725 -0.4592 1.0031 139.4
0.3357 1.6312 1.9445 0.1894 -1.0951 -1.1240 197.2
0.1120 -0.3702 0.2667 -0.6690 0.2933 0.6062 218.0
-0.9717 -1.1936 -1.5904 0.6061 0.6363 -0.3143 182.1
-0.7718 -0.2807 -0.8398 -0.5374 0.9663 -0.2113 159.8
0.0324 0.7587 -1.5532 -1.2658 0.0750 -0.3341 155.3
0.2575 1.8429 -0.9976 0.8823 -0.8008 0.3803 268.4
-0.0028 -1.0067 -1.1570 1.9727 -2.0726 -0.6394 243.3
0.0548 1.5514 0.3334 -1.5011 -0.4695 0.3055 143.4
0.2139 0.6385 0.9757 0.5581 1.2699 0.1265 115.1
1.0972 -0.4269 0.3775 0.4217 0.6268 -0.7700 238.2
0.2726 0.9013 -0.0723 0.2844 1.6478 -0.5158 219.1
;
proc orthoreg data=preci;
model y=z1-z6;
run;
三 回归结果及分析
(1) PROC PRINCOMP过程结果及分析:
1结果
长江中下游6月份降水量的Y与X1,X2,X3,X4,X5,X6的关系之组合因子分析 13
20:41 Sunday, May 3, 1998
Principal Component Analysis
29 Observations
6 Variables
Simple Statistics
X1 X2 X3
Mean -0.000172414 -0.000068966 0.000000000
StD 1.017775307 1.017647635 1.017603979
X4 X5 X6
Mean 0.000344828 0.0311724138 -0.003448276
StD 1.017639014 0.9970482532 1.022061348
长江中下游6月份降水量的Y与X1,X2,X3,X4,X5,X6的关系之组合因子分析 14
20:41 Sunday, May 3, 1998
Principal Component Analysis
Correlation Matrix
X1 X2 X3 X4 X5 X6
X1 1.0000 0.0118 -.0270 0.0636 0.1135 -.0144
X2 0.0118 1.0000 0.2637 -.0867 -.3734 0.2639
X3 -.0270 0.2637 1.0000 -.4187 0.1888 0.1039
X4 0.0636 -.0867 -.4187 1.0000 0.0568 0.1658
X5 0.1135 -.3734 0.1888 0.0568 1.0000 -.0929
X6 -.0144 0.2639 0.1039 0.1658 -.0929 1.0000
长江中下游6月份降水量的Y与X1,X2,X3,X4,X5,X6的关系之组合因子分析 15
20:41 Sunday, May 3, 1998
Principal Component Analysis
Eigenvalues of the Correlation Matrix
Eigenvalue Difference Proportion Cumulative
PRIN1 1.60037 0.191548 0.266728 0.26673
PRIN2 1.40882 0.307611 0.234804 0.50153
PRIN3 1.10121 0.148079 0.183535 0.68507
PRIN4 0.95313 0.373649 0.158855 0.84392
PRIN5 0.57948 0.222502 0.096581 0.94050
PRIN6 0.35698 . 0.059497 1.00000
长江中下游6月份降水量的Y与X1,X2,X3,X4,X5,X6的关系之组合因子分析 16
20:41 Sunday, May 3, 1998
Principal Component Analysis
Eigenvectors
PRIN1 PRIN2 PRIN3 PRIN4 PRIN5 PRIN6
X1 -.122551 -.031311 0.585940 -.779020 -.164388 0.082339
X2 0.612744 0.278584 0.101450 -.193306 0.538883 -.457036
X3 0.500793 -.492572 0.262576 0.158830 0.208034 0.607560
X4 -.399154 0.526026 0.302465 0.199008 0.528608 0.391735
X5 -.330872 -.504786 0.470005 0.322640 0.236233 -.504883
X6 0.299887 0.383758 0.514827 0.432242 -.551426 -.082743
长江中下游6月份降水量的Y与X1,X2,X3,X4,X5,X6的关系之组合因子分析 17
20:41 Sunday, May 3, 1998
OBS _TYPE_ _NAME_ X1 X2 X3 X4 X5 X6
1 MEAN -0.0002 -0.0001 0.0000 0.0003 0.0312 -0.0034
2 STD 1.0178 1.0176 1.0176 1.0176 0.9970 1.0221
3 N 29.0000 29.0000 29.0000 29.0000 29.0000 29.0000
4 CORR X1 1.0000 0.0118 -0.0270 0.0636 0.1135 -0.0144
5 CORR X2 0.0118 1.0000 0.2637 -0.0867 -0.3734 0.2639
6 CORR X3 -0.0270 0.2637 1.0000 -0.4187 0.1888 0.1039
7 CORR X4 0.0636 -0.0867 -0.4187 1.0000 0.0568 0.1658
8 CORR X5 0.1135 -0.3734 0.1888 0.0568 1.0000 -0.0929
9 CORR X6 -0.0144 0.2639 0.1039 0.1658 -0.0929 1.0000
10 EIGENVAL 1.6004 1.4088 1.1012 0.9531 0.5795 0.3570
11 SCORE PRIN1 -0.1226 0.6127 0.5008 -0.3992 -0.3309 0.2999
12 SCORE PRIN2 -0.0313 0.2786 -0.4926 0.5260 -0.5048 0.3838
13 SCORE PRIN3 0.5859 0.1015 0.2626 0.3025 0.4700 0.5148
14 SCORE PRIN4 -0.7790 -0.1933 0.1588 0.1990 0.3226 0.4322
15 SCORE PRIN5 -0.1644 0.5389 0.2080 0.5286 0.2362 -0.5514
16 SCORE PRIN6 0.0823 -0.4570 0.6076 0.3917 -0.5049 -0.0827
2结果分析
在表EIGENVALUES OF THE COVARIANCE MATRIX中给出协方差阵的六个特征值分别为:λ1=1.6004 , λ2= 1.4088,λ3=1.1012,λ4=0.9531,λ5= 0.5795,λ6= 0.3570,
对应的特征向量分别为:
t1 =(-0.1226 , 0.6127 , 0.5008 ,-0.3992 , -0.3309 , 0.2999)
t2 =(-0.0313 , 0.2786 ,-0.4926 , 0.5260 , -0.5048 , 0.3838)
t3 =(0.5859 ,0.1015 , 0.2626 , 0.3025 , 0.4700 ,0.5148)
t4 =(-0.7790 ,-0.1933 ,0.1588 ,0.1990 , 0.3226 ,0.4322 )
t5 =(-0.1644 ,0.5389 ,0.2080 ,0.5286 ,0.2362 ,-0.5514 )
t6 =(0.0823 ,-0.4570 , 0.6076 , 0.3917 ,-0.5049 ,-0.0827)
由此可以得到第一主分量
PRIN 1= -0.1226X1+ 0.6127X2+0.5008X3-0.3992X4-0.3309X5+0.2999X6
第二、三个等主分量类似可得。再将观测值代入主分量即得到新的组合因子(组合因子即主成份)。
(2) PROC ORTHOREG过程结果:
1结果
The SAS System 21:49 Sunday, May 3, 1998 4
ORTHOREG Regression Procedure
Dependent Variable Y
Sum of Squared Errors 38506.064031
Degrees of Freedom 22
Mean Squared Error 1750.2756378
Root Mean Sqr Error 41.836295699
R-square 0.7803532718
Variable DF Parameter Estimate Std Error T-Ratio Prob>|t|
INTERCEP 1 189.405371732589 7.8019879614 24.28 0.0001
z1 1 50.6878249217834 6.0551134365 8.37 0.0001
z2 1 1.37407801045017 6.704907377 0.20 0.8395
z3 1 -14.6819367504118 7.5323693735 -1.95 0.0641
z4 1 11.2190820462235 8.1047311281 1.38 0.1802
z5 1 -13.7989832122319 10.340164965 -1.33 0.1957
z6 1 -0.21598899881973 13.242646135 -0.02 0.9871
2 SAS结果分析
由上述程序的结果我们可以看出正交回归统计分析ORTHOREG过程给出两个表:方差分析表、参数估计表。
方差分析表 :残差平方和 SE = 38506.064031 也给出自由度为22,则平均平方和 MSE=SE/22= 1750.2756378 。表中还给出复相关系数平方R2 =0.7803532718。
参数估计表 :它给出常数项和回归系数分别是189.405371732589 50.6878249217834 1.37407801045017 -14.6819367504118 11.2190820462235 -13.7989832122319 -0.21598899881973;它们的标准误差是7.8019879614 6.0551134365 6.704907377 7.5323693735 8.1047311281 10.340164965 13.242646135 。t值为:24.28 8.37 0.20 -1.95 1.38 -1.33 -0.02并给出t分布随机变量大于它们的概率是0.0001 0.0001 0.8395 0.0641 0.1802 0.1957 0.9871。前两个概率小于0.05,说明截距和Z1的作用明显,不能舍去。
由参数估计表的参数估计值结果可以得到回归方程:
Y=189.41+ 50.69Z1+1.38Z2 -14.69Z3+11.22Z4 -13.80Z5 -0.22Z6
3 物理意义分析
原预报因子与预报对象的相关系数Riy为-.3641 0.5492 0.4458 -.4197 -.5036 0.3761(Riy的求法:将原预报因子X1-X6与预报对象Y建立一个数据集从而可用PROC PRINCOMP过程来计算,该过程的其中一个结果为各个因子之间的相关系数,这样可从中提取出X1-X6与Y的相关系数。)可见,预报因子与预报对象的关系为:X1: 12月份北半球极涡强度,X4:12月份180°E-160°W三点500hPa高度平均及X5:2月份80°N、150°W-130°W;70°N、150°W-130°W六点500hPa高度平均这三个因子与预报对象Y是反相关。而X2:9月份40°N,90°W、100°W两点500hPa高度平均, X3:10月份65°N,150°E-170°W五点500hPa高度平均,X6:9月份月平均温度减20度这三个因子与预报对象正相关。组合因子与预报对象的相关系数R*iy为 0.8394 0.0203 -.1961 0.1393 -.1332 -.0422(R*iy 的求法与Riy的求法相同,不同的是原预报因子变为组合因子Z1-Z6来建立数据集)组合因子Z1、Z2、Z4与预报因子呈正相关,组合因子Z3、Z5及Z6与预报对象呈负相关。(注:由于篇幅有限这里不打印出这两个相应的相关矩阵)。
由上述的结果可以看出:原预报因子与预报对象的关系都非常密切,相关系数相差不多,这是因为原预报因子之间有正交、彼此不相互独立。而组合因子与预报对象的关系则大不相同,第一个组合因子(即第一主分量)明显与预报因子的相关系数达到0.839,所以用组合因子,只要选取一个组合因子就可以预报六月份长江中下游的降水。可见用组合因子建立回归方程,预报方程的显著性大大提高了。
根据参数估计的t检验及气象中应用的相关系数R*iy可以看出组合因子中的Z1对预报方程的贡献最大,故方程可简化为:
Y=189.41+ 50.69Z1
根据预报对象Y及预报因子Z1画图,则图象如下:
The SAS System 11
21:31 Wednesday, May 6, 1998
Plot of Y*Z1. Symbol is value of Y.
Y
86.3 ? 8
79.5 ? 7
71.6 ? 7
55.5 ? 5
389.2 ? 3
305.5 ? 3
298.9 ? 2
293.7 ? 2
274.3 ? 2
268.4 ? 2
243.3 ? 2
238.2 ? 2
221.8 ? 2
219.1 ? 2
218.0 ? 2
211.9 ? 2
197.2 ? 1
191.9 ? 1
182.1 ? 1
172.7 ? 1
159.8 ? 1
156.5 ? 1
155.3 ? 1
143.4 ? 1
139.4 ? 1
137.5 ? 1
137.4 ? 1
130.4 ? 1
115.1 ? 1
妶儍儍儍儍儍儍儓儍儍儍儍儍儍儓儍儍儍儍儍儍儓儍儍儍儍儍儍儓儍儍儍儍儍儍儓
-2 -1 0 1 2 3
Z1
如上图,Y与Z1近似呈线性且斜率为正值。故可以看出用一个组合因子Z1就可以很好的建立这个线性回归方程。
四 总结
回归分析是气象研究和预报中经常使用的相关方法。我们可以用回归方法来建立相关预报方程(降水、气温、体感温度、湿度等预报方程)。方法可采用逐步回归法、选择最优的变量子集法、正交组合的多元回归等,本文采用正交组合的多元回归。该方法简化了方程、提高了预报的显著性。
本文采用了1951---1980年(计29年)的气象资料(五个高度长资料一个温度长资料) 及正交组合的多元回归方法调用SAS中的PROC PRINCOMP 和PROC ORTHOREG来建立相应的回归方程。从而经分析简化最终得到方程:
Y=189.41+ 50.69Z1,
其中,Z1用原始变量表示为:Z1=-0.1226X1+0.6127X2+0.5008X3-0.3992X4-0.3309X5+0.2999X6
� EMBED Excel.Chart.8 \s ���
� EMBED Excel.Chart.8 \s ���
_1219301885.unknown
_1233767390.unknown
_1238088899.unknown
_1238090798.unknown
_1240159317.unknown
_1288936923.unknown
_1320497401.xls
Chart1
-2661
3046
-429
2529
-137
-1223
-990
880
645
-10
575
-1410
478
396
-1192
-2691
3156
-610
-4125
11
298
1013
1166
-575
-803
-180
-2678
157
226
275
733
898
1584
1759
-923
-475
-1268
880
-1202
389
334
1764
-160
-1123
-1517
-354
-897
2525
1918
图2 江苏省沿江地区梅雨降水第二主分量的时间系数
Sheet1
1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002
-2661 3046 -429 2529 -137 -1223 -990 880 645 -10 575 -1410 478 396 -1192 -2691 3156 -610 -4125 11 298 1013 1166 -575 -803 -180 -2678 157 226 275 733 898 1584 1759 -923 -475 -1268 880 -1202 389 334 1764 -160 -1123 -1517 -354 -897 2525 1918
_1320497407.xls
Chart1
7076
491
2745
2168
-3701
-2781
-2818
-2200
247
-1433
-2745
-3933
-2377
-1420
-162
5648
2931
-1973
3614
-1683
2093
5950
-705
-4645
-5231
-281
6137
-1896
3543
1365
123
-1306
2285
4135
-3687
5
-1679
14102
-3205
-1139
-4072
551
4939
-2692
-464
-4654
-2955
-2686
-1625
图1 江苏省沿江地区梅雨降水第一主分量的时间系数
Sheet1
1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002
7076 491 2745 2168 -3701 -2781 -2818 -2200 247 -1433 -2745 -3933 -2377 -1420 -162 5648 2931 -1973 3614 -1683 2093 5950 -705 -4645 -5231 -281 6137 -1896 3543 1365 123 -1306 2285 4135 -3687 5 -1679 14102 -3205 -1139 -4072 551 4939 -2692 -464 -4654 -2955 -2686 -1625
_1240159411.unknown
_1240159294.unknown
_1238090949.unknown
_1238089965.unknown
_1238090137.unknown
_1238090309.unknown
_1238089995.unknown
_1238090056.unknown
_1238090071.unknown
_1238090017.unknown
_1238089983.unknown
_1238089717.unknown
_1238089899.unknown
_1238089403.unknown
_1238089492.unknown
_1238089025.unknown
_1238088221.unknown
_1238088771.unknown
_1238088807.unknown
_1238088300.unknown
_1238087971.unknown
_1238088057.unknown
_1238086951.unknown
_1238087467.unknown
_1233891821.unknown
_1220114499.unknown
_1233767152.unknown
_1233767337.unknown
_1233767352.unknown
_1233767168.unknown
_1233767179.unknown
_1233767213.unknown
_1233767073.unknown
_1233767136.unknown
_1233767095.unknown
_1221655798.unknown
_1223714519.unknown
_1233422632.unknown
_1233759817.unknown
_1233422557.unknown
_1221661032.unknown
_1220114524.unknown
_1219323325.unknown
_1219323899.unknown
_1220114338.unknown
_1220114365.unknown
_1220114389.unknown
_1219339133.unknown
_1220114294.unknown
_1219339193.unknown
_1219339022.unknown
_1219323397.unknown
_1219323479.unknown
_1219323362.unknown
_1219301992.unknown
_1219322115.unknown
_1219322222.unknown
_1219303393.unknown
_1219303442.unknown
_1219302602.unknown
_1219301951.unknown
_1219301973.unknown
_1219301926.unknown
_951112282.unknown
_951317898.unknown
_1172250695.unknown
_1172251142.unknown
_1176539915.unknown
_1180543837.unknown
_1219301663.unknown
_1219301826.unknown
_1194624731.unknown
_1176540098.unknown
_1180543630.unknown
_1176540040.unknown
_1172251342.unknown
_1173249968.unknown
_1176178923.unknown
_1176436564.unknown
_1172251369.unknown
_1172251088.unknown
_1172251105.unknown
_1172250900.unknown
_1172250955.unknown
_1172250712.unknown
_1172250723.unknown
_1172250419.unknown
_1172250594.unknown
_1172250607.unknown
_951320099.unknown
_951320174.unknown
_1172241793.unknown
_1172250127.unknown
_951320348.unknown
_951399058.unknown
_951320127.unknown
_951318866.unknown
_951319722.unknown
_951319990.unknown
_951319623.unknown
_951318509.unknown
_951153353.unknown
_951221261.unknown
_951225969.unknown
_951226600.unknown
_951311171.unknown
_951231672.unknown
_951226551.unknown
_951221298.unknown
_951221213.unknown
_951221247.unknown
_951153371.unknown
_951153310.unknown
_951112515.unknown
_951144983.unknown
_951112456.unknown
_951107733.unknown
_951111029.unknown
_951112089.unknown
_951112007.unknown
_951112077.unknown
_951107933.unknown
_951110960.unknown
_951110998.unknown
_951107878.unknown
_951107821.unknown
_802411687.unknown
_951106069.unknown
_951106826.unknown
_951107030.unknown
_951107074.unknown
_951107230.unknown
_951106985.unknown
_951106679.unknown
_951106744.unknown
_951106184.unknown
_802411782.unknown
_951105813.unknown
_802824627.unknown
_802411759.unknown
_802179423.unknown
_802325095.unknown
_802331923.unknown
_802332116.unknown
_802332868.unknown
_802332006.unknown
_802325118.unknown
_802179578.unknown
_802325053.unknown
_802179544.unknown
_802179463.unknown
_801214756.unknown
_802164121.unknown
_802164899.unknown
_802164922.unknown
_802166875.unknown
_802164841.unknown
_801215108.unknown
_801599933.unknown
_801599962.unknown
_801215208.unknown
_801215043.unknown
_801213926.unknown
_801214096.unknown
_801135021.unknown
_801213704.unknown
_801134984.unknown