主成分分析源代码

主成分分析也称主分量分析，旨在利用降维的思想，把多指标转化为少数几个综合指标。在实证问题研究中，为了全面、系统地分析问题，我们必须考虑众多影响因素。这些涉及的因素一般称为指标，在多元统计分析中也称为变量。因为每个变量都在不同程度上反映了所研究问题的某些信息，并且指标之间彼此有一定的相关性，因而所得的统计数据反映的信息在一定程度上有重叠。在用统计方法研究多变量问题时，变量太 多会增加计算量和增加分析问题的复杂性，人们希望在进行定量分析的过程中，涉及的变量较少，得到的信息量较多。主成分分析正是适应这一要求产生的。
主成分分析法是一种数学变换的方法, 它把给定的一组相关变量通过线性变换转成另一组不相关的变量，这些新的变量按照方差依次递减的顺序排列。在数学变换中保持变量的总方差不变，使第一变量具有最大的方差，称为第一主成分，第二变量的方差次大，并且和第一变量不相关，称为第二主成分。依次类推，I个变量就有I个主成分。
其中Li为p维正交化向量（Li＊Li＝1），Zi之间互不相关且按照方差由大到小排列，则称Zi为X的第I个主成分。设X的协方差矩阵为Σ，则Σ必为半正定对称矩阵，求特征值λi（按从大到小排序）及其特征向量，可以证明，λi所对应的正交化特征向量，即为第I个主成分Zi所数据标准化;

分析步骤

　求相关系数矩阵;
一系列正交变换，使非对角线上的数置0，加到主对角上;
得特征根xi（即相应那个主成分引起变异的方差),并按照从大到小的顺序把特征根排列;
求各个特征根对应的特征向量;
用下式计算每个特征根的贡献率Vi;
Vi=xi/(x1+x2+........)
根据特征根及其特征向量解释主成分物理意义。对应的系数向量Li，而Zi的方差贡献率定义为λi／Σλj,通常要求提取的主成分的数量k满足Σλk／Σλj>0.85。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = pca(mixedsig)

%程序说明：y = pca(mixedsig)，程序中mixedsig为 n*T 阶混合数据矩阵，n为信号个数，T为采样点数
% y为 m*T 阶主分量矩阵。
% n是维数，T是样本数。

if nargin == 0
error('You must supply the mixed data as input argument.');
end
if length(size(mixedsig))>2
error('Input data can not have more than two dimensions. ');
end
if any(any(isnan(mixedsig)))
error('Input data contains NaN''s.');
end

%——————————————去均值————————————
meanValue = mean(mixedsig')';
[m,n] = size(mixedsig);
%mixedsig = mixedsig - meanValue*ones(1,size(meanValue)); %当数据本身维数很大时

容易出现Out of memory
for s = 1:m
for t = 1:n
mixedsig(s,t) = mixedsig(s,t) - meanValue(s);
end
end
[Dim,NumofSampl] = size(mixedsig);
oldDimension = Dim;
fprintf('Number of signals: %d ',Dim);
fprintf('Number of samples: %d ',NumofSampl);
fprintf('Calculate PCA...');
firstEig = 1;
lastEig = Dim;
covarianceMatrix = corrcoef(mixedsig'); %计算协方差矩阵
[E,D] = eig(covarianceMatrix); %计算协方差矩阵的特征值和特征向量

%———计算协方差矩阵的特征值大于阈值的个数lastEig———
%rankTolerance = 1;
%maxLastEig = sum(diag(D) >= rankTolerance);
%lastEig = maxLastEig;
lastEig = 10;

%——————————降序排列特征值——————————
eigenvalues = flipud(sort(diag(D)));

%—————————去掉较小的特征值——————————
if lastEig < oldDimension
lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig + 1))/2;
else
lowerLimitValue = eigenvalues(oldDimension) - 1;
end
lowerColumns = diag(D) > lowerLimitValue;

%—————去掉较大的特征值(一般没有这一步)——————
if firstEig > 1
higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig))/2;
else
higherLimitValue = eigenvalues(1) + 1;
end
higherColumns = diag(D) < higherLimitValue;

%—————————合并选择的特征值——————————
selectedColumns =lowerColumns & higherColumns;

%—————————输出处理的结果信息—————————
fprintf('Selected [%d] dimensions. ',sum(selectedColumns));
fprintf('Smallest remaining (non-zero) eigenvalue[ %g ] ',eigenvalues(lastEig));
fprintf('Largest remaining (non-zero) eigenvalue[ %g ] ',eigenvalues(firstEig));
fprintf('Sum of removed eigenvalue[ %g ] ',sum(diag(D) .* (~selectedColumns)));

%———————选择相应的特征值和特征向量———————
E = selcol(E,selectedColumns);
D = selcol(selcol(D,selectedColumns)',selectedColumns);

%——————————计算白化矩阵———————————
whiteningMatrix = inv(sqrt(D)) * E';
dewhiteningMatrix = E * sqrt(D);

%——————————提取主分量————————————
y = whiteningMatrix * mixedsig;

%——————————行选择子程序———————————
function newMatrix = selcol(oldMatrix,maskVector)
if size(maskVector,1)~= size(oldMatrix,2)
error('The mask vector and matrix are of uncompatible size.');
end
numTaken = 0;
for i = 1:size(maskVector,1)
if maskVector(i,1) == 1
takingMask(1,numTaken + 1) = i;
numTaken = numTaken + 1;
end
end
newMatrix = oldMatrix(:,takingMask);