MATLAB实现主成分分析

0.利用Matlab编程实现主成分分析1.概述Matlab语言是当今国际上科学界（尤其是自动控制领域）最具影响力、也是最有活力的软件。它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、与其他程序和语言的便捷接口的功能。Matlab语言在各国高校与研究单位起着重大的作用。主成分分析是把原来多个变量划为少数几个综合指标的一种统计分析方法,从数学角度来看，这是一种降维处理技术。1.1主成分分析计算步骤计算相关系数矩阵在（3.5.3）式中,rij工1r12rp2rpp(1)rij(i，j=1为原变量的xi与xj之间的相关系数，其计算公式为nzkdn(XkiXi)(XkjXj)、(Xki-Xi)2(Xkj-Xj)2(2)因为R是实对称矩阵（即k壬rij=rji），所以只需计算上三角元素或下三角元素即可。计算特征值与特征向量首先解特征方程刘一R=0，通常用雅可比法（Jacobi）求出特征值丸i（i二1,2；八，P），并使其按大小顺序排列，即P人X&2二0；然后分别求出对应于特征值人的特征向量e/i=1,2,p）。这里要求|e|=1，即送ej=1，其中e0表示向量e的第j个分量。计算主成分贡献率及累计贡献率主成分z的贡献率为ip二kk4(i=1,2;,p)累计贡献率为iVkk4Pkk4(i=1,2/,p)一般取累计贡献率达8595%的特征值打，入2，扎m所对应的第一、第二，第m（mwp）个主成分。计算主成分载荷其计算公式为lij-p(zi,Xj)=.，iej(i,j二1,2,p)(3)得到各主成分的载荷以后，还可以按照（3.5.2）式进一步计算，得到各主成分的得分Z12Z229Zn2Z1mz2m9znm(4)2.程序结构及函数作用在软件Matlab中实现主成分分析可以采取两种方式实现：一是通过编程来实现；二是直接调用要介绍利用Matlab的矩阵计算功能编程实现主成分分析。2.1程序结构Matlab种自带程序实现。下面主要主2.2函数作用Cwstd.m用总和标准化法标准化矩阵Cwfac.m计算相关系数矩阵；计算特征值和特征向量；对主成分进行排序；计算各特征值贡献率；挑选主成分（累计贡献率大于85%）,输岀主成分个数；计算主成分载荷Cwscore.m计算各主成分得分、综合得分并排序Cwprint.m读入数据文件；调用以上三个函数并输出结果3. 源程序3.1 cwstd.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.2 cwfac.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(%gn,newy)rate=y/sum(y);fprintf(n贡献率：n)newrate=newy/sum(newy)fork=length(y):-1:1sumrate=sumrate+rate(k);newi(length(y)+1-k)=i(k);ifsumrate0.85break;endend%记下累积贡献率大85%的特征值的序号放入newi中fprintf(主成分数：%gnn,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.3 cwscore.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.4 cwprint.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）;4. 程序测试4.1原始数据中国大陆35个大城市某年的10项社会经济统计指标数据见下表。年底非农农业工业客运货运地方财政城乡居民在岗职在岗职工城市名称总人口（万人）业人口比（%）总产值（万元）总产值（万元）总量（万人）总量（万吨）预算内收入（万元）年底储蓄余额（万元）工人数（万人）工资总额（万元）北京1249.900.59781843427199997062032345562279086326806646410.805773301天津910.170.580915013622645502325926317112807311301931202.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.261890338南宁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.02520.01830.01390.03910.00560.00930.00530.02900.00870.01740.02340.01580.00970.02630.00860.00280.00640.00640.00450.00620.01110.00750.03150.03750.03050.01980.02130.03760.02430.03980.03570.02780.02530.02950.04430.02860.02950.04680.03040.03340.02480.02330.03210.02420.04370.02030.01320.02330.01530.02120.02700.02130.04310.02760.06280.01420.01840.01840.02060.02850.04550.03160.06100.04400.04880.18530.01760.10860.18480.11480.08880.13520.02500.03180.02330.04440.03910.02730.02840.02510.03000.03270.02860.02120.03340.04080.04900.02850.01920.03280.02550.02850.02500.01520.03370.03610.06090.02510.02150.02320.01640.01990.02000.01900.01480.00850.01340.00370.01000.00720.01250.00890.02710.01630.05080.02230.02430.01750.02000.02220.01830.01640.00600.02900.00790.01950.01020.00630.01790.00930.01240.01590.01970.02370.01620.00780.01010.00780.00720.01170.01640.01160.02590.02430.03500.02140.01620.02870.01970.01820.02200.01820.03270.02200.05620.03910.03670.04160.02820.02200.02730.02320.02860.02040.01600.01800.02860.01650.01660.02270.02230.01680.03440.03490.02860.02550.02680.03770.02590.02540.03930.03170.02710.01850.02700.01050.02390.01400.01390.01530.01830.01440.03180.03700.03770.07930.06030.05820.07540.09010.04820.07350.00560.04720.00710.06920.02400.01040.07910.04210.02400.04560.01330.02420.01700.00390.01410.00800.00640.00970.01190.00900.00250.04970.00110.00240.01460.00570.00490.00720.00500.00480.14280.01230.09830.02920.14370.06130.03850.04020.05900.03870.04660.01990.04560.02000.11000.04790.02400.03310.03500.02900.01490.02710.00850.00760.04300.01010.00850.00790.01460.01010.02200.02300.01870.01230.01540.02940.02240.01820.02320.02030.03130.02440.01740.01250.02830.02380.01750.02590.03000.02130.01340.03240.00610.01000.00500.01160.00730.01170.01730.01330.00620.03110.00160.00240.00480.00360.00210.00380.00720.00530.00440.03400.00400.00220.00580.00290.00320.00360.00630.00430.00740.04910.00190.00630.00730.02210.01090.01050.01460.0125相关系数矩阵std=1.0000-0.34440.84250.36030.73900.62150.40390.49670.67610.4689-0.34441.0000-0.47500.3096-0.35390.19710.35710.26000.15700.30900.8425-0.47501.00000.33580.58910.50560.32360.44560.55750.37420.36030.30960.33581.00000.15070.76640.94120.84800.73200.86140.7390-0.35390.58910.15071.00000.42940.19710.31820.38930.25950.62150.19710.50560.76640.42941.00000.83160.89660.93020.90270.40390.35710.32360.94120.19710.83161.00000.92330.83760.95270.49670.26000.44560.84800.31820.89660.92331.00000.92010.97310.67610.15700.55750.73200.38930.93020.83760.92011.00000.93960.46890.30900.37420.86140.25950.90270.95270.97310.93961.0000特征向量（vec）:vec=-0.13670.2282-0.26280.19390.6371-0.21630.3176-0.1312-0.41910.2758-0.0329-0.02170.00090.0446-0.1447-0.44370.4058-0.55620.54870.0593-0.0522-0.02800.2040-0.0492-0.5472-0.42250.34400.3188-0.44380.24010.0067-0.4176-0.2856-0.23890.1926-0.4915-0.41890.27260.20650.34030.04040.14080.08960.0380-0.1969-0.0437-0.4888-0.6789-0.44050.1861-0.03430.23600.0640-0.82940.03770.26620.1356-0.12900.02780.37820.29810.47390.56850.23580.1465-0.1502-0.26310.12450.21520.36440.15670.3464-0.64850.2489-0.40430.2058-0.07040.04620.12140.38120.4879-0.57070.12170.17610.09870.35500.3280-0.01390.00710.3832-0.7894-0.16280.19250.2510-0.04220.26940.03960.04560.16680.3799特征值（val）val=0.003900000000000.024000000000000.030700000000000.099100000000000.123200000000000.256600000000000.320700000000000.530000000000002.351400000000006.2602特征根排序：6.260222.351380.5300470.3206990.2566390.1232410.09909150.03070880.02403550.00393387各主成分贡献率:newrate=0.62600.23510.05300.03210.02570.01230.00990.00310.00240.0004第一、二主成分的载荷:0.69010.14830.60070.85150.46560.94630.91170.9537-0.64270.8414-0.68050.3167-0.67540.04260.32990.18620.95890.01090.95060.2558第一、二、三、四主成分的得分:score=0.71850.04990.76842.00000.38060.03860.41924.00000.1848-0.04330.141421.00000.11860.03110.149720.00000.05490.01150.066433.00000.22880.00700.23587.00000.2364-0.00810.228310.00000.1778-0.01670.161116.00000.2292-0.03370.195514.00000.83820.13390.97211.00000.22760.00640.23408.00000.2279-0.02220.205612.00000.1989-0.03820.160718.00000.0789-0.00610.072832.00000.1711-0.03170.139423.00000.09260.02660.119225.00000.0900-0.00000.089928.00000.1692-0.00820.161017.00000.2441-0.03180.212411.00000.1507-0.01080.139922.00000.23160.00120.23289.00000.1294-0.02110.108327.00000.47160.03280.50453.00000.27370.08340.35705.00000.0754-0.00130.074131.00000.04480.03490.079730.00000.4759-0.20280.27316.00000.2907-0.08830.202413.00000.0944-0.01180.082629.00000.15460.00350.158119.00000.1718-0.00920.162615.00000.08650.02300.109526.00000.03490.02160.056635.00000.03430.02280.057234.00000.08890.04220.131024.00001249.90,0.5978,1843427,19999706,20323,45562,2790863,26806646,410.80,5773301158.92, 910.17,0.5809,1501136,22645502,3259,26317,1128073,11301931,202.68,2254343875.40,0.2332,2918680,6885768,2929,1911,352348,7095875,95.60,7588770.6563,236038,2737750,1937,11895,203277,3943100,88.65,654023207.78,0.4412,365343,816452,2351,2623,105783,1396588,42.11,309337677.08,0.6299,1295418,5826733,7782,15412,567919,9016998,135.45,1152811545.31,0.4946,1879739,8426385,10780,19187,709227,7556796,94.15,965922691.23,0.4068,1853210,5966343,4810,9532,357096,4803744,102.63,884447927.09,0.4627,2663855,4186123,6720,7520,481443,6450020,172.79,13091511313.12,0.7384,2069019,54529098,6406,44485,4318500,25971200,336.84,5605445537.44,0.5341,989199,13072737,14269,11193,664299,5680472,113.81,1357861616.05,0.3556,1414737,12000796,17883,11684,449593,7425967,96.90,1180947538.41,0.2547,1428235,10622866,22215,10298,501723,5246350,62.15,8240340.3184,628764,2514125,4893,1517,233628,1622931,47.27,369577583.13,0.2733,2152288,6555351,8851,7190,467524,5030220,69.59,680607128.99,0.4865,333374,5751124,3728,2570,418758,2108331,46.93,6574840.3988,688289,2305881,3674,3189,167714,2640460,62.08,479,555557.63,0.4085,1486302,6285882,5915,11775,460690,4126970,83.31,756696702.97,0.3693,2382320,11492036,13408,17038,658435,4978045,103.52,961704615.36,0.3424,677425,5287601,10433,6768,387252,5135338,84.66,6968480.5869,1211291,7506085,9793,15442,604658,5748055,149.20,1314766582.47,0.3107,1146367,3098179,8706,5718,323660,3461244,69.57,596986685.00,0.6214,1600738,23348139,22007,23854,1761499,20401811,182.81,3047594119.85,0.7931,299662,20368295,8754,4274,1847908,9519900,91.26,1890338285.87,0.4064,720486,1149691,5130,3293,149700,2190918,45.09,3718090.8354,44815,717461,5345,2356,115174,1626800,19.01,1981383072.34,0.2067,4168780,8585525,52441,25124,898,912,9090969,223.73,16068041003.56,0.335,1935590,5894289,40140,19632,561189,7479684,132.89,12006710.4557,362061,2247934,15703,4143,197908,1787748,55.28,419681473.39,0.3865,793356,3605729,5604,12042,524216,4127900,88.11,8423210.4094,739905,3665942,10311,9766,408896,5863980,114.01,885169287.59,0.5445,259444,2940884,1832,4749,169540,2641568,65.83,5508900.5227,65848,711310,1746,1469,49134,855051,27.21,2192510.5709,171603,661226,2106,1193,74758,814103,23.72,1786210.8244,78513,1847241,2668,9041,254870,2365508,55.27,517622