MATLAB - 主成分数据处理

主成分分析

主成分分析（principal component Analysis），是由皮尔逊（pearson）于1901年首先引入，后来由霍特林（hotelling）于1933年进行了发展。

在实际问题中，为了尽可能完整的获取有关的信息，往往需要考虑众多的变量，这虽然可以避免重要信息的疏漏，但也增加了分析的复杂性，一般来说，当研究的问题涉及很多变量，并且变量间相关性明显，即包含的信息有所重叠时，我们就会很自然地想到，能否在各个变量之间相关关系研究的基础上，用较少的新变量代替原来较多的变量，而且使这些较少的新变量尽可能多地保留原来较多的变量所反映的信息？事实上，这种想法是可以实现的，本节拟介绍的主成分分析方法就是综合处理这种问题的一种强有力的方法。这样容易抓住事物的主要矛盾，使得问题得到简化。

主成分分析是一种通过降维技术把多个变量化为少数几个主成分（即综合变量）的多元统计方法，这些主成分能够反映原始变量的大部分信息，通常表示为原始变量的线性组合，为使得这些主成分所包含的信息互不重叠，要求各主成分之间互不相关。

本章主要内容包括：主成分分析的理论简介，主成分分析的MATLAB实现，主成分分析的主要具体案例。

11.1主成分分析简介

11.1.1主成分分析的几何意义

假设从二元总体x?(x1,x2)' （EX=0）中抽取容量为n的样本，绘出样本观测值的散点图，散点大致分布在一个椭圆内x1与x2呈现出明显的线性相关。这n个样品在x1轴方向和x2方向具有相似的离散度，离散度可以用x1和x2包含了近视相等的信息量，丢掉其中任意一个变量，都会损失比较多的信息。逆时针旋转一个角度?，使得x1轴旋转到椭圆的长轴方向

y1，x2轴旋转到椭圆的短轴y2，则有

?y1?x1cos??x2sin?? （11.1） y??xsin??xcos??212此时可以看到，n个点在新坐标系下的坐标

y1和y2几乎不相关，并且

y1的方差要比y2的方差大得多，也就是说y1包含了原始数据中大部

分的信息，此时丢掉变量第一主成分

y2，信息的损失是比较小的。这里称y1为

y2为第二主成分。

主成分分析的目的就是对原变量加以改造，在不致损失原变量太多信息的情况下尽可能地降低原变量的维数，即用较少的新变量代替原来的各变量。

主成分分析的过程其实就是坐标系旋转的过程，新坐标系的各个坐标系的轴的方向是原始数据变差最大的方向，各主成分表达式就是新旧坐标转换关系式。

11.1.2 总体的主成分

1、从总体协方差矩阵出发求解主成分

'x?(x,x,?,x)设p为一个p维总体，假定x期望和协方差矩阵均存

TE(x)??var(x)???E[(X?EX)(X?EX)]，在并已知，记，

考虑如下线性变换

?y1?a11x1?a12x2?...?a1pxp?a'1x?y?ax?ax?...?ax?a'x2112222pm2?2?? ??yp?ap1x1?ap2x2?...?appxp?a'px?其中，a1,a2,?,ap均为单位向量。 下面求a1，使得

y1的方差达到最大。

p个特征值，t1,t2,?tp为相应的正

设?1??2????p?0为?的交单位特征向量，即?ti??iti，ti'ti?1，ti'tj?0，i?j，

i,j?1,2,?,p

由矩阵知识可知

??T?T'???itit'ii?1p

其中T?(t1,t2,?tp)为正交矩阵，角矩阵。 考虑

?是对角元素为?1,?2,?,?p的对

y1的方差

var(y1)?var(a'1x)?c'1var(x)a1???ia'1tit'ia1

i?1p?p????i(a'1ti)??1?(a'1ti)??1a'1??tit'i?a1 i?1i?1?i?1?22pp??1a'1TT'a1??1a'1a1??1

（11.3）

由式（11.3）可知，当a1?t1时，y1?t'1x的方差达到最大，最大

值为?1。称y1?t'1x为第一主成分。如果第一主成分从数据中提取的信息还不够多，还应考虑第二主成分。下面求a2，在cov(y1,y2)?0条件下使得

y2的方差达到最大。由

cov(y1,y2)?cov(t'1x,a'2x)?t'1?a2?a'2?t1??1a'2t1?0

可得a'2t1?0，于是

p

var(y1)?var(a'2x)?a'2var(x)a2???ia'2tit'ia2

i?1ppp??22???i(a'2ti)??2?(a'2ti)??2a'2??tit'i?a2 i?1i?1?i?1???2a'2TT'a2??2a'2a2??2

（11.4）

由式（11.4）可知，当a2?t2时，值为?2。称

y2?t'2x的方差达到最大，最大

y2?t'2x为第二主成分。类似的，在约束

cov(yk,yi)?0(k?1,2,?,i?1)下可得，当ai?ti时yi?t'ix的方差达到最大，最大值为。?i称yi?t'ix(i?1,2,?,p) 为第i主成分。 2主成分的性质

（1）主成分向量的协方差矩阵为对角阵 记

?y1??t'1x?????y2??t'2x??y???(t1,t2,?,tp)'x?T'x?????? ?????y??t'x??p??p?（11.5） 则 E(y)?E(T'x)?T'?，

var(y)?var(T'x)?T'var(x)T?T'?T??

即主成分向量的协方差矩阵为对角矩阵。 （2）主成分的总方差等于原始变量的总方差：

设协方差矩阵??(?ij)，则var(xi)??ij(i?1,2,?,p)，于是

?var(y)???ii?1i?1ppi?tr(?)???ij??var(xi)

i?1i?1pp由此可见，原始数据的总方差等于和，也就是说

p个互不相关的主成分的方差之

p个互不相关的主成分包含了原始数据中的全部信息，

p但是主成分所包含的信息更为集中。

总方差中第i个主成分称为主成分

yi的方差所占的比例?i??j(i?1,2,?,p)j?1yi的贡献率。主成分的贡献率反映了主成分综合原始变量

p信息的能力，也可理解为解释原始变量的能力。由贡献率定义知，

个主成分的贡献率依次递减，即综合原始变量信息的能力依次递减。第一个主成分的贡献率最大，即第一个主成分综合原始变量信息的能力强。

前m(m?p)个主成分的贡献率之和??ii?1m??j?1pj称为前m个主成分的累积贡献率，它反映了前m个主成分综合原始变量信息（或解释原始变量）的能力。由于主成分分析的主要目的是降维，所以需要在信息损失不太多的情况下，用少数几个主成分来代替原始变量x1,x2,?,xp，以进行后续的分析，究竟用几个主成分来代替原始变量才合适呢？通常的做法是取较小的m，使得恰前m个主成分的累积贡献率不低于某一水平（如85%以上），这样就达到了降维的目的。

（3）原始变量xi与主成分

yi之间的相关系数?(xi,yi)

由式（11.5）可知x?Ty于是 （11.6） 从而

cov(xi,yi)?cov(tijyi,yj)?tijcov(yj,yj)?tij?j

xi?ti1y1?ti2y2???tipyp

?(xi,yi)?cov(xi,yi)var(xi)var(yi)??j?iitij,i,j?1,2,?,p

（4）前m个主成分对变量xi的贡献率 称

??(xi,yi)?2j?1m1?ii2?t?jij j?1m为前m个主成分对变量xi的贡献率。这个贡献率反映了前m个主成分从变量xi中提取的信息的多少。由式（11.6）可知

2?ii??1ti21??2ti22????ptip，固所有p个主成分对变量xi的贡献

率为

2?(xi,yi)? ?j?1p1?ii2?t?jij?1 j?1p（5）原始变量对主成分

yi的贡献

主成分

yi的表达式为

yj?t'j?t1jx1?t2jx2???tpjxp,j?1,2,?,p

称tij为第j个主成分yj在第i个原始变量xi上的载荷，它反映了xi对yj的重要程度。在实际问题中，通常根据载荷tij解释主成分的实际意义。

3，从总体相关系数矩阵出发求解主成分

当总体各变量取值的单位或数量级不同时，从总体协方差矩阵出发求解主成分就显得不合适了，此时应将每个变量标准化。记标准化变量为

x??xi?E(xi)var(xi),i?1,2,?,p

**x??(x1,x2,?x*p)'则可以从标准化总体的协方差矩阵求解主

?xx成分，即从总体的相关系数矩阵出发求解主成分，因为总体协方

差矩阵就是x的相关系数矩阵。

设总体x的相关系数矩阵为R，从R出发求解主成分的步骤与从

***?出发求解主成分的步骤一样，设?1??2???p?0为R的p个特征

***t,t,?,t值，12p为相应的正交单位特征向量。则p个主成分为

y（11.7）

*i?t'x,i?1,2,?,p

*i??y1???t?'1x?????????y2??t'2x?**??*?y?????(t1,t2,?,tp)'x????记

????????y??t'px???p??（11.8） 则有以下结论

?**** E(y?)?0，var(y)???diag(t1,t2,?,tp)

p??i?1*i?tr(R)?P

?(xi*,y*j)?cov(xi*,y*j)var(xi*)var(y*j)*??*tjij,i,j?1,2,?,p

1m*此时前m个主成分的累积贡献率为??i。

pi?111.1.3 样本的主成分

在实际生活问题中，总体x的协方差矩阵?或相关系数矩阵R往往是未知的，需要由样本进行估计。设x1,x2,?,xn为取自总体x的样本，其中xi?(xi1,xi2,?,xip)'(i?1,2,?,n)。记样本观测值矩阵为

?x11x12?x1p???xx?x21222p?X??????? ???x?x?xn1n2np??X的每一行对应一个样品，每一列对应一个变量。记样本协方差矩阵

和样本相关系数矩阵分别为

1n(xi?x)(xi?x)'?(sij) S??n?1i?1 R?(rij)，rij??sijsiisjj

1n?其中x??xi为样本均值。将S作为?的估计，R作为R的估

ni?1计，从S或R出发可求得样本的主成分

1.从样本协方差矩阵S出发求解主成分 设?1??2???p?0为S的交单位特征向量，则样本的

???p个特征值，t1,t2,?,tp为相应的正

?????p个主成分为

yi?t'ix,i?1,2,?,p （11.9）

将样品xi的观测值带入第

??j个主成分，称得到的值

为样品xi的观测值带入第j主yij?t'jx,(i?1,2,?,n;j?1,2,?,p)成分得分。

从样本相关系数矩阵R出发求解主成分 设???1???????0为R的p个特征值，t?,t?,?,t????2??p???12p为相应

的正交单位特征向量，则样本的

p个主成分为

y?t'x?,i?1,2,?,p??i??i

（11.10）

将样品xi标准化后的观测值xi带入第的第

j个主成分，即可得样品xij主成分得分

y?t'x?,i?1,2,?,n;j?1,2,?,p

??ij??j3，由主成分得分重建（恢复）原始数据

假定从样本协方差矩阵S出发求解主成分，记Y为样本的主成分得分值矩阵，则

?y12?y1p??x11????xy22?y2p??21?????????x??n1?yn2?ynp??（11.11）

?????y11???yY??21?????y?n1??x12?x1p????x22?x2p???(t1,t2,?,tp)?XT????xn2?xnp??注意到T为正交矩阵，则有T?T'，于是由式（11.11）可得

??1?X?YT'，也就是说根据主成分得分和主成分表达式，可以重建（恢

复）原始数据，这在数据压缩与解压缩中有着重要的作用。当然在实际应用中，可能不会得到全部的

??p个主成分，假定只用前

m(m?p)个主成分记样本的前m个主成分的得分矩阵为

?????yyLy121p??11??????yy22Ly2p?Ym??21?MMM? ??????y?yLyn2np??n1当前m个主成分的累积贡献率达到一个比较高的水平时，由

Xm?YmT'得到的矩阵Xm可以作为原始样本观测值矩阵X的一个很好的

??近视，此时X?Xm为样本的残差，MATLAB统计工具箱中提供了重建数据和求残差的函数pcares。若Ym和T'的数据量小于原始样本观测值矩阵X的数据量，就能起到数据压缩的目的。

以上讨论的是从样本协方差矩阵S出发求解主成分，然后由样本的主成分得分重建原始数据。若从样本的相关系数矩阵R出发求解主成分，同样可以由样本的主成分得分重建原始数据，只是此时需要进行逆标准化变换，这里不再作详细讨论。

11.1.4关于主成分表达式的两点说明

这里需要说明的是，即使限定了方差矩阵或相关系数矩阵的

???p个

特征值对应的特征向量为正交单位向量，它们也是不唯一的，从而主成分的表达式也是不唯一的，假如若y?t'x是总体或样本的一个主成分，则y??t'x也是总体或样本的一个主成分。主成分表达式的不唯一对后续分析没有太大影响。

若第

p个主成分的贡献率非常非常小，可认为第p个主成分yp的

方差var(yp)?0，即yp?c（c为一个常数），这揭示了变量之间的一

't个共线性关系：px?c。

11.2 主成分分析的MATLAB函数

与主成分相关的MATLAB函数主要有pcacov,princomp和pcares,下面分别介绍。 11.2．1 pcacov函数

pcacov函数用来根据协方差矩阵或相关系数矩阵进行主成分分析，其调用格式如下：

COEFF=pacov(v)

[COEFF,latent]= pcacov(v)

[COEFF,latent,explained]=pcacov(v)

以上调用的输入参数V是总体或样本的协方差矩阵或相关系数矩阵，对于

p维总体，V是p?p的矩阵。输出参数COEFF是p个主成

p个主成分的方差构成的列向量，即V的

分的系数矩阵，它是p?p的矩阵，它的第i列是第i个主成分的系数向量。输出参数latent是

p个特征值（从大到小）构成的向量。输出参数explained是p个

主成分的贡献率向量，已经转化为百分比。

11.2.2 princomp函数

princomp函数用来根据样本观测值矩阵进行主成分分析，其调用格式如下：

1）[COEFF,SCORE]=princomp(x)

根据样本观测值矩阵X进行主成分分析。输入参数X是n行p列的矩阵，每一行对应一个观测（样品），每一列对应一个变量。输出参数COEFF是

p个主成分的系数矩阵，它是p?p的矩阵，它的第i

p个

列是第i个主成分的系数向量。输出参数SCORE是n个样品的

主成分的得分矩阵，它是n行p列的矩阵，每一行对应一个观测值，每一列对应一个主成分，第i行第j列元素是i个样品的第j个主成分得分。

2）[COEFF,SCORE,latent]=princomp(x)

返回样本协方差矩阵的特征向量latent，它是由p个特征值构成的列向量，其中特征值按降序排列。

3）[COEFF,SCORE,latent,tsqure]=princomp(x) 返回一个包含p个元素的列向量tsqure，它的第i个元素是第i个观测对应的霍特林（Hotelling）T2统计量，描述了第i个观测与数据集（样本观测矩阵）的中心之间的距离，可用来寻找远离中心的极端数据。

设?1??2???p?0为样本协方差矩阵的p个特征值，并设第i个样品的第j个主成分得分yij(i?1,2,?,n;j?1,2,?,p)，则第i个样品对应的（Hotelling）T2统计量为

???Ti2??j?1p2yij?j,i?1,2,?,n

注意：princomp函数对样本数据进行了中心化处理，即把X中的每一个元素减去其所在列的均值，相应地，princomp函数返回的主成分就是中心化的主成分得分。

当n?p，即观测的个数小于或等于维数时，SCORE矩阵的第n列到第p列元素均为0，latent第n到第p个元素均为0。

4）[...]=princomp(x,'econ')

通过设置'econ'参数，使得当n?p时，只返回latent中的前n-1个元素（去掉不必要的0元素）及COEFF和SCORE矩阵中相应的列。

11.2.3 pcares函数

在11.1.3节曾讨论过由样本的主成分得分重建（恢复）原始数据的问题，若只用前m(m?p)个主成分的得分来重建原始数据，

则可能会有一定的误差，前面称之为残差。MATLAB统计工具箱中提供了pcares函数，用来重建数据，并求样本观测值矩阵中每个观测的每一个分量所对应的残差，其调用格式如下：

residuals=pcares(x,ndim)

[residuals,reconstructed]=pcares(x,ndim) 上述调用中X是n行p列的样本观测值矩阵，它的每一行对应一个观测（样品），每一列对应一个变量，ndim参数用来指定所用的主成分的个数，它是一个小于或等于p的正的标量，最好取为正整数。输出参数residuals是一个与X同样大小的矩阵，其元素为X中相应元素所对应的残差。输出参数reconstructed为用前ndim个主成分的得分重建的观测数据，它是X的一个近似。

注意：pcares调用了 princomp函数，它只能接受原始样本观测数据作为他的输入，并且它不会自动对数据作标准化变换，若需要对数据作标准化变换，可以先用zscore函数将数据标准化，然后调

3 4 5 6 7 8 9 10 11 12 13 14 13.2 22.3 34.3 35.6 22.0 48.4 40.6 24.8 12.5 1.8 32.3 38.5 3.3 6.7 11.8 12.5 7.8 13.4 19.1 8.0 9.7 0.6 13.9 9.1 3.9 5.6 7.1 16.4 9.9 10.9 19.8 9.8 4.2 0.7 9.4 11.3 4.3 3.7 7.1 16.7 10.2 9.9 19.0 8.9 4.2 0.7 8.3 9.5 4.4 6.0 8.0 22.8 12.6 10.9 29.7 11.9 4.6 0.8 9.8 12.2 5.5 7.4 8.9 29.3 17.6 13.9 39.6 16.2 6.5 1.1 13.3 16.4 0.578 0.176 1.726 3.017 0.847 1.772 2.449 0.789 0.874 0.056 2.126 1.327 3.6 7.3 27.5 26.6 10.6 17.8 35.8 13.7 3.9 1.0 17.1 11.6 解：样本均值向量为：

x?(27.97910.9509.1008.54311.06414.6141.55214.686)T,

样本协方差矩阵为：

?168.33360.35745.75741.21557.90671.6728.602101.620???37.20716.82515.50523.53529.0294.78544.023???24.84324.33536.47849.2783.62939.410???24.42336.28349.1463.67538.718? S???56.04675.4045.00259.723???103.0186.82174.523???1.1376.722???102.707?????168.3360.35745.75841.21657.90671.6728.602101.62??60.35737.20716.82515.50523.53529.0294.784644.023????45.75816.82524.84324.33536.47849.2783.62939.41???41.21615.50524.33524.42336.28349.1463.674738.718? S???57.90623.53536.47836.28356.04675.4045.002259.723???71.67229.02949.27849.14675.404103.026.821574.523???8.6024.78463.6293.67475.00226.82151.1376.7217?????101.6244.02339.4138.71859.72374.5236.7217102.71??由于S中主对角线元素差异较大，因此我们样本相关矩阵R出发

进行主成分分析。样本相关矩阵R为：

?1 0.76266 0.70758 0.64281 0.59617 0.54426 0.62178 0.77285???1 0.553410.51434 0.51538 0.468880.73562 0.71214???10.98793 0.9776 0.974090.68282 0.78019??? 1 0.98071 0.97980.69735 0.77306? R??? 1 0.992350.62663 0.78718??? 10.6303 0.72449??? 1 0.62202?????? 1?矩阵R的特征值及相应的特征向量分别为： 特征值 6.1366 1.0421 0.43595 0.22037 0.15191 0.0088274 0.0029624 0.0012238 特征向量 0.32113 0.29516 0.38912 0.38472 0.37955 0.37087 0.31996 0.35546 -0.4151 -0.59766 0.22974 0.27869 0.31632 0.37151 -0.27814 -0.15684 -0.45123 0.10303 -0.039895 0.053874 -0.037292 0.075186 0.77059 -0.42478 -0.66817 0.36336 -0.22596 -0.11081 0.14874 0.069353 -0.13495 0.55949 -0.038217 0.62435 0.12273 -0.036909 0.15928 0.21062 -0.43006 -0.58105 -0.10167 0.13584 -0.15811 0.86226 -0.25204 -0.34506 -0.13934 -0.026557 0.1596 -0.061134 -0.53966 0.046606 0.7609 -0.27809 0.06203 -0.13126 0.19295 -0.031987 -0.64176 0.11002 -0.25397 0.68791 -0.006045 -0.0054031 特征值 6.1366 1.0421 0.43595 0.22037 0.15191 0.0088274 0.0029624 0.0012238 R的特征值及贡献率见下表 贡献率（%） 累计贡献率（%） 0.76708 0.13027 0.054494 0.027547 0.018988 0.0011034 0.0003703 0.00015297 0.76708 0.89734 0.95184 0.97938 0.99837 0.99948 0.99985 1 前3个标准化样本主成分类及贡献率已达到95.184%，故只需取前三个主成分即可。

前3个标准化样本主成分中各标准化变量 xi*?xi?xisii(i?1,2,...,8)前的系数即为对应特征向量，由此得到3个标准化样本主成分为

*******?y1?0.32113x1+0.29516x*2+0.38912x3+0.38472x4+0.37955x5+0.37087x6+0.31996x7+0.35546x8?********?y2?-0.4151x1-0.59766x2+0.22974x3+0.27869x4+0.31632x5+0.37151x6-0.27814x7-0.15684x8?********y?-0.45123x+0.10303x-0.039895x+0.053874x-0.037292x+0.075186x+0.77059x-0.42478x312345678?

注意到，y1近似是8个标准化变量xi*?xi?xisii(i?1,2,...,8)的等权重之

和，是反映各企业总效应大小的综合指标，y1的值越大，则企业的效益越好。由于y1的贡献率高达76.708%，故若用y1的得分值对各企业进行排序，能从整体上反映企业之间的效应差别。将S中sii的值及x中各xi的值以及各企业关于xi的观测值代入y1的表达式中，可求得各企业y1的得分及其按其得分由大到小的排序结果。

企业序号 得分 12 -0.97354 4 -0.64856 3 -0.62743 11 -0.48558 10 -0.21949 7 -0.189 14 -0.004803 5 0.016879 8 0.17711 13 0.18925 1 0.29351 2 0.65315 6 0.85566 9 0.96285

所以，第9家企业的效益最好，

第12家企业的效益最差。

Matlab程序：[coeff,score,latent]=princomp(X) 注：该函数使用协方差阵作主成分分析。

案例33：从样本观测值矩阵出发求解主成分

表11-2列出了2007年我国31个省，市，自治区和直辖市的农村居民家庭平均每人全年消费支出的8个主要变量数据。数据来源：中华人民共和国国家统计局网站，2008年《中国统计年鉴》。数据格式如表11-2所列，是根据这8个主要变量的观测数据，进行主成分分析

家庭设地 区 食 品 衣 着 居 住 备 及 服 务 北 京 天 津 河 北 山 西 内蒙古 辽 宁 吉 林 黑龙江 上 海 江 苏 浙 江 安 徽 福 建 江 西 2132.51 513.44 1023.21 1367.75 286.33 674.81 1025.72 185.68 627.98 1033.68 260.88 392.78 1280.05 228.40 473.98 1334.18 281.19 513.11 1240.93 227.96 399.11 1077.34 254.01 691.02 3259.48 475.51 2097.21 1968.88 251.29 752.73 2430.60 405.32 1498.50 1192.57 166.31 479.46 1870.32 235.61 660.55 1492.02 147.71 474.49 文教娱交通和 乐 通 讯 用品及服务 医疗保健 其他商品 及 服 务 340.15 778.52 870.12 629.56 111.75 126.74 400.11 312.07 306.19 64.30 140.45 318.19 243.30 188.06 57.40 120.86 268.75 370.97 170.85 63.81 117.64 375.58 423.75 281.46 75.29 142.07 361.77 362.78 265.01 108.05 120.95 337.46 339.77 311.37 87.89 104.99 335.28 312.32 272.49 69.98 451.40 883.71 857.47 571.06 249.04 228.51 543.97 642.52 263.85 134.41 338.80 782.98 750.69 452.44 142.26 144.23 258.29 283.17 177.04 52.98 184.21 465.40 356.26 174.12 107.00 121.54 277.15 252.78 167.71 61.08 山 东 河 南 湖 北 湖 南 广 东 广 西 海 南 重 庆 四 川 贵 州 云 南 西 藏 陕 西 甘 肃 青 海 宁 夏 新 疆 1369.20 224.18 682.13 1017.43 189.71 615.62 1479.04 168.64 434.91 1675.16 161.79 508.33 2087.58 162.33 763.01 1378.78 86.90 1430.31 86.26 554.14 305.90 195.99 422.36 424.89 230.84 71.98 136.37 269.46 212.36 173.19 62.26 166.25 281.12 284.13 178.77 97.13 152.60 278.78 293.89 219.95 86.88 163.85 443.24 254.94 199.31 128.06 112.24 245.97 172.45 149.01 47.98 93.26 248.08 223.98 95.55 73.23 1376.00 136.34 263.73 1435.52 156.65 366.45 998.39 99.44 329.64 138.34 208.69 195.97 168.57 39.06 142.64 241.49 177.19 174.75 52.56 70.93 154.52 147.31 79.31 34.16 1226.69 112.52 586.07 1079.83 245.00 418.83 941.81 944.14 161.08 512.40 112.20 295.23 107.15 216.67 181.73 167.92 38.43 133.26 156.57 65.39 50.00 68.74 106.80 254.74 304.54 222.51 55.71 91.40 186.17 208.90 149.82 29.36 1069.04 191.80 359.74 1019.35 184.26 450.55 939.03 218.18 445.02 122.17 292.10 135.13 229.28 47.23 109.27 265.76 192.00 239.40 68.17 91.45 234.70 166.27 210.69 45.25

主成分分析实例

对于某区域地貌-水文系统，其57个流域盆地的九项地理要素：x1为流域盆地总高度（m)x2为流域盆地山口的海拔高度（m)，x3为流域盆地周长（m)，x4为河道总长度（km)，x5为河

表2-14 某57个流域盆地地理要素数据

道总数，x6为平均分叉率，x7为河谷最大坡度（度)，x8为河源数及x9为流域盆地面积（km2)的原始数据如表2-14所示。张超先生（1984)曾用这些地理要素的原始数据对该区域地貌-水文系统作了主成分分析。下面，我们将其作为主成分分析方法在地理学研究中的一个应用实例介绍给读者，以供参考。

表2-15相关系数矩阵

（1)首先将表2-14中的原始数据作标准化处理，由公式（4)计算得相关系数矩阵（见表2-15)。

（2)由相关系数矩阵计算特征值，以及各个主成分的贡献率与累计贡献率（见表2-16)。由表2-16可知，第一，第二，第三主成分的累计贡献率已高达86.5％，故只需求出第一，第二，第三主成分z1，z2，z3即可。

表2-16 特征值及主成分贡献率

（3)对于特征值λ1=5.043，λ2=1.746，λ3=0.997分别求出其特征向量e1，e2，e3，并计算各变量x1，x2，??，x9在各主成分上的载荷得到主成分载荷矩阵（见表2-17)。

表2-17 主成分载荷矩阵

从表2-17可以看出，第一主成分z与x1，x3，x4，x5，x8，x9有较大的正相关，这是由于这六个地理要素与流域盆地的规模有关，因此第一主成分可以被认为是流域盆地规模的代表：第二主成分z2与x2有较大的正相关，与x7有较大的负相关，而这两个地理要素是与流域切割程度有关的，因此第二主成分可以被认为是流域侵蚀状况的代表；第三主成分z3与x6有较大的正相关，而地理要素x6是流域比较独立的特性——河系形态的表征，因此，第三主成分可以被认为是代表河系形态的主成分。

1

以上分析结果表明，根据主成分载荷，该区域地貌-水文系统的九项地理要素可以被归为三类，即流域盆地的规模，流域侵蚀状况和流域河系形态。如果选取其中相关系数绝对值最大者作为代表，则流域面积，流域盆地出口的海拔高度和分叉率可作为这三类地理要素的代表，利用这三个要素代替原来九个要素进行区域地貌-水文系统分析，可以使问题大大地简化。

利用Matlab编程实现主成分分析

1.概述

Matlab语言是当今国际上科学界 (尤其是自动控制领域) 最具影响力、也是最有

活力的软件。它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、与其他程序和语言的便捷接口的功能。Matlab 语言在各国高校与研究单位起着重大的作用。主成分分析是把原来多个变量划为少数几个综合指标的一种统计分析方法，从数学角度来看，这是一种降维处理技术。

1.1主成分分析计算步骤

① 计算相关系数矩阵

?r11?r21R???????rp1

r12r22?rp2?r1p??r2p???????rpp?? （1）

在（3.5.3）式中，rij（i，j=1，2，…，p）为原变量的xi与xj之间的相关系数，其计算公式为

rij??(xk?1nki?xi)(xkj?xj)2

?(xk?1nki?xi)?(xk?1nkj?xj)2 （2）

因为R是实对称矩阵（即rij=rji），所以只需计算上三角元素或下三角元素即可。

② 计算特征值与特征向量

首先解特征方程?I?R?0，通常用雅可比法（Jacobi）求出特征值

?i(i?1,2,?,p)，并使其按大小顺序排列，即?????,???0；然后分别求出

12p2?1，其中对应于特征值?i的特征向量ei(i?1,2,?,p)。这里要求ei=1，即?eijj?1peij表示向量ei的第j个分量。

③ 计算主成分贡献率及累计贡献率 主成分zi的贡献率为

?i??k?1p(i?1,2,?,p)

k累计贡献率为

????k?1k?1pik(i?1,2,?,p)k

一般取累计贡献率达85—95%的特征值?1,?2,?,?m所对应的第一、第二，…，第m（m≤p）个主成分。

④ 计算主成分载荷 其计算公式为

lij?p(zi,xj)??ieij(i,j?1,2,?,p) （3）

得到各主成分的载荷以后，还可以按照（3.5.2）式进一步计算，得到各主成分的得分

?z11?zZ??21????zn1

z12z22?zn2?z1m??z2m???????znm? （4）

2.程序结构及函数作用

在软件Matlab中实现主成分分析可以采取两种方式实现：一是通过编程来实现；二是直接调用Matlab种自带程序实现。下面主要主要介绍利用Matlab的矩阵计算功能编程实现主成分分析。

2.1程序结构

主函数

子函数

Cwstd.m

Cwprint.m Cwfac.m Cwscore.m 2.2函数作用

Cwstd.m——用总和标准化法标准化矩阵

Cwfac.m——计算相关系数矩阵；计算特征值和特征向量；对主成分进行排序；计算各特征值贡献率；挑选主成分（累计贡献率大于85%），输出主成分个数；计算主成分载荷

Cwscore.m——计算各主成分得分、综合得分并排序 Cwprint.m——读入数据文件；调用以上三个函数并输出结果

3.源程序

3.1 cwstd.m总和标准化法标准化矩阵

%cwstd.m,用总和标准化法标准化矩阵 function std=cwstd(vector)

cwsum=sum(vector,1); %对列求和

[a,b]=size(vector); %矩阵大小,a为行数,b为列数 for i=1:a for j=1:b

std(i,j)= vector(i,j)/cwsum(j); end end

3.2 cwfac.m计算相关系数矩阵

%cwfac.m

function result=cwfac(vector); fprintf('相关系数矩阵:\\n')

std=CORRCOEF(vector) %计算相关系数矩阵 fprintf('特征向量(vec)及特征值(val)：\\n')

[vec,val]=eig(std) %求特征值(val)及特征向量(vec) newval=diag(val) ;

[y,i]=sort(newval) ; %对特征根进行排序，y为排序结果，i为索引 fprintf('特征根排序：\\n') for z=1:length(y)

newy(z)=y(length(y)+1-z); end

fprintf('%g\\n',newy) rate=y/sum(y);

fprintf('\\n贡献率：\\n') newrate=newy/sum(newy) sumrate=0; newi=[];

for k=length(y):-1:1

sumrate=sumrate+rate(k); newi(length(y)+1-k)=i(k); if sumrate>0.85 break;

end

end %记下累积贡献率大85%的特征值的序号放入newi中 fprintf('主成分数：%g\\n\\n',length(newi)); fprintf('主成分载荷：\\n') for p=1:length(newi) for q=1:length(y)

result(q,p)=sqrt(newval(newi(p)))*vec(q,newi(p)); end

end %计算载荷 disp(result)

3.3 cwscore.m

%cwscore.m,计算得分

function score=cwscore(vector1,vector2); sco=vector1*vector2; csum=sum(sco,2);

[newcsum,i]=sort(-1*csum); [newi,j]=sort(i);

fprintf('计算得分：\\n') score=[sco,csum,j]

%得分矩阵：sco为各主成分得分；csum为综合得分；j为排序结果

3.4 cwprint.m

%cwprint.m

function print=cwprint(filename,a,b);

%filename为文本文件文件名，a为矩阵行数(样本数)，b为矩阵列数(变量指标数)

fid=fopen(filename,'r')

vector=fscanf(fid,'%g',[a b]); fprintf('标准化结果如下：\\n') v1=cwstd(vector) result=cwfac(v1); cwscore(v1,result);

4.程序测试

4.1原始数据

中国大陆35个大城市某年的10项社会经济统计指标数据见下表。

城 市 名 称 北 京 天 津 石 家 庄 太 原 呼和浩特 年底 总人口 非农业 农 业 总产值 (万元) 城乡居民客运货运地方财政 年底储蓄总产值 总量 总量 预算内收余额 (万人) (万吨) 入(万元) (万元) (万元) 工业 在岗职在岗职工工人数工资总额 (万人) (万元) 人口比(万人) (%) 1 249.90 0.597 8 1 843 427 19 999 706 20 323 45 562 2 790 863 26 806 646 410.80 5 773 301 910.17 0.580 9 1 501 136 22 645 502 3 259 26 317 1 128 073 11 301 931 202.68 2 254 343 875.40 0.233 2 2 918 680 6 885 768 2 929 1 911 352 348 299.92 0.656 3 236 038 207.78 0.441 2 365 343 2 737 750 1 937 11 895 203 277 816 452 2 351 2 623 105 783 7 095 875 3 943 100 1 396 588 95.60 88.65 42.11 758 877 654 023 309 337 沈 阳 大 连 长 春 哈 尔 滨 上 海 南 京 杭 州 宁 波 合 肥 福 州 厦 门 南 昌 济 南 青 岛 郑 州 武 汉 长 沙 广 州 深 圳 南 宁 海 口 重 庆 成 都 贵 阳 昆 明 西 安 兰 州 西 宁 银 川 乌鲁木齐 677.08 0.629 9 1 295 418 5 826 733 7 782 15 412 567 919 545.31 0.494 6 1 879 739 8 426 385 10 780 19 187 709 227 691.23 0.406 8 1 853 210 5 966 343 4 810 9 532 357 096 927.09 0.462 7 2 663 855 4 186 123 6 720 7 520 481 443 537.44 0.534 1 989 199 13 072 737 14 269 11 193 664 299 616.05 0.355 6 1 414 737 12 000 796 17 883 11 684 449 593 538.41 0.254 7 1 428 235 10 622 866 22 215 10 298 501 723 429.95 0.318 4 628 764 128.99 0.486 5 333 374 424.20 0.398 8 688 289 2 514 125 4 893 1 517 233 628 5 751 124 3 728 2 570 418 758 2 305 881 3 674 3 189 167 714 583.13 0.273 3 2 152 288 6 555 351 8 851 7 190 467 524 9 016 998 135.45 1 152 811 7 556 796 94.15 965 922 884 447 4 803 744 102.63 6 450 020 172.79 1 309 151 5 680 472 113.81 1 357 861 7 425 967 5 246 350 1 622 931 5 030 220 2 108 331 2 640 460 4 126 970 5 135 338 3 461 244 96.90 62.15 47.27 69.59 46.93 62.08 83.31 84.66 69.57 91.26 45.09 19.01 1 180 947 824 034 369 577 680 607 657 484 479 ,555 756 696 961 704 696 848 596 986 1 890 338 371 809 198 138 1 313.12 0.738 4 2 069 019 54 529 098 6 406 44 485 4 318 500 25 971 200 336.84 5 605 445 557.63 0.408 5 1 486 302 6 285 882 5 915 11 775 460 690 702.97 0.369 3 2 382 320 11 492 036 13 408 17 038 658 435 615.36 0.342 4 677 425 5 287 601 10 433 6 768 387 252 740.20 0.586 9 1 211 291 7 506 085 9 793 15 442 604 658 582.47 0.310 7 1 146 367 3 098 179 8 706 5 718 323 660 4 978 045 103.52 5 748 055 149.20 1 314 766 685.00 0.621 4 1 600 738 23 348 139 22 007 23 854 1 761 499 20 401 811 182.81 3 047 594 119.85 0.793 1 299 662 20 368 295 8 754 4 274 1 847 908 9 519 900 285.87 0.406 4 720 486 54.38 0.835 4 44 815 1 149 691 5 130 3 293 149 700 717 461 5 345 2 356 115 174 2 190 918 1 626 800 3 072.34 0.206 7 4 168 780 8 585 525 52 441 25 124 898,912 1 003.56 0.335 1 935 590 5 894 289 40 140 19 632 561 189 321.50 0.455 7 362 061 473.39 0.386 5 793 356 674.50 0.409 4 739 905 287.59 0.544 5 259 444 133.95 0.522 7 95.38 158.92 0.824 4 65 848 78 513 0.570 9 171 603 2 247 934 15 703 4 143 197 908 3 605 729 5 604 12 042 524 216 3 665 942 10 311 9 766 408 896 2 940 884 1 832 4 749 169 540 711 310 661 226 1 746 1 469 2 106 1 193 49 134 74 758 9 090 969 223.73 1 606 804 7 479 684 132.89 1 200 671 1 787 748 4 127 900 2 641 568 855 051 814 103 2 365 508 55.28 88.11 65.83 27.21 23.72 55.27 419 681 842 321 885 169 550 890 219 251 178 621 517 622 5 863 980 114.01 1 847 241 2 668 9 041 254 870

4.2运行结果

>> cwprint('cwbook.txt',35,10)

fid =

6

数据标准化结果如下：

v1 =

0.0581 0.0356 0.0435 0.0680 0.0557 0.1112 0.1194 0.1184 0.1083 0.1392 0.0423 0.0346 0.0354 0.0770 0.0089 0.0642 0.0483 0.0499 0.0534 0.0544

0.0407 0.0139 0.0688 0.0234 0.0080 0.0047 0.0151 0.0314 0.0252 0.0183 0.0139 0.0391 0.0056 0.0093 0.0053 0.0290 0.0087 0.0174 0.0234 0.0158 0.0097 0.0263 0.0086 0.0028 0.0064 0.0064 0.0045 0.0062 0.0111 0.0075 0.0315 0.0375 0.0305 0.0198 0.0213 0.0376 0.0243 0.0398 0.0357 0.0278 0.0253 0.0295 0.0443 0.0286 0.0295 0.0468 0.0304 0.0334 0.0248 0.0233 0.0321 0.0242 0.0437 0.0203 0.0132 0.0233 0.0153 0.0212 0.0270 0.0213 0.0431 0.0276 0.0628 0.0142 0.0184 0.0184 0.0206 0.0285 0.0455 0.0316

0.0610 0.0440 0.0488 0.1853 0.0176 0.1086 0.1848 0.1148 0.0888 0.1352 0.0250 0.0318 0.0233 0.0444 0.0391 0.0273 0.0284 0.0251 0.0300 0.0327 0.0286 0.0212 0.0334 0.0408 0.0490 0.0285 0.0192 0.0328 0.0255 0.0285 0.0250 0.0152 0.0337 0.0361 0.0609 0.0251 0.0215 0.0232 0.0164 0.0199 0.0200 0.0190 0.0148 0.0085 0.0134 0.0037 0.0100 0.0072 0.0125 0.0089 0.0271 0.0163 0.0508 0.0223 0.0243 0.0175 0.0200 0.0222 0.0183 0.0164 0.0060 0.0290 0.0079 0.0195 0.0102 0.0063 0.0179 0.0093 0.0124 0.0159 0.0197 0.0237 0.0162 0.0078 0.0101 0.0078 0.0072 0.0117 0.0164 0.0116 0.0259 0.0243 0.0350 0.0214 0.0162 0.0287 0.0197 0.0182 0.0220 0.0182 0.0327 0.0220 0.0562 0.0391 0.0367 0.0416 0.0282 0.0220 0.0273 0.0232 0.0286 0.0204 0.0160 0.0180 0.0286 0.0165 0.0166 0.0227 0.0223 0.0168 0.0344 0.0349 0.0286 0.0255 0.0268 0.0377 0.0259 0.0254 0.0393 0.0317 0.0271 0.0185 0.0270 0.0105 0.0239 0.0140 0.0139 0.0153 0.0183 0.0144 0.0318 0.0370 0.0377 0.0793 0.0603 0.0582 0.0754 0.0901 0.0482 0.0735 0.0056 0.0472 0.0071 0.0692 0.0240 0.0104 0.0791 0.0421 0.0240 0.0456 0.0133 0.0242 0.0170 0.0039 0.0141 0.0080 0.0064 0.0097 0.0119 0.0090 0.0025 0.0497 0.0011 0.0024 0.0146 0.0057 0.0049 0.0072 0.0050 0.0048 0.1428 0.0123 0.0983 0.0292 0.1437 0.0613 0.0385 0.0402 0.0590 0.0387 0.0466 0.0199 0.0456 0.0200 0.1100 0.0479 0.0240 0.0331 0.0350 0.0290 0.0149 0.0271 0.0085 0.0076 0.0430 0.0101 0.0085 0.0079 0.0146 0.0101 0.0220 0.0230 0.0187 0.0123 0.0154 0.0294 0.0224 0.0182 0.0232 0.0203 0.0313 0.0244 0.0174 0.0125 0.0283 0.0238 0.0175 0.0259 0.0300 0.0213 0.0134 0.0324 0.0061 0.0100 0.0050 0.0116 0.0073 0.0117 0.0173 0.0133 0.0062 0.0311 0.0016 0.0024 0.0048 0.0036 0.0021 0.0038 0.0072 0.0053 0.0044 0.0340 0.0040 0.0022 0.0058 0.0029 0.0032 0.0036 0.0063 0.0043 0.0074 0.0491 0.0019 0.0063 0.0073 0.0221 0.0109 0.0105 0.0146 0.0125

相关系数矩阵:

std =

1.0000 -0.3444 0.8425 0.3603 0.7390 0.6215 0.4039 0.4967 0.6761 0.4689 -0.3444 1.0000 -0.4750 0.3096 -0.3539 0.1971 0.3571 0.2600 0.1570 0.3090 0.8425 -0.4750 1.0000 0.3358 0.5891 0.5056 0.3236 0.4456 0.5575 0.3742 0.3603 0.3096 0.3358 1.0000 0.1507 0.7664 0.9412 0.8480 0.7320 0.8614 0.7390 -0.3539 0.5891 0.1507 1.0000 0.4294 0.1971 0.3182 0.3893 0.2595 0.6215 0.1971 0.5056 0.7664 0.4294 1.0000 0.8316 0.8966 0.9302 0.9027 0.4039 0.3571 0.3236 0.9412 0.1971 0.8316 1.0000 0.9233 0.8376 0.9527 0.4967 0.2600 0.4456 0.8480 0.3182 0.8966 0.9233 1.0000 0.9201 0.9731 0.6761 0.1570 0.5575 0.7320 0.3893 0.9302 0.8376 0.9201 1.0000 0.9396 0.4689 0.3090 0.3742 0.8614 0.2595 0.9027 0.9527 0.9731 0.9396 1.0000

特征向量(vec)： vec =

-0.1367 0.2282 -0.2628 0.1939 0.6371 -0.2163 0.3176 -0.1312 -0.4191 0.2758

-0.0329 -0.0217 0.0009 0.0446 -0.1447 -0.4437 0.4058 -0.5562 0.5487 0.0593

-0.0522 -0.0280 0.2040 -0.0492 -0.5472 -0.4225 0.3440 0.3188 -0.4438 0.2401 0.0067 -0.4176 -0.2856 -0.2389 0.1926 -0.4915 -0.4189 0.2726 0.2065 0.3403 0.0404 0.1408 0.0896 0.0380 -0.1969 -0.0437 -0.4888 -0.6789 -0.4405 0.1861

-0.0343 0.2360 0.0640 -0.8294 0.0377 0.2662 0.1356 -0.1290 0.0278 0.3782 0.2981 0.4739 0.5685 0.2358 0.1465 -0.1502 -0.2631 0.1245 0.2152 0.3644 0.1567 0.3464 -0.6485 0.2489 -0.4043 0.2058 -0.0704 0.0462 0.1214 0.3812 0.4879 -0.5707 0.1217 0.1761 0.0987 0.3550 0.3280 -0.0139 0.0071 0.3832 -0.7894 -0.1628 0.1925 0.2510 -0.0422 0.2694 0.0396 0.0456 0.1668 0.3799

特征值(val) val =

0.0039 0 0 0 0 0 0 0 0 0 0 0.0240 0 0 0 0 0 0 0 0 0 0 0.0307 0 0 0 0 0 0 0 0 0 0 0.0991 0 0 0 0 0 0 0 0 0 0 0.1232 0 0 0 0 0 0 0 0 0 0 0.2566 0 0 0 0 0 0 0 0 0 0 0.3207 0 0 0

0 0 0 0 0 0 0 0.5300 0 0

0 0 0 0 0 0 0 0 2.3514 0 0 0 0 0 0 0 0 0 0 6.2602

特征根排序： 6.26022 2.35138 0.530047 0.320699 0.256639 0.123241 0.0990915 0.0307088 0.0240355 0.00393387

各主成分贡献率：

newrate =

0.6260 0.2351 0.0530 0.0321 0.0257 0.0123 0.0099 0.0031 0.0024 0.0004

第一、二主成分的载荷： 0.690 1 -0.6427 0.148 3 0.8414 0.600 7 -0.6805 0.851 5 0.3167 0.465 6 -0.6754 0.946 3 0.0426 0.911 7 0.3299 0.953 7 0.1862 0.958 9 0.0109 0.950 6 0.2558

第一、二、三、四主成分的得分：

score =

0.718 5 0.049 9 0.768 4 2.0000 0.380 6 0.038 6 0.419 2 4.0000 0.184 8 -0.043 3 0.141 4 21.0000 0.118 6 0.031 1 0.149 7 20.0000 0.054 9 0.011 5 0.066 4 33.0000 0.228 8 0.007 0 0.235 8 7.000 0 0.2364 -0.0081 0.2283 10.0000 0.1778 -0.0167 0.1611 16.0000 0.2292 -0.0337 0.1955 14.0000 0.8382 0.1339 0.9721 1.0000 0.2276 0.0064 0.2340 8.0000 0.2279 -0.0222 0.2056 12.0000 0.1989 -0.0382 0.1607 18.0000 0.0789 -0.0061 0.0728 32.0000 0.1711 -0.0317 0.1394 23.0000 0.0926 0.0266 0.1192 25.0000 0.0900 -0.0000 0.0899 28.0000 0.1692 -0.0082 0.1610 17.0000 0.2441 -0.0318 0.2124 11.0000 0.1507 -0.0108 0.1399 22.0000 0.2316 0.0012 0.2328 9.0000 0.1294 -0.0211 0.1083 27.0000 0.4716 0.0328 0.5045 3.0000 0.2737 0.0834 0.3570 5.0000 0.0754 -0.0013 0.0741 31.0000 0.0448 0.0349 0.0797 30.0000 0.4759 -0.2028 0.2731 6.0000 0.2907 -0.0883 0.2024 13.0000 0.0944 -0.0118 0.0826 29.0000 0.1546 0.0035 0.1581 19.0000 0.1718 -0.0092 0.1626 15.0000 0.0865 0.0230 0.1095 26.0000 0.0349 0.0216 0.0566 35.0000 0.0343 0.0228 0.0572 34.0000 0.0889 0.0422 0.1310 24.0000

主成分分析Matlab源码分析

function [pc, score, latent, tsquare] = princomp(x);

% PRINCOMP Principal Component Analysis (centered and scaled data).

% [PC, SCORE, LATENT, TSQUARE] = PRINCOMP(X) takes a data matrix X and % returns the principal components in PC, the so-called Z-scores in SCORES, % the eigenvalues of the covariance matrix of X in LATENT, and Hotelling's

% T-squared statistic for each data point in TSQUARE.

% Reference: J. Edward Jackson, A User's Guide to Principal Components % John Wiley & Sons, Inc. 1991 pp. 1-25. % B. Jones 3-17-94

% Copyright 1993-2002 The MathWorks, Inc. % $Revision: 2.9 $ $Date: 2002/01/17 21:31:45 $ [m,n] = size(x); % 得到矩阵的规模，m行，n列 r = min(m-1,n); % max possible rank of x % 该矩阵最大的秩不能超过列数， % 也不能超过行数减1

avg = mean(x); % 求每一列的均值，付给一个n维行向量 centerx = (x - avg(ones(m,1),:));

% x的每个元素减去该列的均值， % 使样本点集合重心与坐标原点重合 [U,latent,pc] = svd(centerx./sqrt(m-1),0); % “经济型”的奇异值分解

score = centerx*pc; % 得分矩阵即为原始矩阵乘主成分矩阵 if nargout < 3, return; end

latent = diag(latent).^2; % 将奇异值矩阵转化为一个向量 if (r

latent = [latent(1:r); zeros(n-r,1)]; score(:,r+1:end) = 0; end

if nargout < 4, return; end

tmp = sqrt(diag(1./latent(1:r)))*score(:,1:r)'; tsquare = sum(tmp.*tmp)';

主成分分析[Matlab版] function main()

%*************主成份分析************ %读入文件数据

X=load('data.txt');

%==========方法1:求标准化后的协差矩阵,再求特征根和特征向量================= %标准化处理 [p,n]=size(X); for j=1:n

mju(j)=mean(X(:,j));

sigma(j)=sqrt(cov(X(:,j))); end

for i=1:p for j=1:n

Y(i,j)=(X(i,j)-mju(j))/sigma(j); end end

sigmaY=cov(Y);

%求X标准化的协差矩阵的特征根和特征向量 [T,lambda]=eig(sigmaY);

disp('特征根(由小到大):'); disp(lambda);

disp('特征向量:'); disp(T);

%方差贡献率;累计方差贡献率 Xsum=sum(sum(lambda,2),1); for i=1:n

fai(i)=lambda(i,i)/Xsum; end

for i=1:n

psai(i)= sum(sum(lambda(1:i,1:i),2),1)/Xsum; end

disp('方差贡献率:'); disp(fai);

disp('累计方差贡献率:'); disp(psai);

%综合评价....略

%+============方法2:求X的相关系数矩阵,再求特征根和特征向量================

%X的标准化的协方差矩阵就是X的相关系数矩阵 R=corrcoef(X);

%求X相关系数矩阵的特征根和特征向量 [TR,lambdaR]=eig(R);

disp('特征根(由小到大):'); disp(lambdaR);

disp('特征向量:'); disp(TR);

本文来源：https://www.bwwdw.com/article/llkr.html