注册地址 登录
数学建模社区-数学中国 返回首页

酌水知源的个人空间 http://www.madio.net/?992353 [收藏] [复制] [分享] [RSS]

日志

主成分分析

已有 425 次阅读2014-4-23 22:28

利用Matlab编程实现主成分分析1.概述 

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

 

     计算相关系数矩阵

                                     (1

在(3.5.3)式中,rijij=12p为原变量的xixj之间的相关系数,其计算公式为

                             (2

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

 计算特征值与特征向量

首先解特征方程,通常用雅可比法(Jacobi)求出特征值,并使其按大小顺序排列,即;然后分别求出对应于特征值的特征向量。这里要求=1,即,其中表示向量的第j个分量。

 计算主成分贡献率及累计贡献率

主成分的贡献率为

累计贡献率为

一般取累计贡献率达8595%的特征值所对应的第一、第二,,第mmp个主成分。

 计算主成分载荷

其计算公式为

                       (3

得到各主成分的载荷以后,还可以按照(3.5.2)式进一步计算,得到各主成分的得分

                                    (4

2.程序结构及函数作用

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

2.1程序结构

 

                                                       主函数

 

                    子函数

 

                 

 

 

 

 

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

7 095 875

95.60

758 877

    

299.92

0.656 3

236 038

2 737 750

1 937

11 895

203 277

3 943 100

88.65

654 023

呼和浩特

207.78

0.441 2

365 343

816 452

2 351

2 623

105 783

1 396 588

42.11

309 337

    

677.08

0.629 9

1 295 418

5 826 733

7 782

15 412

567 919

9 016 998

135.45

1 152 811

    

545.31

0.494 6

1 879 739

8 426 385

10 780

19 187

709 227

7 556 796

94.15

965 922

    

691.23

0.406 8

1 853 210

5 966 343

4 810

9 532

357 096

4 803 744

102.63

884 447

  

927.09

0.462 7

2 663 855

4 186 123

6 720

7 520

481 443

6 450 020

172.79

1 309 151

    

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

    

537.44

0.534 1

989 199

13 072 737

14 269

11 193

664 299

5 680 472

113.81

1 357 861

    

616.05

0.355 6

1 414 737

12 000 796

17 883

11 684

449 593

7 425 967

96.90

1 180 947

    

538.41

0.254 7

1 428 235

10 622 866

22 215

10 298

501 723

5 246 350

62.15

824 034

    

429.95

0.318 4

628 764

2 514 125

4 893

1 517

233 628

1 622 931

47.27

369 577

    

583.13

0.273 3

2 152 288

6 555 351

8 851

7 190

467 524

5 030 220

69.59

680 607

    

128.99

0.486 5

333 374

5 751 124

3 728

2 570

418 758

2 108 331

46.93

657 484

    

424.20

0.398 8

688 289

2 305 881

3 674

3 189

167 714

2 640 460

62.08

479 ,555

    

557.63

0.408 5

1 486 302

6 285 882

5 915

11 775

460 690

4 126 970

83.31

756 696

    

702.97

0.369 3

2 382 320

11 492 036

13 408

17 038

658 435

4 978 045

103.52

961 704

    

615.36

0.342 4

677 425

5 287 601

10 433

6 768

387 252

5 135 338

84.66

696 848

    

740.20

0.586 9

1 211 291

7 506 085

9 793

15 442

604 658

5 748 055

149.20

1 314 766

    

582.47

0.310 7

1 146 367

3 098 179

8 706

5 718

323 660

3 461 244

69.57

596 986

广    

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

91.26

1 890 338

    

285.87

0.406 4

720 486

1 149 691

5 130

3 293

149 700

2 190 918

45.09

371 809

    

54.38

0.835 4

44 815

717 461

5 345

2 356

115 174

1 626 800

19.01

198 138

    

3 072.34

0.206 7

4 168 780

8 585 525

52 441

25 124

898,912

9 090 969

223.73

1 606 804

    

1 003.56

0.335

1 935 590

5 894 289

40 140

19 632

561 189

7 479 684

132.89

1 200 671

    

321.50

0.455 7

362 061

2 247 934

15 703

4 143

197 908

1 787 748

55.28

419 681

    

473.39

0.386 5

793 356

3 605 729

5 604

12 042

524 216

4 127 900

88.11

842 321

西    

674.50

0.409 4

739 905

3 665 942

10 311

9 766

408 896

5 863 980

114.01

885 169

    

287.59

0.544 5

259 444

2 940 884

1 832

4 749

169 540

2 641 568

65.83

550 890

西    

133.95

0.522 7

65 848

711 310

1 746

1 469

49 134

855 051

27.21

219 251

    

95.38

0.570 9

171 603

661 226

2 106

1 193

74 758

814 103

23.72

178 621

乌鲁木齐

158.92

0.824 4

78 513

1 847 241

2 668

9 041

254 870

2 365 508

55.27

517 622

 

 

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.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

                                                             


路过

雷人

握手

鲜花

鸡蛋

全部作者的其他最新日志

评论 (0 个评论)

facelist doodle 涂鸦板

您需要登录后才可以评论 登录 | 注册地址

qq
收缩
  • 电话咨询

  • 04714969085

关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

手机版|Archiver| |繁體中文 手机客户端  

蒙公网安备 15010502000194号

Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

GMT+8, 2025-8-19 11:06 , Processed in 0.235767 second(s), 27 queries .

回顶部