%程序说明：PCA= pca(X)，程序中X为 n*m阶混合数据矩阵，m为信号个数即变量的个数，n为采样点数即样本个数
% T为主分量矩阵,P为得分向量矩阵，Q和T2为统计量。
%function PCA= pca(X)
tic;
load('高维特征.mat')
X=features(:,1:30);
[m n]=size(X);
if X == 0
error('You must supply the mixed data as input argument.');%判断是否输入数据
end
if length(size(X))>2
error('Input data can not have more than two dimensions. ');%判断输入数据维数
end
if any(any(isnan(X)))
error('Input data contains NaN''s.');%判断输入数据是否含有非数值的数据
end
%——————————————标准化数据（首先去均值，然后标定到单位方差）————————————
a=0.99;
meanValue = mean(X); %求采集数据的样本均值
e=ones(size(X,1),1);
s=std(X);
d=diag(s.^-1,0);
fprintf('Normalizing the input data X：\n');
normX = (X - e*meanValue)*d; %将采集数据去均值并标准化
[NumofSampl,Dim] = size(X); %最初采集数据的维数及每个变量样本的个数
oldDim = Dim; %最初变量个数，即数据的维数
fprintf('Number of process variables: %d\n',Dim);
fprintf('Number of samples of every variable: %d\n',NumofSampl);
fprintf('Calculate PCA：\n');
%———计算协方差矩阵的特征值大于阈值的个数numEig(实际求出特征值中非零的特征值个数)———
firstEig=1;
numEig = 0;
covarianceMatrix = (1/(NumofSampl-1))*(normX'*normX); %计算normX的协方差矩阵
[V,D] = eig(covarianceMatrix) ; %计算协方差矩阵的特征值和特征向量,全部特征值构成对角阵D，特征向量为V的列向量
threshold = 1e-5;
for n=1:Dim
if D(n,n)>threshold
numEig=numEig+1;
else if D(n,n)<=threshold
numEig=numEig; %实际求出的特征值大于零的非零特征值的个数
end
end
end
%——————————降序排列各特征值——————————
eigenvalues = flipud(sort(diag(D))); %按从大到小顺序排列的所有的特征值
t=eigenvalues';
%fprintf('The eigenvalues are as follows: %g \n',eigenvalues);
%fprintf('Number of non-zero eigenvalues: %d\n',numEig);
%—————————去掉较小的特征值（去掉较小值的阈值是去掉最大和最小特征值的和求平均）——————————
if numEig < oldDim %特征值的个数小于最初数据的维数的情况
lowerValue=(sum(eigenvalues)-eigenvalues(1)-eigenvalues(2))/numEig;
else
lowerValue=(sum(eigenvalues)-eigenvalues(1)-eigenvalues(numEig))/numEig; %特征值个数等于最初数据维数的情况
end
columns = diag(D) >=lowerValue ;%求得去掉较小特征值后的特征值向量（逻辑向量，大于设定值的特征值为1，小于设定值的为0）
%—————————合并选择的特征值，确定所选择则的特征值——————————
selectedColumns =columns ;
newDim=sum(selectedColumns);
%—————————输出处理的结果信息—————————
fprintf('Selected[ %d ] dimensions.\n',newDim);
fprintf('Smallest remaining (non-zero) eigenvalue: %g \n',eigenvalues(numEig));
fprintf('Largest remaining (non-zero) eigenvalue : %g \n',eigenvalues(firstEig));
%———————选择相应的特征值和特征向量———————
B=V*fliplr(diag(selectedColumns)); %根据特征值大小，选择相应的特征向量
y=find(sum(B)~=0); %计算经过换算过的矩阵的列向量是否为0
selectedEvector=B(:,y); %根据所选择的特征值，确定所选择的单位特征向量组成的新矩阵
%——————————计算负荷向量构成的矩阵———————————
disp('Calculate the loading vectors:');
P=selectedEvector ; %计算负荷向量P
%——————————提取主分量T————————————
disp('Calculate the score vectors:');
T= normX*P ; %计算得分向量，即主分量矩阵
%——————————计算Q统计量，即主元模型的SPE———————————
disp('Calculate the model prediction values:');
normXpre=T*P'; %计算归一化后的主元模型预测值矩阵
error=normX-normXpre; %计算归一化的相对于主元模型的预测的误差
squaredError=error.^2 ;%计算归一化后的平方预测误差
disp('Calculate the Q statistic:');
Q=sum(squaredError,2) ; %计算Q统计量
%——————————计算T统计量———————————
selectedEigvalue=diag(eigenvalues(1:newDim)); %计算所选择的特征值矩阵
disp('Calculate the T2 statistic:');
T2=diag(T*inv(selectedEigvalue)*T'); %计算T统计量
CPV=sum(selectedEigvalue(:))/sum(eigenvalues(:))
%——————————计算控制限Qa和Ta———————————
sparedEigvalue=diag(eigenvalues(newDim+1:Dim)); %提取被删除的特征值
q1=sum(sparedEigvalue); %计算θ1
q1=sum(q1,2);
q2=sum(sparedEigvalue.^2); %计算θ2
q2=sum(q2,2);
q3=sum(sparedEigvalue.^3); %计算θ3
q3=sum(q3,2);
ho=1-(2*q1*q3)/(3*q2^2) ;%计算h0
Ca=icdf('norm',a,0,1); %计算正态分布在a检验水平下的临界值Ca
disp('SPE control limits:');
Qa=q1*(((Ca*(2*q2*ho^2)^0.5)/q1+1+(q2*ho*(ho-1))/q1^2)^(1/ho)); %计算Q统计量
F=icdf('f',a,newDim,NumofSampl-newDim);%计算F分布的值
disp('Statistic T control limits:');
Ta=(newDim*(NumofSampl-1)/(NumofSampl-newDim))*F ;%计算T统计量的控制限Ta
%——————————画出Q与T统计量与控制线的图———————————
Qamatrix=repmat(Qa,NumofSampl,1); %Q控制限转换成矩阵形式
Tamatrix=repmat(Ta,NumofSampl,1); %T控制限转换成矩阵形式
k=0:NumofSampl-1;
H1=figure;
subplot(2,1,1);
plot(k,Q,'b-',k,Qamatrix,'r:');%画图显示Q统计量与控制限
legend('Q','99% Confidence limit');
title('Conventional PCA: Q statistic');
xlabel('k');
ylabel('Q');
H2=figure
subplot(2,1,2);
plot(k,T2,'b-',k,Tamatrix,'r:');%画图显示T统计量与控制限
legend('T2','99% Confidence limit');
title('Conventional PCA: T2 statistic');
xlabel('k');
ylabel('T2');
CPVmatrix=repmat(CPV,NumofSampl,1); %主元累积贡献率转换成矩阵
cpvk=0:1:NumofSampl-1;
H2=figure;
CPVlimits=repmat(0.85,NumofSampl,1);
subplot(2,1,1);
plot(cpvk,CPVmatrix,'b',cpvk,CPVlimits,'r:');%画出主元分析累积贡献率的柱状图
title('Conventional PCA:CPV');
axis([0,NumofSampl,0,1]);
xlabel('k');
ylabel('CPV');
legend('CPV','CPV=85%');
subplot(2,1,2);
plot(cpvk,newDim,'b*');%画出每次移动窗计算的主元个数
title('Conventional PCA: The number of Principal Components(PCs)');
axis([0,NumofSampl+1,0,Dim]);
xlabel('k');
ylabel('The number of PCs');
legend('The number of PCs');
toc %计算程序运行的时间
