主成分分析和MATLAB应用经典版,推荐.doc

function pc, score, latent, tsquare = princomp(x);% PRINCOMP Principal Component Analysis (centered and scaled data).% PC, SCORE, LATENT, TSQUARE = PRINCOMP(X) takes a data matrix X and% returns the principal components in PC, the so-called Z-scores in SCORES,% the eigenvalues of the covariance matrix of X in LATENT, and Hotellings% T-squared statistic for each data point in TSQUARE.% Reference: J. Edward Jackson, A Users Guide to Principal Components% John Wiley & Sons, Inc. 1991 pp. 1-25.% B. Jones 3-17-94% Copyright 1993-2002 The MathWorks, Inc. % $Revision: 2.9 $ $Date: 2002/01/17 21:31:45 $m,n = size(x); % 得到矩阵的规模，m行，n列r = min(m-1,n); % max possible rank of x % 该矩阵最大的秩不能超过列数， % 也不能超过行数减1avg = mean(x); % 求每一列的均值，付给一个n维行向量centerx = (x - avg(ones(m,1),:); % x的每个元素减去该列的均值， % 使样本点集合重心与坐标原点重合U,latent,pc = svd(centerx./sqrt(m-1),0); % “经济型”的奇异值分解score = centerx*pc; % 得分矩阵即为原始矩阵乘主成分矩阵if nargout 3, return; endlatent = diag(latent).2; % 将奇异值矩阵转化为一个向量if (rN) latent = latent(1:r); zeros(n-r,1); score(:,r+1:end) = 0;endif nargout 4, return; endtmp = sqrt(diag(1./latent(1:r)*score(:,1:r);tsquare = sum(tmp.*tmp);主成分分析Matlab版function main()%*主成份分析*%读入文件数据X=load(data.txt);%=方法1:求标准化后的协差矩阵,再求特征根和特征向量=%标准化处理p,n=size(X);for j=1:n mju(j)=mean(X(:,j); sigma(j)=sqrt(cov(X(:,j);endfor i=1:p for j=1:n Y(i,j)=(X(i,j)-mju(j)/sigma(j); endendsigmaY=cov(Y);%求X标准化的协差矩阵的特征根和特征向量T,lambda=eig(sigmaY);disp(特征根(由小到大):);disp(lambda);disp(特征向量:);disp(T);%方差贡献率;累计方差贡献率Xsum=sum(sum(lambda,2),1);for i=1:n fai(i)=lambda(i,i)/Xsum;endfor i=1:n psai(i)= sum(sum(lambda(1:i,1:i),2),1)/Xsum;enddisp(方差贡献率:);disp(fai);disp(累计方差贡献率:);disp(psai);%综合评价.略%+=方法2:求X的相关系数矩阵,再求特征根和特征向量=%X的标准化的协方差矩阵就是X的相关系数矩阵R=corrcoef(X);%求X相关系数矩阵的特征根和特征向量TR,lambdaR=eig(R);disp(特征根(由小到大):);disp(lambdaR);disp(特征向量:);disp(TR);主成分分析类型：一种处理高维数据的方法。降维思想：在实际问题的研究中，往往会涉及众多有关的变量。但是，变量太多不但会增加计算的复杂性，而且也会给合理地分析问题和解释问题带来困难。一般说来，虽然每个变量都提供了一定的信息，但其重要性有所不同，而在很多情况下，变量间有一定的相关性，从而使得这些变量所提供的信息在一定程度上有所重叠。因而人们希望对这些变量加以“改造”，用为数极少的互补相关的新变量来反映原变量所提供的绝大部分信息，通过对新变量的分析达到解决问题的目的。一、总体主成分1.1 定义设 X1，X2，Xp 为某实际问题所涉及的 p 个随机变量。记X=(X1，X2，,Xp)T，其协方差矩阵为 它是一个 p 阶非负定矩阵。设 （1）则有 （2）第 i 个主成分：一般地，在约束条件 及下，求 li 使 Var(Yi)达到最大，由此 li 所确定的 称为 X1，X2，Xp 的第 i 个主成分。1.2 总体主成分的计算设 是的协方差矩阵，的特征值及相应的正交单位化特征向量分别为 及则 X 的第 i 个主成分为 （3）此时 1.3 总体主成分的性质1.3.1 主成分的协方差矩阵及总方差记 为主成分向量，则 Y=PTX，其中，且 由此得主成分的总方差为即主成分分析是把 p 个原始变量 X1，X2，Xp 的总方差分解成 p 个互不相关变量 Y1，Y2，Yp的方差之和，即而 。第 k 个主成分的贡献率：；前m个主成分累计贡献率：，它表明前 m 个主成分Y1，Y2，Ym综合提供 X1，X2，Xp中信息的能力。1.3.2 主成分 Yi 与变量 Xj 的相关系数由于 Y=PTX，故 X=PY，从而由此可得 Yi 与 Xj 的相关系数为 （4）1.4 标准化变量的主成分在实际问题中，不同的变量往往有不同的量纲，由于不同的量纲会引起各变量取值的分散程度差异较大，这时总体方差则主要受方差较大的变量的控制。为了消除由于量纲的不同可能带来的影响，常采用变量标准化的方法，即令 （5）其中 这时 的协方差矩阵便是的相关矩阵 ，其中 （6）利用 X 的相关矩阵 作主成分分析，有如下结论：设 为标准化的随机向量，其协方差矩阵（即 X 的相关矩阵）为 ，则 的第 i 个主成分为 （7）并且 （8）其中 为 的特征值，为相应于特征值 的正交单位特征向量。第 i 个主成分的贡献率：；前 m 个主成分的累计贡献率：；与的相关系数为 。二、样本主成分前面讨论的是总体主成分，但在实际问题中，一般 （或）是未知的，需要通过样本来估计。设为取自 的一个容量为n的简单随机样本，则样本协方差矩阵及样本相关矩阵分别为 （9）其中分别以 S 和 R 作为 和的估计，然后按总体主成分分析的方法作样本主成分分析。三、 例 题某市为了全面分析机械类个企业的经济效益，选择了8个不同的利润指标，14企业关于这8个指标的统计数据如下表所示，试进行主成分分析。表1 14家企业的利润指标的统计数据变量企业序号净产值利润率（%）xi1 固定资产利润率（%）xi2 总产值利润率（%）xi2 销售收入利润率（%）xi3 产品成本利润率(%)xi5物耗利润率(%)xi6人均利润率xi7(千元/人)流动资金利润率(%)xi8140.424.77.26.18.38.72.44220.0225.012.711.211.012.920.23.5429.1313.23.33.94.34.45.50.5783.6422.36.75.63.76.07.40.1767.3534.311.87.17.18.08.91.72627.5635.612.516.416.722.829.33.01726.6722.07.89.910.212.617.60.84710.6848.413.410.99.910.913.91.77217.8940.619.119.819.029.739.62.44935.81024.88.09.88.911.916.20.78913.71112.59.74.24.24.66.50.8743.9121.80.60.70.70.81.10.0561.01332.313.99.48.39.813.32.12617.11438.59.111.39.512.216.41.32711.6解：样本均值向量为：,样本协方差矩阵为：由于S中主对角线元素差异较大，因此我们样本相关矩阵R出发进行主成分分析。样本相关矩阵R为：矩阵R的特征值及相应的特征向量分别为：特征值特征向量6.13660.32113 0.29516 0.38912 0.38472 0.37955 0.37087 0.31996 0.355461.0421-0.4151 -0.59766 0.22974 0.27869 0.31632 0.37151 -0.27814 -0.156840.43595-0.45123 0.10303 -0.039895 0.053874 -0.037292 0.075186 0.77059 -0.424780.22037-0.66817 0.36336 -0.22596 -0.11081 0.14874 0.069353 -0.13495 0.559490.15191-0.038217 0.62435 0.12273 -0.036909 0.15928 0.21062 -0.43006 -0.581050.0088274-0.10167 0.13584 -0.15811 0.86226 -0.25204 -0.34506 -0.13934 -0.0265570.00296240.1596 -0.061134 -0.53966 0.046606 0.7609 -0.27809 0.06203 -0.131260.00122380.19295 -0.031987 -0.64176 0.11002 -0.25397 0.68791 -0.006045 -0.0054031R的特征值及贡献率见下表特征值贡献率（%）累计贡献率（%）6.13660.767080.767081.04210.130270.897340.435950.0544940.951840.220370.0275470.979380.151910.0189880.998370.00882740.00110340.999480.00296240.00037030.999850.00122380.000152971前3个标准化样本主成分类及贡献率已达到95.184%，故只需取前三个主成分即可。前3个标准化样本主成分中各标准化变量 前的系数即为对应特征向量，由此得到3个标准化样本主成分为注意到，y1近似是8个标准化变量的等权重之和，是反映各企业总效应大小的综合指标，y1的值越大，则企业的效益越好。由于y1的贡献率高达76.708%，故若用y1的得分值对各企业进行排序，能从整体上反映企业之间的效应差别。将S中sii的值及中各的值以及各企业关于xi的观测值代入y1的表达式中，可求得各企业y1的得分及其按其得分由大到小的排序结果。企业序号得分12-0.973544-0.648563-0.6274311-0.4855810-0.219497-0.18914-0.00480350.01687980.17711130.1892510.2935120.6531560.8556690.96285所以，第9家企业的效益最好，第12家企业的效益最差。Matlab程序：coeff,score,latent=princomp(X)注：该函数使用协方差阵作主成分分析。主成分分析程序a=; b=corrcoef(zscore(a)%计算相关系数矩阵 D=tril(b)%得到三角矩阵 d,v=eig(b)%计算特征值和特征向量 y1=zscore(a)*d(:,7)%计算第一主成分数值 f1,i1=sort(y1); f2,i2=sort(i1); flipud(i1),flipud(f1),f2%第一主成分得分排序 y2=zscore(a)*d(:,6)%计算第二主成分数值 f1,i1=sort(y2); f2,i2=sort(i1); flipud(i1),flipud(f1),f2%第二主成分得分排序 y3=zscore(a)*d(:,5)%计算第三主成分数值 f1,i1=sort(y3); f2,i2=sort(i1); flipud(i1),flipud(f1),f2%第三主成分得分排序 y4=zscore(a)*d(:,4)%计算第三主成分数值 f1,i1=sort(y4); f2,i2=sort(i1); flipud(i1),flipud(f1),f2%第si主成分得分排序 y5=zscore(a)*d(:,3)%计算第一主成分数值 f1,i1=sort(y1); f2,i2=sort(