主成分分析源代码)
- 格式:doc
- 大小:35.00 KB
- 文档页数:3
R语言进行主成分分析实例1、基于princomp函数进行实例说明:(中学生身体四项指标的主成分分析)在某中学随机抽取某年级30名学生,测量其身高(X1)、体重(X2)、胸围(X3)和坐高(X4),数据如下。
试对这30名中学生身体四项指标数据做主成分分析将上面这些数据保存在students_data.csv中data <- read.csv('D:/students_data.csv', header = T)注:header = T表示将students_data.csv中的第一行数据设置为列名,这种情况下,students_data.csv中的第二行到最后一行数据作为data中的有效数据。
header = F表示不将students_data.csv中的第一行数据设置为列名,这种情况下,students_data.csv 中的第一行到最后一行数据作为data中的有效数据。
第二步:进行主成分分析student.pr <- princomp(data, cor = T)注:cor = T的意思是用相关系数进行主成分分析。
Screeplot(student.pr,type=”line”,main=”碎石图”,lwd=2)第三步:观察主成分分析的详细情况summary(student.pr, loadings = T)执行完这一步的具体结果如下:说明:结果中的Comp.1、Comp.2、Comp.3和Comp.4是计算出来的主成分,Standard deviation代表每个主成分的标准差,Proportion of Variance代表每个主成分的贡献率,Cumulative Proportion代表各个主成分的累积贡献率。
每个主成分都不属于X1、X2、X 3和X4中的任何一个。
第一主成分、第二主成分、第三主成分和第四主成分都是X1、X2、X3和X4的线性组合,也就是说最原始数据的成分经过线性变换得到了各个主成分。
主成分分析类型:一种处理高维数据的方法。
降维思想:在实际问题的研究中,往往会涉及众多有关的变量。
但是,变量太多不但会增加计算的复杂性,而且也会给合理地分析问题和解释问题带来困难。
一般说来,虽然每个变量都提供了一定的信息,但其重要性有所不同,而在很多情况下,变量间有一定的相关性,从而使得这些变量所提供的信息在一定程度上有所重叠。
因而人们希望对这些变量加以“改造”,用为数极少的互补相关的新变量来反映原变量所提供的绝大部分信息,通过对新变量的分析达到解决问题的目的。
一、总体主成分1.1 定义设 X 1,X 2,…,X p 为某实际问题所涉及的 p 个随机变量。
记 X=(X 1,X 2,…,Xp)T ,其协方差矩阵为()[(())(())],T ij p p E X E X X E X σ⨯∑==--它是一个 p 阶非负定矩阵。
设1111112212221122221122Tp p Tp pT pp p p pp p Y l X l X l X l X Y l X l X l X l X Y l X l X l X l X⎧==+++⎪==+++⎪⎨⎪⎪==+++⎩ (1) 则有()(),1,2,...,,(,)(,),1,2,...,.T T i i i i TT T i j ijij Var Y Var l X l l i p Cov Y Y Cov l X l X l l j p ==∑===∑= (2)第 i 个主成分: 一般地,在约束条件1T i i l l =及(,)0,1,2,..., 1.T i k i k Cov Y Y l l k i =∑==-下,求 l i 使 Var(Y i )达到最大,由此 l i 所确定的T i i Y l X =称为 X 1,X 2,…,X p 的第 i 个主成分。
1.2 总体主成分的计算设 ∑是12(,,...,)T p X X X X =的协方差矩阵,∑的特征值及相应的正交单位化特征向量分别为120p λλλ≥≥≥≥及12,,...,,p e e e则 X 的第 i 个主成分为1122,1,2,...,,T i i i i ip p Y e X e X e X e X i p ==+++= (3)此时(),1,2,...,,(,)0,.Ti i i i Ti k i k Var Y e e i p Cov Y Y e e i k λ⎧=∑==⎪⎨=∑=≠⎪⎩ 1.3 总体主成分的性质1.3.1 主成分的协方差矩阵及总方差记 12(,,...,)T p Y Y Y Y = 为主成分向量,则 Y=P T X ,其中12(,,...,)p P e e e =,且12()()(,,...,),T T p Cov Y Cov P X P P Diag λλλ==∑=Λ=由此得主成分的总方差为111()()()()(),p ppTTiii i i i Var Y tr P P tr PP tr Var X λ=====∑=∑=∑=∑∑∑即主成分分析是把 p 个原始变量 X 1,X 2,…,X p 的总方差1()pii Var X =∑分解成 p 个互不相关变量 Y 1,Y 2,…,Y p 的方差之和,即1()pii Var Y =∑而 ()k k Var Y λ=。
主成分分析源程序代码分析function[pc, score, latent, tsquare] = princomp(x);% PRINCOMP Principal Component Analysis (centered and scaled data).% [PC, SCORE, LATENT, TSQUARE] = PRINCOMP(X) takes a data matri x X and% returns the principal components in PC, the so-called Z-scores in SCORE S,% the eigenvalues of the covariance matrix of X in LATENT, and Hotelling 's% T-squared statistic for each data point in TSQUARE.% Reference: J. Edward Jackson, A User's 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 (r<N)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:nmju(j)=mean(X(:,j));sigma(j)=sqrt(cov(X(:,j)));endfor i=1:pfor j=1:nY(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:nfai(i)=lambda(i,i)/Xsum;endfor i=1:npsai(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);。
目录主成分分析和主成分回归(附实际案例和sas代码) (2)1 主成分分析的主要思想 (2)2 主成分分析的定义 (2)3 案例基本情况介绍餐饮业零售额相关因素 (3)4 案例相关因素的介绍相关因素的具体数据 (3)5 影响餐饮业零售额因素的主成分分析 (4)6 主成分回归 (9)主成分分析和主成分回归(附实际案例和sas 代码)1 主成分分析的主要思想在进行高维数据系统分析时,通过主成分分析,可以在纷繁的指标变量描述下,了解影响这个系统存在与发展的主要因素。
主成分分析是1933年由霍特林首先提出来的。
在信息损失最小的前提下,将描述某一系统的多个变量综合成少数几个潜变量,从而迅速揭示系统形成的主要因素,并把原来高维空间降到低维子空间。
主成分分析是研究如何通过少数几个主成分来解释多变量的方差的分析方法,也就是求出少数几个主成分,使他们尽可能多地保留原始变量的信息,且彼此不相关它是一种数学变换方法,即把给定的一组变量通过线性变换,转换为一组不相关的变量,在这种变换中保持变量的总方差不变,同时具有最大总方差,称为第一主成分;具有次大方差,成为第二主成分。
依次类推。
若共有p 个变量,实际应用中一般不是找p 个主成分,而是找出个)(p m m <主成分就够了,只要这m 个主成分能够反映原来所有变量的绝大部分的方差。
2 主成分分析的定义设研究对象涉及P 个指标,分别用p X X X ,,21表示,这个指标构成P 维随机向量为)',,,(21p X X X X =。
设随机向量的均值为u ,协方差矩阵为Σ。
主成分分析就是对随机向量进行线性变换以形成新的综合变量,用i Z 表示,满足下式:1212,1,2,,i i i ip P Z u X u X u X i p =++⋅⋅⋅+= (1)为了使新的综合变量能够充分反映原来变量的信息,则i Z 的方差尽可能大且各个i Z 之间不相关。
由于没有限制条件方差可以任意大,设有线面的约束条件:222121,(1,2,)i i ip u u u i p ++⋅⋅⋅== (2)主成分则为满足条件的i Z 。
Python数据挖掘实战:PCA算法PCA 算法也叫主成分分析(principal components analysis),主要是用于数据降维的。
为什么要进行数据降维?因为实际情况中我们的训练数据会存在特征过多或者是特征累赘的问题,比如:∙一个关于汽车的样本数据,一个特征是”km/h的最大速度特征“,另一个是”英里每小时“的最大速度特征,很显然这两个特征具有很强的相关性∙拿到一个样本,特征非常多,样本缺很少,这样的数据用回归去你和将非常困难,很容易导致过度拟合PCA算法就是用来解决这种问题的,其核心思想就是将n 维特征映射到k维上(k < n),这k 维是全新的正交特征。
我们将这k 维成为主元,是重新构造出来的k 维特征,而不是简单地从n 维特征中取出其余n-k 维特征。
PCA 的计算过程假设我们得到2 维数据如下:其中行代表样例,列代表特征,这里有10个样例,每个样例有2个特征,我们假设这两个特征是具有较强的相关性,需要我们对其进行降维的。
第一步:分别求x 和y 的平均值,然后对所有的样例都减去对应的均值这里求得x 的均值为1.81 ,y 的均值为1.91,减去均值后得到数据如下:注意,此时我们一般应该在对特征进行方差归一化,目的是让每个特征的权重都一样,但是由于我们的数据的值都比较接近,所以归一化这步可以忽略不做第一步的算法步骤如下:本例中步骤3、4没有做。
第二步:求特征协方差矩阵公式如下:第三步:求解协方差矩阵的特征值和特征向量第四步:将特征值从大到小进行排序,选择其中最大的k 个,然后将其对应的k 个特征向量分别作为列向量组成特征矩阵这里的特征值只有两个,我们选择最大的那个,为:1.28402771 ,其对应的特征向量为:注意:matlab 的eig 函数求解协方差矩阵的时候,返回的特征值是一个特征值分布在对角线的对角矩阵,第i 个特征值对应于第i 列的特征向量第五步:将样本点投影到选取的特征向量上假设样本列数为m ,特征数为n ,减去均值后的样本矩阵为DataAdjust(m*n),协方差矩阵为n*n ,选取k 个特征向量组成后的矩阵为EigenVectors(n*k),则投影后的数据FinalData 为:FinalData (m*k)= DataAdjust(m*n) X EigenVectors(n*k)得到的结果是:这样,我们就将n 维特征降成了k 维,这k 维就是原始特征在k 维上的投影。
主成份分析PCA源代码```pythonimport numpy as np#数据中心化mean = np.mean(data, axis=0)data_centered = data - mean#计算协方差矩阵cov_matrix = np.cov(data_centered, rowvar=False)#计算特征值和特征向量eigenvalues, eigenvectors = np.linalg.eig(cov_matrix) #对特征向量按特征值从大到小排序sorted_indices = np.argsort(eigenvalues)[::-1]sorted_eigenvectors = eigenvectors[:, sorted_indices] #转换数据到新的低维空间#计算方差贡献比例total_variance = np.sum(eigenvalues)explained_variance_ratio = eigenvalues / total_variance return transformed_data, explained_variance_ratio```以上是一个简单的PCA源代码实现,下面对代码进行详细解释。
1. 首先导入所需的库,numpy用于数值计算。
3. 数据中心化:计算原始数据的均值mean,然后将数据减去均值得到数据的中心化版本data_centered。
4. 计算协方差矩阵:使用numpy函数cov计算数据的协方差矩阵,设置rowvar=False表示每一列代表一个特征。
5. 计算特征值和特征向量:使用numpy函数linalg.eig计算协方差矩阵的特征值和特征向量。
6. 对特征向量按特征值从大到小排序:使用numpy函数argsort找到特征值从小到大排序的索引,[::-1]表示将索引逆序得到从大到小的排序,然后对特征向量按照这个排序重新排序。
[stata代码模板]主成分分析及因子分析1. 主成分分析黄色字体为自己填写部分,红色字体为可缺省部分。
————————————模板————————————factor 变量名,pc factor(#) covariance means mineigen(#)————————————模板————————————pc代表是主成分分析,如果没有pc,则为因子分析。
factor(#)指定保留因子的个数,可缺省。
covariance指定主成分是从协方差阵计算,而不是从相关阵,也就是说,不加covariance意味着变量被标准化了,可缺省。
means给出各变量的均数、标准差、最小值、最大值,可缺省。
mineigen(#)指定保留的最小特征根。
2. 因子分析主成分分析是将原指标的综合,因子分析是将原指标分解。
(1)因子载荷估计黄色字体为自己填写部分,红色字体为可缺省部分。
————————————模板————————————factor 变量名, factor(#) covariance means 因子提取的方法————————————模板————————————factor(#)、covariance、means与前面意义一样。
因子提取的方法有:Pf 主因子法(缺省时默认)pcf 主成分因子法ipf 迭代因子法ml 极大似然法mineigen(#)指定保留的最小特征根,用主成分提取因子时,缺失值为1,其他情况缺失值为0。
(2)因子旋转当因子估计的模型中的公共因子含义不清或没有合理解释时,可对因子载荷阵进行旋转,使因子载荷的结构简化,以便于对公共因子进行解释。
其原理很像调节显微镜的焦点,以便看清楚观察物的细微之处。
————————————模板————————————rotate,因子旋转的方法————————————模板————————————因子旋转的方法可以缺省,常有以下三种:正交方差极大旋转(varimax),默认为此斜交旋转(promax(#),括号内数为参加旋转的因子数),一般取2或3个因子参加旋转,stata中promax(3)为缺省值。
R语言主成分分析案例Question1Q1.1:> print(eigen_values)[1] 2.4802416 0.9897652 0.3565632 0.1734301Q1.2> print(eigen_vectors)[,1] [,2] [,3] [,4][1,] -0.5358995 0.4181809 -0.3412327 0.64922780[2,] -0.5831836 0.1879856 -0.2681484 -0.74340748[3,] -0.2781909 -0.8728062 -0.3780158 0.13387773[4,] -0.5434321 -0.1673186 0.8177779 0.08902432Q1.3> print('variance for each eigen_values')[1] "variance for each eigen_values"> print(scores)Comp.1 Comp.2 Comp.3 Comp.40.9655342206 0.027******* 0.0057995349 0.0008489079Question2:Q2.1:See in codeQ2.2:The result of ordinary linear regression:> OLSCall:lm(formula = Apps ~ ., data = collegeTrainData)Coefficients:(Intercept) Private Accept Enroll Top10perc Top25perc F.Undergrad-8.753e+02 -6.409e+02 1.345e+00 -2.841e-01 4.792e+01 -1.465e+01 1.980e-02P.Undergrad Outstate Room.Board Books Personal PhD Terminal-1.612e-03 -4.370e-02 2.831e-01 2.356e-01 8.284e-02 1.552e-01 -9.877e+00S.F.Ratio perc.alumni Expend Grad.Rate1.547e+01 -6.582e+00 6.118e-02 4.944e+00And the result in terms of MSE and r-squared is;> print(mse)[1] 1454941> print(rsqured)[1] 0.9162122Q2.3:Use the lambda of seq(0, 1, 0.05) in r, which means from 0 to 1 by 0.05,The result by ridge regression of cross validation is:> print(mse)[1] 1464329> print(ridgeRsquared)[1] 0.9156716Which is slightly worse than the ordinary linear regression.Q2.3:Use the lambda of seq(0, 1, 0.05) in r, which means from 0 to 1 by 0.05,The result by lasso regression of cross validation is:> mse[1] 1471047> LassoRsquared[1] 0.9152847And I make the following table to compare the parameters by the three different models:It can found that Lasso set the parameter of “Phd” to 0. Then it can be inferred that the adjusted r-square of Lasso regression is the best among the three models.Question3:Q3.1:> h_1 = sd(F12)*(4/3/length(F12))^(1/5)> h_1[1] 0.3101212Q3.2:> min(F12)[1] -2.995732> max(F12)[1] 7.930889The min value of log_F12 is -2.99, the maximum value is 7.93. Therefore, I choose the sample from -3 to 8 by 0.05, the following is the plot of the estimated density.Q3.3:I choose 4 different bandwidth:h_2 <- 0.1h_3 <- 0.2h_4 <- 0.5h_5 <- 0.7And the following plot can be get:The middle one is the plot by question b.And the numerical summary of the simulated density for the five different bandwidthWe can see that the larger bandwidth will cause a evener gentler distribution.。
(完整版)主成分分析m a t l a b源程序代码-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN263.862 1.61144 2.75468 0.266575268.764 2.07218 2.61756 0.182597261.196 1.59769 2.35037 0.182114248.708 2.09609 2.85279 0.257724253.365 1.69457 2.9492 0.189702268.434 1.56819 2.78113 0.13252258.741 2.14653 2.69111 0.136469244.192 2.02156 2.22607 0.298066219.738 1.61224 1.88599 0.166298244.702 1.91477 2.25945 0.187569245.286 2.12499 2.35282 0.161602251.96 1.83714 2.53519 0.240271251.164 1.74167 2.62961 0.211887251.824 2.00133 2.62665 0.211991257.68 2.14878 2.65686 0.203846]stdr=std(dataset); %求个变量的标准差[n,m]=size(dataset); %定义矩阵行列数sddata=dataset./stdr(ones(n,1),:); %将原始数据采集标准化sddata %输出标准化数据[p,princ,eigenvalue,t2]=princomp(sddata);%调用前三个主成分系数p3=p(:,1:3); %提取前三个主成分得分系数,通过看行可以看出对应的原始数据的列,每个列在每个主成分的得分p3 %输出前三个主成分得分系数sc=princ(:,1:3); %提取前三个主成分得分值sc %输出前三个主成分得分值e=eigenvalue(1:3)'; %提取前三个特征根并转置M=e(ones(m,1),:).^0.5; %输出前三个特征根并转置compmat=p3.*M; %利用特征根构造变换矩阵per=100*eigenvalue/sum(eigenvalue); %求出成分载荷矩阵的前三列per%求出各主成分的贡献率cumsum(per); %列出各主成分的累积贡献率figure(1)pareto(per); %将贡献率绘成直方图t2figure(2)%输出各省与平局距离plot(eigenvalue,'r+'); %绘制方差贡献散点图hold on%保持图形plot(eigenvalue,'g-'); %绘制方差贡献山麓图%关闭图形plot(princ(:,1),princ(:,2),'+'); %绘制2维成份散点图%gname%,(rowname) %标示个别散点代表的省data市[st2,index]=sort(t2);%st2=flipud(st2);%index=flipud(index);%extreme=index(1);。
主成分分析(principal component analysis)是将多指标化为少数几个综合指标的一种统计分析方法,这种降维的技术而生成的主成分,能够反映原始变量的绝大部分信息,通常表示为原始变量的线性组合。
下面主要介绍在R中的主成分分析(1)概念:①主成分的均值和协方差阵②主成分的总方差贡献率及累计贡献率③原始变量与主成分变量之间的相关系数④m个主成分对原始变量的贡献率⑤原始变量对主成分的影响(2)从相关矩阵或者协方差矩阵出发求主成分①变量的标准化scale()(3)在R中,可以用stats包中的prcomp函数及princmp()函数进行主成分分析。
## 类'formula'的S3方法prcomp(formula, data = NULL, subset, na.action, ...)## Default S3 method:prcomp(x, retx = TRUE, center = TRUE, scale = FALSE,tol = NULL, ...)参数介绍:formula:在公式方法中设定的没有因变量的公式,用来指明数据分析用到的数据框汇中的列data:包含在formula中指定的数据的数据框对象,subset:向量对象,用来指定分析时用到的观测值,其为可选参数na.action:指定处理缺失值的函数x:在默认的方法下,指定用来分析的数值型或者复数矩阵retx:逻辑变量,指定是否返回旋转变量center:逻辑变量,指定是否将变量中心化scale:逻辑变量,指定是否将变量标准化tol:数值型变量,用来指定精度,小于该数值的值将被忽略。
princomp(formula, data = NULL, subset, na.action, ...)## Default S3 method:princomp(x, cor = FALSE, scores = TRUE, covmat = NULL,subset = rep_len(TRUE, nrow(as.matrix(x))), ...)## S3 method for class 'princomp'predict(object, newdata, ...)参数介绍:formula:是没有相应变量的公式cor:逻辑变量,若为cor=T表示用样本的相关矩阵R作主成分分析,cor=F,表示用样本的协方差矩阵s作为主成分分析covmat:协方差矩阵,如果数据不用x提供,可由协方差提供。
1.function y = pca(mixedsig)
2.
3.%程序说明:y = pca(mixedsig),程序中mixedsig为 n*T 阶混合数据矩阵,
n为信号个数,T为采样点数
4.% y为 m*T 阶主分量矩阵。
5.% n是维数,T是样本数。
6.
7.if nargin == 0
8. error('You must supply the mixed data as input argument.');
9.end
10.if length(size(mixedsig))>2
11. error('Input data can not have more than two dimensions. ');
12.end
13.if any(any(isnan(mixedsig)))
14. error('Input data contains NaN''s.');
15.end
16.
17.%——————————————去均值————————————
18.meanValue = mean(mixedsig')';
19.[m,n] = size(mixedsig);
20.%mixedsig = mixedsig - meanValue*ones(1,size(meanValue)); %当数据本
身维数很大时容易出现Out of memory
21.for s = 1:m
22. for t = 1:n
23.mixedsig(s,t) = mixedsig(s,t) - meanValue(s);
24. end
25.end
26.[Dim,NumofSampl] = size(mixedsig);
27.oldDimension = Dim;
28.fprintf('Number of signals: %d\n',Dim);
29.fprintf('Number of samples: %d\n',NumofSampl);
30.fprintf('Calculate PCA...');
31.firstEig = 1;
stEig = Dim;
33.covarianceMatrix = corrcoef(mixedsig'); %计算协方差矩阵
34.[E,D] = eig(covarianceMatrix); %计算协方差矩阵的特征
值和特征向量
35.
36.%———计算协方差矩阵的特征值大于阈值的个数lastEig———
37.%rankTolerance = 1;
38.%maxLastEig = sum(diag(D) >= rankTolerance);
39.%lastEig = maxLastEig;
stEig = 10;
41.
43.eigenvalues = flipud(sort(diag(D)));
44.
45.%—————————去掉较小的特征值——————————
46.if lastEig < oldDimension
47. lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig
+ 1))/2;
48.else
49. lowerLimitValue = eigenvalues(oldDimension) - 1;
50.end
51.lowerColumns = diag(D) > lowerLimitValue;
52.
53.%—————去掉较大的特征值(一般没有这一步)——————
54.if firstEig > 1
55. higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(f
irstEig))/2;
56.else
57. higherLimitValue = eigenvalues(1) + 1;
58.end
59.higherColumns = diag(D) < higherLimitValue;
60.
61.%—————————合并选择的特征值——————————
62.selectedColumns =lowerColumns & higherColumns;
63.
64.%—————————输出处理的结果信息—————————
65.fprintf('Selected [%d] dimensions.\n',sum(selectedColumns));
66.fprintf('Smallest remaining (non-zero) eigenvalue[ %g ]\n',eigenval
ues(lastEig));
67.fprintf('Largest remaining (non-zero) eigenvalue[ %g ]\n',eigenvalu
es(firstEig));
68.fprintf('Sum of removed eigenvalue[ %g ]\n',sum(diag(D) .* (~select
edColumns)));
69.
70.%———————选择相应的特征值和特征向量———————
71.E = selcol(E,selectedColumns);
72.D = selcol(selcol(D,selectedColumns)',selectedColumns);
73.
74.%——————————计算白化矩阵———————————
75.whiteningMatrix = inv(sqrt(D)) * E';
76.dewhiteningMatrix = E * sqrt(D);
77.
78.%——————————提取主分量————————————
79.y = whiteningMatrix * mixedsig;
80.
82.function newMatrix = selcol(oldMatrix,maskVector)
83.if size(maskVector,1)~= size(oldMatrix,2)
84. error('The mask vector and matrix are of uncompatible size.
');
85.end
86.numTaken = 0;
87.for i = 1:size(maskVector,1)
88. if maskVector(i,1) == 1
89.takingMask(1,numTaken + 1) = i;
90.numTaken = numTaken + 1;
91. end
92.end
93.newMatrix = oldMatrix(:,takingMask);。