MatLAB在主成分分析中的应用

58§10.利用Matlab编程实现主成分分析1.概述Matlab语言是当今国际上科学界(尤其是自动控制领域)最具影响力、也是最有活力的软件。它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、与其他程序和语言的便捷接口的功能。Matlab语言在各国高校与研究单位起着重大的作用。主成分分析是把原来多个变量划为少数几个综合指标的一种统计分析方法，从数学角度来看，这是一种降维处理技术。1.1主成分分析计算步骤①计算相关系数矩阵pppppprrrrrrrrrR212222111211（1）在（3.5.3）式中，rij（i，j=1，2，…，p）为原变量的xi与xj之间的相关系数，其计算公式为nknkjkjikinkjkjikiijxxxxxxxxr11221)()())((（2）因为R是实对称矩阵（即rij=rji），所以只需计算上三角元素或下三角元素即可。59②计算特征值与特征向量首先解特征方程0RI，通常用雅可比法（Jacobi）求出特征值),,2,1(pii，并使其按大小顺序排列，即0,21p；然后分别求出对应于特征值i的特征向量),,2,1(piei。这里要求ie=1，即112pjije，其中ije表示向量ie的第j个分量。③计算主成分贡献率及累计贡献率主成分iz的贡献率为),,2,1(1pipkki累计贡献率为),,2,1(11pipkkikk一般取累计贡献率达85—95%的特征值m,,,21所对应的第一、第二，…，第m（m≤p）个主成分。④计算主成分载荷其计算公式为),,2,1,(),(pjiexzplijijiij（3）60得到各主成分的载荷以后，还可以按照（3.5.2）式进一步计算，得到各主成分的得分nmnnmmzzzzzzzzzZ212222111211（4）2.程序结构及函数作用在软件Matlab中实现主成分分析可以采取两种方式实现：一是通过编程来实现；二是直接调用Matlab种自带程序实现。下面主要主要介绍利用Matlab的矩阵计算功能编程实现主成分分析。2.1程序结构主函数子函数2.2函数作用Cwstd.m——用总和标准化法标准化矩阵Cwfac.m——计算相关系数矩阵；计算特征值和特征向量；对主成分进行排序；计算各特征值贡献率；挑选主成分（累计贡献率大于85%），输出主成分个数；计算主成分载荷Cwprint.mCwstd.mCwfac.mCwscore.m61Cwscore.m——计算各主成分得分、综合得分并排序Cwprint.m——读入数据文件；调用以上三个函数并输出结果3.源程序3.1cwstd.m%cwstd.m,用总和标准化法标准化矩阵functionstd=cwstd(vector)cwsum=sum(vector,1);%对列求和[a,b]=size(vector);%矩阵大小,a为行数,b为列数fori=1:aforj=1:bstd(i,j)=vector(i,j)/cwsum(j);endend3.2cwfac.m%cwfac.mfunctionresult=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')forz=1:length(y)newy(z)=y(length(y)+1-z);endfprintf('%g\n',newy)rate=y/sum(y);fprintf('\n贡献率：\n')newrate=newy/sum(newy)62sumrate=0;newi=[];fork=length(y):-1:1sumrate=sumrate+rate(k);newi(length(y)+1-k)=i(k);ifsumrate0.85break;endend%记下累积贡献率大85%的特征值的序号放入newi中fprintf('主成分数：%g\n\n',length(newi));fprintf('主成分载荷：\n')forp=1:length(newi)forq=1:length(y)result(q,p)=sqrt(newval(newi(p)))*vec(q,newi(p));endend%计算载荷disp(result)3.3cwscore.m%cwscore.m,计算得分functionscore=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.4cwprint.m%cwprint.mfunctionprint=cwprint(filename,a,b);%filename为文本文件文件名，a为矩阵行数(样本数)，b为矩阵列数(变量指标数)fid=fopen(filename,'r')vector=fscanf(fid,'%g',[ab]);fprintf('标准化结果如下：\n')v1=cwstd(vector)result=cwfac(v1);cwscore(v1,result);634.程序测试4.1原始数据中国大陆35个大城市某年的10项社会经济统计指标数据见下表。城市名称年底总人口(万人)非农业人口比(%)农业总产值(万元)工业总产值(万元)客运总量(万人)货运总量(万吨)地方财政预算内收入(万元)城乡居民年底储蓄余额(万元)在岗职工人数(万人)在岗职工工资总额(万元)北京1249.900.59781843427199997062032345562279086326806646410.805773301天津910.170.5809150113622645502325926317112807311301931202.682254343石家庄875.400.23322918680688576829291911352348709587595.60758877太原299.920.65632360382737750193711895203277394310088.65654023呼和浩特207.780.441236534381645223512623105783139658842.11309337沈阳677.080.6299129541858267337782154125679199016998135.451152811大连545.310.4946187973984263851078019187709227755679694.15965922长春691.230.406818532105966343481095323570964803744102.63884447哈尔滨927.090.462726638554186123672075204814436450020172.791309151上海1313.120.7384206901954529098640644485431850025971200336.845605445南京537.440.53419891991307273714269111936642995680472113.811357861杭州616.050.35561414737120007961788311684449593742596796.901180947宁波538.410.25471428235106228662221510298501723524635062.15824034合肥429.950.3184628764251412548931517233628162293147.27369577福州583.130.27332152288655535188517190467524503022069.59680607厦门128.990.4865333374575112437282570418758210833146.93657484南昌424.200.3988688289230588136743189167714264046062.08479,555济南557.630.408514863026285882591511775460690412697083.31756696青岛702.970.369323823201149203613408170386584354978045103.52961704郑州615.360.34246774255287601104336768387252513533884.66696848武汉740.200.5869121129175060859793154426046585748055149.201314766长沙582.470.31071146367309817987065718323660346124469.57596986广州685.000.62141600738233481392200723854176149920401811182.813047594深圳119.850.793129966220368295875442741847908951990091.26189033864南宁285.870.4064720486114969151303293149700219091845.09371809海口54.380.83544481571746153452356115174162680019.01198138重庆3072.340.2067416878085855255244125124898,9129090969223.731606804成都1003.560.3351935590589428940140196325611897479684132.891200671贵阳321.500.45573620612247934157034143197908178774855.28419681昆明473.390.38657933563605729560412042524216412790088.11842321西安674.500.409473990536659421031197664088965863980114.01885169兰州287.590.5445259444294088418324749169540264156865.83550890西宁133.950.522765848711310174614694913485505127.21219251银川95.380.5709171603661226210611937475881410323.72178621乌鲁木齐158.920.824478513184724126689041254870236550855.275176224.2运行结果cwprint('cwbook.txt',35,10)fid=6数据标准化结果如下：v1=0.05810.03560.04350.06800.05570.11120.11940.11840.10830.13920.04230.03460.03540.07700.00890.06420.04830.04990.05340.05440.04070.01390.06880.02340.00800.00470.01510.03140