MATLAB
结课作业
指导老师：张肃 班级：信管121 姓名：桂亚东
学号：201200654118
58
利用Matlab编程实现主成分分析
概述
Matlab语言是当今国际上科学界 (尤其是自动控制领域) 最具影响力、也是
最有活力的软件。它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、与其他程序和语言的便捷接口的功能。Matlab 语言在各国高校与研究单位起着重大的作用。主成分分析是把原来多个变量划为少数几个综合指标的一种统计分析方法，从数学角度来看，这是一种降维处理技术。
1.1主成分分析计算步骤
① 计算相关系数矩阵
?r11?r21R???????rp1r12r22?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），所以只需计算上三角元素或下三角元素即可。
② 计算特征值与特征向量
59
首先解特征方程?I?R?0，通常用雅可比法（Jacobi）求出特征值
?i(i?1,2,?,p)，并使其按大小顺序排列，即?????,???0；然后分别求
12p2出对应于特征值?i的特征向量ei(i?1,2,?,p)。这里要求ei=1，即?eij?1，其
j?1p中eij表示向量ei的第j个分量。
③ 计算主成分贡献率及累计贡献率 主成分zi的贡献率为
?i??k?1p(i?1,2,?,p)
k累计贡献率为
i????k?1k?1pk(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）式进一步计算，得到各主成分的得分
60
?z11?zZ??21????zn1
z12z22?zn2?z1m??z2m???????znm? （4）
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计算相关系数矩阵
61
%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);
62