第五章 主成分分析(1)(主成分模型).doc

第五章 主成分分析与经验正交分解5.1 主分量分析的数学模型当存在若干个随机变量时，寻求它们的少量线性组合（即主成分），用以解释这些随机变量，是很必要的。首先我们看一个例子。几个数据集1、（1） 身材情况能否用单个指标刻画（2） 男女身材之间有什么异同 chest waist hips gender chest waist hips gender34 30 32 male 36 24 35 female37 32 37 male 36 25 37 female38 30 36 male 34 24 37 female36 33 39 male 33 22 34 female38 29 33 male 36 26 38 female43 32 38 male 37 26 37 female40 33 42 male 34 25 38 female38 30 40 male 36 26 37 female40 30 37 male 38 28 40 female41 32 39 male 35 23 35 female2、subject maths english history geography chemistry physics1 60 70 75 58 53 422 80 65 66 75 70 763 53 60 50 48 45 434 85 79 71 77 68 795 45 80 80 84 44 463、air pollution in cities in the USA. The following variables were obtained for 1 US cities:SO2: SO2 content of air in micrograms per cubic metre;temp: average annual temperature in degrees Fahrenheit;manu: number of manufacturing enterprises employing 20 or more workers;popul: population size (1970 census) in thousands;wind: average annual wind speed in miles per hour;precip: average annual precipitation in inches;predays: average number of days with precipitation per year.例51 为了调查学生的身材状况，可以测量他们的身高()、体重()、胸围()和坐高()。可是用这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的内积。对第三，第四主成分同样可证。例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.2 相关系数阵和协方差阵的主分量分析在实际问题中，X的每一分量可取不同单位，单位取小时（例如长度单位取毫米，甚至微米）该分量的方差会变大，从而在主成分中变得突出；而单位选取不应影响主成分。为了避免量纲对主成分的影响。常常将随机变量都标化，即令，它就是无量纲量，令再求X*的主成分，即标准化后的主成分。将代入，可求随机向量X的主成分。容易证明定理5.2 设随机向量X的相关阵为，特征值为,对应的彼此正交单位特征向量为，则标准化后X的第j个主成分是。因此，标准化后的主成分称为由相关阵决定的主成分。直接由随机向量的协方差阵算出的主成分称为由协差阵决定的主成分。同样一组随机变量，用它们的协差阵和相关阵求出的主成分是不一样的。这是因为优化的准则（目标函数）不同：前者要求=最大，而后者要求=最大，其中。例5.3 （协差阵和相关阵决定的主成分不同）设随机变量；其协方差阵是，特征值和特征向量是， 。因而由协方差阵决定的主成分是： ,。但随机变量标准化后得到；其中 。X*的协差阵即X的相关阵是，其特征值和特征向量是，从而由相关阵决定的主成分是：。由于主成分由方差决定，可以略去常数，因而由相关阵得到的主成分可写为：，可见由协方差阵与相关阵决定的主成分不同。分别从相关系数阵和协方差阵计算主成分的例子。该例取自Jolliffe (2002). 该数据包含72位病人的八项血液化验指标。data blood_corr(type=cov); input rblood plate wblood neut lymph bilir sodium potass;cards; 1.000 0.290 0.202 -0.055 -0.105 -0.252 -0.229 0.0580.290 1.000 0.415 0.285 -0.376 -0.349 -0.164 -0.1290.202 0.415 1.000 0.419 -0.521 -0.441 -0.145 -0.076-0.055 0.285 0.419 1.000 -0.877 -0.076 0.023 -0.131-0.105 -0.376 -0.521 -0.877 1.000 0.206 0.034 0.151-0.252 -0.349 -0.441 -0.076 0.206 1.000 0.192 0.077-0.229 -0.164 -0.145 0.023 0.034 0.192 1.000 0.4230.058 -0.129 -0.076 -0.131 0.151 0.077 0.423 1.000;proc princomp; /*用相关系数计算主成分*/var rblood plate wblood neut lymph bilir sodium potass; run;data blood_cov(type=cov); input rblood plate wblood neut lymph bilir sodium potass;cards; 0.137641000 4.4384103 0.14501277 -0.001571185 -0.002765805 -0.37742720 -0.232107988 0.0063908464.438410270 1701.8100090 33.12719033 0.905297085 -1.101290088 -58.12188799 -18.483324144 -1.5805261890.145012770 33.1271903 3.74422500 0.062428905 -0.071577585 -3.44491339 -0.766530900 -0.043676820-0.001571185 0.9052971 0.06242890 0.005929000 -0.004794559 -0.02362452 0.004838372 -0.002995839-0.002765805 -1.1012901 -0.07157758 -0.004794559 0.005041000 0.05904516 0.006595048 0.003184137-0.377427204 -58.1218880 -3.44491339 -0.023624524 0.059045162 16.29736900 2.117584128 0.092322153-0.232107988 -18.4833241 -0.76653090 0.004838372 0.006595048 2.11758413 7.463824000 0.3432238920.006390846 -1.5805262 -0.04367682 -0.002995839 0.003184137 0.09232215 0.343223892 0.088209000;proc princomp cov; /*用协方差计算主成分*/var rblood plate wblood neut lymph bilir sodium potass; run;5.3 主成分个数的确定下面讨论对变量很多时，选择多少个主成分。它称为第i个主成分的方差贡献，表示第i个主成分变化大小，从而反映第i个主成分提供的信息的大小。定义5.3 称为主成分的方差贡献率；称为前k个主成分的累计方差贡献率；与X第k个分量的相关系数称为因子负荷量。当某个主成分的方差贡献率很小时，认为它提供的信息很少，可以略去此主成分。通常取q,使前q个主成分的累计方差贡献率达到70%-80%,然后只考虑前q个主分量，用它们解释随机向量X的特性，其余主成分认为是观测误差等随机因素造成的。Jolliffe (2002), Rencher (2002).给出了许多实例，一些常用的选择主成分个数的法则如下：1. Retain just enough components to explain some speci_ed large percentage of the total variation of the original variables. Values between 70% and 90% are usually suggested, although smaller values might be appropriate as q or n, the sample size, increases.2 . Exclude those principal components whose eigenvalues are less than the average, the average eigenvalue is also the average variance of the original variables. This method then retains·those components that account for more variance than the average for the observed variables.3. When the components are extracted from the correlation matrix, trace(R) = q, and the average variance is therefore one, so applying the rule in the previous bullet point, components with eigenvalues less than one are excluded. This rule was originally suggested by Kaiser (1958), but Jolliffe (1972), on the basis of a number of simulation studies, proposed that a more appropriate procedure would be to exclude components extracted from a correlation matrix whose associated eigenvalues are less than 0.7.4. Cattell (1966) suggests examination of the plot of the against i, the socalled scree diagram. The number of components selected is the value of I corresponding to anelbow"in the curve, i.e., a change of slope fromsteep" to shallow". In fact, Cattell was more speci_c than this, recommending to look for a point on the plot beyond which the scree diagram de_nes a more or less straight line, not necessarily horizontal. The _rst point on the straight line is then taken to be the last component to be retained. And it should also be remembered that Cattell suggested the scree diagram in the context of factor analysis rather than applied to principal components analysis.5. A modification of the scree digram described by Farmer (1971) is the log-eigenvalue diagram consisting of a plot of against i.Returning to the results of the principal components analysis of the blood chemistry data given in Section 3.3, we find that the first four components account for nearly 80% of the total variance, but it takes a further two components to push this figure up to 90%. A cutoffb of one for the eigenvalues leads to retaining three components, and with a cuto_ of 0.7 four components are kept. Figure 3.1 shows the scree diagram and log-eigenvalue diagram for the data5.4 样本主成分实际问题中随机向量的协差阵、相关阵都是未知的，只能得到样品。这时总用样本协差阵与样本相关阵代替协差阵、相关阵求主成分。定义5.4 样本协差阵与样本相关阵的特征向量，计算主成分。所得的主成分称为样本主成分。这样求主成分是有道理的：若总体，的特征值和正交单位特征向量是和；是的极大似然估计，即。的特征值为，相应正交单位特征向量为，则可证定理5.3 若X服从正态分布，则是的极大似然估计；是的极大似然估计。因此，若X服从正态分布，应当用第j个样本主成分作为总体主成分的估计值。从样本协差阵或样本相关阵出发，做主成分分析，所得样本主成分通常简称为主成分。通常取为样本协差阵（的无偏估计），由或R算出的样本相关阵是相同的，所产生（相关差阵决定）的主成分当然相同。而R与有相同的特征向量，R的特征值是特征值的n/(n-1)倍。因而由R与所产生的（协方差阵决定的）主成分相同。若X不一定服从正态分布，这时仍可由样本协差阵R或相关阵出发，计算主成分。同上节指出的一样：样本相关阵和样本协差阵决定的主成分是不同的。5.5 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 北京19511976年冬季的气温资料如表5-1，第一列为年度，第二列为该年12月的月平均温度。第三、四列为次年1、2月的月平均温度。试做主成分分析。表 5-1 北京19511976年冬季月平均气温yearx1x2x319511.0-2.7-4.31952-5.3-5.9-3.51953-2.0-3.4-0.81954-5.7-4.7-1.11955-0.9-3.8-3.11956-5.7-5.3-5.91957-2.1-5.0-1.619580.6-4.3-0.21959-1.7-5.72.01960-3.6-3.61.31961-3.0-3.1-0.819620.1-3.9-1.11963-2.6-3.0-5.21964-1.4-4.9-1.71965-3.9-5.7-2.51966-4.7-4.8-3.31967-6.0-5.6-4.91968-1.7-6.4-5.11969-3.4-5.6-2.01970-3.1-4.2-2.91971-3.8-4.9-3.91972-2.0-4.1-2.41973-1.7-4.2-2.01974-3.6-3.3-2.01975-2.7-3.70.11976-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.31952 -5.3 -5.9 -3.51953 -2.0 -3.4 -0.81954 -5.7 -4.7 -1.11955 -0.9 -3.8 -3.11956 -5.7 -5.3 -5.91957 -2.1 -5.0 -1.61958 0.6 -4.3 -0.21959 -1.7 -5.7 2.01960 -3.6 -3.6 1.31961 -3.0 -3.1 -0.81962 0.1 -3.9 -1.11963 -2.6 -3.0 -5.21964 -1.4 -4.9 -1.71965 -3.9 -5.7 -2.51966 -4.7 -4.8 -3.31967 -6.0 -5.6 -4.91968 -1.7 -6.4 -5.11969 -3.4 -5.6 -2.01970 -3.1 -4.2 -2.91971 -3.8 -4.9 -3.91972 -2.0 -4.1 -2.41973 -1.7 -4.2 -2.01974 -3.6 -3.3 -2.01975 -2.7 -3.7 0.11976 -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月温度距平所得值的反差。例55 美国各州犯罪率情况如表5-2。试以murder（谋杀）,rape（强奸）,robbery（抢劫）,assult（斗殴）,burglary（夜盗）,larceny（偷窃）,auto（汽车犯罪）为7元随机向量，做主成分分析。表 5-2 美国各州犯罪率（十万人中犯罪人数）murderraperobberyassultburglarylarcenyautoAlbama14.225.296.8278.31135.51881.9280.7Alaska10.851.696.8284.01331.73369.8753.3Arirona9.534.2138.2312.32346.14467.4439.5Arkansas8.834.2138.2312.32346.14467.4439.5Califonia11.549.4287.0358.02139.43499.8663.5Colorado6.342.0170.7292.91935.23903.2477.1Conecticat4.216.8129.5131.81346.02620.7593.2Delaware6.024.9157.0194.21682.63678.4467.0Florida10.239.6187.9449.11859.93840.5351.4Geogia11.731.1140.5256.51351.12170.2297.9Hawaii7.225.5128.064.11911.53920.4489.4Idaho5.519.439.6172.51050.82599.6237.6Illinois9.921.8211.3209.01085.02828.5528.6Indiana7.426.5123.2153.51086.22498.7377.4Iowa2.310.641.289.8812.52685.1219.9Kansas6.622.0100.7180.51270.42739.3244.3Kentaky10.119.181.1123.3872.21662.1245.4Loisana15.530.9142.9335.51165.52469.9337.7Maine2.413.538.7170.01253.12350.7246.9Maryland8.034.8292.1358.91400.03177.7428.5Masschusetts3.120.8169.1231.61532.22311.31140.1Michigan9.338.9261.9274.61522.73159.0545.5Minnesota2.719.585.985.81134.72559.3343.1Mississippi14.319.665.7189.1915.61239.9144.4Missouri9.628.3189.0233.51318.32424.2378.4Montana5.416.739.2156.8804.92773.2309.3Nebraska3.918.164.7112.7760.02316.1249.1Nevada15.849.1323.1355.02453.14212.6559.2Mew Hampashare3.210.723.276.01041.72343.9293.4New Jersey 5.621.0180.4185.11435.82774.5511.5New Maxico8.839.1109.6343.41418.73008.6259.5New York10.729.4472.6319.11728.02782.0745.8North Carolina10.617.061.3318.31154.12037.8192.1North Dakoda100.99.013.343.8446.11843.0