MATLAB在主成分分析中的应用

（一）引言：

主成分分析作为一种重要的多元统计分析方法，广泛地应用于人口统计学、数量地理学、分子动力学模拟、数学建模、数理分析等学科中，它的优越性主要体现在既可以降低数据“维度”，又保留了原数据的大部分信息，能达到分析、简化数据集的效果。本文仅以一个案例，并利用两种不同的方法讲述如何在MATLAB软件中实现主成分分析过程。

（二）案例重述：

城市生态环境化是城市发展的必然趋势，表现为社会、经济、环境与生态全方位的现代化水平，一个符合生态规律的生态城市应该是结构合理、功能高效和城市生态协调的城市生态系统。所谓结构合理是指适度的人口密度，合理的土地应用，良好的环境质量，充足的绿地系统完善的基础设施，有效的自然保护；功能高效是指资源的优化配置、物力的经济投入、人力的充分发挥、物流的畅通有序、信息流的快捷；关系协调是指人和自然协调、社会关系协调、城乡协调、资源利用和更新协调。一个城市要实现生态城市的发展目标，关键是在市场经济的体制下逐步改善城市的生态环境质量、防止生态环境质量恶化，因此，对城市的生态环境水平调查评价很有必要。

我们对江苏省十个城市的生态环境状况进行了解，得到生态环境指标的指标值，见表1.现对生态环境水平分析和评价。

一级指标	结构				功能			协调度
二级指标	人口结构x1	基础设施x2	地理结构x3	城市绿化x4	物质还原x5	资源配置x6	生产效率x7	城市文明x8	可持续性x9
无锡市	0.7883	0.7663	0.4745	0.8246	0.8791	0.9538	0.8785	0.6305	0.8928
常州市	0.7391	0.7287	0.5126	0.7603	0.8736	0.9257	0.8542	0.6187	0.7831
镇江市	0.8111	0.7629	0.881	0.6888	0.8183	0.9285	0.8537	0.6313	0.5608
张家港市	0.6587	0.8552	0.8903	0.8977	9446	0.9434	0.9027	0.7415	0.8419
连云港市	0.6543	0.7564	0.8288	0.7926	0.9202	0.9154	0.8729	0.6398	0.8464
扬州市	0.8259	0.7455	0.785	0.7856	0.9263	0.8871	0.8485	0.6142	0.7616
泰州市	0.8486	0.78	0.8032	0.6509	0.9185	0.9357	0.8473	0.5734	0.8234
徐州市	0.6834	0.949	0.8862	0.8902	0.9505	0.876	0.9044	0.898	0.6384
南京市	0.8495	0.8918	0.3987	0.6799	0.862	0.9575	0.8866	0.6186	0.9604
苏州市	0.7846	0.8954	0.397	0.9877	0.8873	0.9741	0.9035	0.7382	0.8514

(三)利用princomp函数实现主成分分析过程（方法一）

Step1:将原始数据标准化

x=load(‘data.txt’);

x=x’;

a=zscore(x);

注：zscore函数解释

基于原始数据的均值（mean）和标准差（standard deviation）进行数据的标准化。将原始数据x使用z-score标准化到a。

即：新数据=（原数据-均值）/标准差

调用用法为：

y=zscore(x)

x为标准化之前的数据，y为标准化后的数据

Step2:计算主成分系数矩阵、主成分得分和特征根

[p,princ,egenvalue]=princomp(a)

注：princomp函数解释

1) p是变换矩阵，也是矩阵x所对应的协方差阵的所有特征向量组成的矩阵,通过p可以知道x是怎样转换成princ的。

2) princ是主成分得分，即对原数据进行变换，在新的坐标系下获得的数据，并且按照贡献率从大到小排列。

3) egenvalue是一维列向量矩阵，即x所对应的协差阵的所有特征值，并且按照从大到小的顺序排列。

Step3:计算主成分的方差贡献率

per=100*egenvalue/sum(egenvalue)

输出结果为：

per =[43.1236 29.3409 11.7745 7.5527 4.6069 2.9347 0.5123 0.1540 0.0004]’

Step4：由累计方差贡献率确定主成分个数

由于第一主成分的方差贡献率为43.12%，第二主成分的方差贡献率为29.34%，第三主成分的方差贡献率为11.97%，前三个主成分累计贡献率达84.43%。如果按80%以上的信息量选取因子，则可以选取前三个因子。

Step5：输出前三个主成分系数、主成分、特征根和贡献率

p3=p(:,1:3)

sc=princ(:,1:3)

egenvalue=egenvalue(1:3,:)

per=per(1:3,:)

输出结果如下：

(由于特征向量不唯一，一个特征值对应的所有特征向量是所求向量的任意常数倍，但为与论文保持一致，以后的p3、sc均均取为下面的形式)

p3 =

-0.3677 0.1442 -0.3282

0.3702 0.2313 -0.3535

0.1364 -0.5299 0.0498

0.4048 0.1812 0.0582

0.3355 -0.1601 0.5664

-0.1318 0.5273 -0.0270

0.4236 0.3116 -0.0958

0.4815 -0.0267 -0.2804

-0.0643 0.4589 0.5933

sc =

-0.8301 1.3897 0.3946

-1.3364 -0.2159 0.2729

-1.8408 -1.6267 -2.0552

2.3754 0.1623 0.9043

0.3634 -0.9002 1.4326

-0.9266 -1.6636 0.5837

-1.8984 -0.6907 0.5432

3.9332 -1.5024 -0.9765

-1.2134 2.3442 -0.4498

1.3736 2.7034 -0.6496

egenvalue=[3.8811,2.6407,1.0597]’

per=[43.12,29.34,11.97]’

Step6:对主成分进行加权综合，并得出排序结果

F=0.4312*princ(:,1)+0.2934*princ(:,2)+0.1197*princ(:,3)

[Ft,index]=sort(F,'descend')

（注：使得分按照降序排列，同时输出10个城市的代号）

输出结果如下：

F=[0.0970 -0.6069 -1.5170 1.1801 0.0640 -0.8178 -0.9562 1.1383 0.1107 1.3077]’

Ft=[1.3077 1.1801 1.1383 0.1107 0.0970 0.0640 -0.6069 -0.8178 -0.9562 -1.5170]’

index=[10 4 8 9 1 5 2 6 7 3]’

根据前三个主成分得分，用其贡献率加权，即得十个城市各自的总得分，即F,通过降序排列后得到Ft，则十个城市的排序如下：

城市	排序
苏州市	1
张家港市	2
徐州市	3
南京市	4
无锡市	5
连云港市	6
常州市	7
扬州市	8
泰州市	9
镇江市	10

（三）利用eig函数实现主成分分析过程（方法二）

Step1:将原始数据标准化

x=load(‘data.txt’);

x=x’;

a=zscore(x);

Step2：建立变量的相关系数矩阵

b=cov(a)

或b=corrcoef(x)

注：标准化数据的协差阵正好是原指标的相关阵

Step3：计算特征值、特征向量

[c,d]=eig(b)

e=diag(d)

[dt,index]=sort(e,'descend')

注：

1) c是协差阵b的所有特征向量组成的矩阵

2) d是以协差阵b的所有特征值为主对角线的对角阵

3) 利用diag函数提取d的主对角线元素

4) 利用sort对e进行降序排列

Step4:计算主成分的方差贡献，根据累计方差贡献率确定主成分个数，并写出主成分得分。

per=dt/sum(dt)

n=1;

while sum(dt(1:n))/sum(dt)<0.84

n=n+1;

end

输出结果如下：

per =[43.1236 29.3409 11.7745 7.5527 4.6069 2.9347 0.5123 0.1540 0.0004]’

n=3

因此我们可以看出前三个主成分的累计贡献率达84%以上，则可以提取前三个主成分。

p=c(:,index)

（使特征向量按贡献率从大到小排列）

p3=p(:,1:3)

（提取贡献率最大的三个特征值对应的特征向量）

y=a*p

（计算所有主成分得分）

sc=y(:,1:3)

（提取前三个主成分）

Step5: 对主成分进行加权综合，并得出排序结果

F=0.4312*sc(:,1)+0.2934*sc(:,2)+0.1197*sc(:,3)

（结果同第一种方法）

（四）

附录1：MATLAB程序及运算结果（方法一）

x =

0.7883 0.7633 0.4745 0.8246 0.8791 0.9538 0.8785 0.6305 0.8928

0.7391 0.7287 0.5126 0.7603 0.8736 0.9257 0.8542 0.6187 0.7831

0.8111 0.7629 0.8810 0.6888 0.8183 0.9285 0.8537 0.6313 0.5608

0.6587 0.8552 0.8903 0.8977 0.9446 0.9434 0.9027 0.7415 0.8419

0.6543 0.7564 0.8288 0.7926 0.9202 0.9154 0.8729 0.6398 0.8464

0.8259 0.7455 0.7850 0.7856 0.9263 0.8871 0.8485 0.6142 0.7616

0.8486 0.7800 0.8032 0.6509 0.9185 0.9357 0.8473 0.5734 0.8234

0.6834 0.9490 0.8862 0.8902 0.9505 0.8760 0.9044 0.8980 0.6384

0.8495 0.8918 0.3987 0.6799 0.8620 0.9579 0.8866 0.6186 0.9604

0.7846 0.8954 0.3970 0.9877 0.8873 0.9741 0.9035 0.7382 0.8514

>> a=zscore(x)

a =

0.3158 -0.6381 -0.9962 0.2674 -0.4585 0.7826 0.1394 -0.4146 0.8167

-0.3329 -1.0839 -0.8165 -0.3302 -0.5917 -0.1322 -0.8967 -0.5372 -0.1090

0.6164 -0.6432 0.9209 -0.9947 -1.9304 -0.0410 -0.9180 -0.4063 -1.9851

-1.3930 0.5461 0.9648 0.9467 1.1272 0.4441 1.1713 0.7382 0.3872

-1.4510 -0.7270 0.6747 -0.0300 0.5365 -0.4675 -0.0993 -0.3180 0.4252

0.8115 -0.8674 0.4682 -0.0951 0.6841 -1.3888 -1.1397 -0.5839 -0.2905

1.1108 -0.4229 0.5540 -1.3469 0.4953 0.1934 -1.1909 -1.0076 0.2311

-1.0673 1.7546 0.9455 0.8770 1.2700 -1.7502 1.2438 2.3636 -1.3302

1.1227 1.0176 -1.3537 -1.0774 -0.8725 0.9161 0.4848 -0.5382 1.3872

0.2670 1.0640 -1.3617 1.7832 -0.2600 1.4435 1.2054 0.7040 0.4674

>> [p,princ,egenvalue]=princomp(a)

p =

-0.3677 -0.1442 -0.3282 0.5463 -0.3313 0.4538 0.2368 -0.2554 0.0286

0.3702 -0.2313 -0.3535 0.4979 0.2129 -0.0732 -0.1143 0.6093 0.0392

0.1364 0.5299 0.0498 0.0470 0.5365 0.4650 0.3518 0.0265 0.2595

0.4048 -0.1812 0.0582 -0.3933 -0.5199 0.4445 0.2514 0.3147 0.1229

0.3355 0.1601 0.5664 0.4181 -0.1529 0.2725 -0.3914 -0.1230 -0.3166

-0.1318 -0.5273 -0.0270 -0.2566 0.4320 0.5016 -0.4358 -0.0990 0.0066

0.4236 -0.3116 -0.0958 -0.0084 0.2356 -0.0824 0.4800 -0.3932 -0.5167

0.4815 0.0267 -0.2804 0.0289 -0.0997 -0.1063 -0.2694 -0.5329 0.5569

-0.0643 -0.4589 0.5933 0.2350 0.1061 -0.1760 0.3127 -0.0139 0.4877

princ =

-0.8301 -1.3897 0.3946 -0.5109 -0.3448 0.0021 0.1369 -0.2779 0.0040

-1.3364 0.2159 0.2729 -0.8771 -0.5227 -0.6756 -0.3559 0.0253 -0.0071

-1.8408 1.6267 -2.0552 -0.8162 0.5612 0.2345 0.1149 0.0084 0.0001

2.3754 -0.1623 0.9043 -0.3563 0.8660 0.4842 0.0349 -0.0299 -0.0108

0.3634 0.9002 1.4326 -0.6753 0.4730 -0.4261 0.1340 0.1190 0.0101

-0.9266 1.6636 0.5837 0.6376 -1.0987 0.3040 0.2890 0.0633 -0.0047

-1.8984 0.6907 0.5432 1.1447 0.3917 0.5903 -0.3383 -0.0286 0.0047

3.9332 1.5024 -0.9765 0.7153 -0.2558 -0.4346 -0.0903 -0.0704 0.0028

-1.2134 -2.3442 -0.4498 1.1867 0.5229 -0.6784 0.1563 0.0521 -0.0032

1.3736 -2.7034 -0.6496 -0.4487 -0.5929 0.5997 -0.0814 0.1388 0.0041

egenvalue =

3.8811

2.6407

1.0597

0.6797

0.4146

0.2641

0.0461

0.0139

0.0000

>> per=100*egenvalue/sum(egenvalue)

per =

43.1236

29.3409

11.7745

7.5527

4.6069

2.9347

0.5123

0.1540

0.0004

>> p3=p(:,1:3)

p3 =

-0.3677 -0.1442 -0.3282

0.3702 -0.2313 -0.3535

0.1364 0.5299 0.0498

0.4048 -0.1812 0.0582

0.3355 0.1601 0.5664

-0.1318 -0.5273 -0.0270

0.4236 -0.3116 -0.0958

0.4815 0.0267 -0.2804

-0.0643 -0.4589 0.5933

>> sc=princ(:,1:3)

sc =

-0.8301 -1.3897 0.3946

-1.3364 0.2159 0.2729

-1.8408 1.6267 -2.0552

2.3754 -0.1623 0.9043

0.3634 0.9002 1.4326

-0.9266 1.6636 0.5837

-1.8984 0.6907 0.5432

3.9332 1.5024 -0.9765

-1.2134 -2.3442 -0.4498

1.3736 -2.7034 -0.6496

>> egenvalue=egenvalue(1:3,:)

egenvalue=

3.8811

2.6407

1.0597

>>per=per(1:3,:)

per =

43.1236

29.3409

11.7745

>>F=0.4312*sc(:,1)+0.2934*sc(:,2)+0.1197*sc(:,3)

F =

0.0970

-0.6069

-1.5170

1.1801

0.0640

-0.8178

-0.9562

1.1383

0.1107

1.3077

>>[Ft,index]=sort(F,'descend')

Ft =