PCA代码

PCA算法介绍

    PCA算法作为多元统计方法的一种，是通过线性空间变换求取主元变量，将高维数据空间投影到低维主元空间，从而消除观测数据之间的冗余信息，得到主元模型和统计控制限。新的映射空间由原始数据变量的线性组合构成。由于投影空间统计特征向量彼此正交，则消除了变量间的关联性，简化了原始过程数据特性分析的复杂程度。内容包括主元的定义和获取，以及通过主元的数据重构。
    PCA算法适用于线性、高斯分布的数据，实现数据的降维。只涉及二阶统计特性，并没有考虑到数据高阶统计特性，因此变换后的数据可能仍有高阶冗余信息，只解除了数据之间的相关性，对非线性问题并没有进行相应分析。同时，也没有考虑到数据类之间的信息，只是对数据进行重构而不是分类。因此，对PCA算法的改进一直以来都被广为探索。
    PCA假定变量服从高斯分布，则在此基础上计算得到的监控指标T2统计量和Q统计量分布服从一定的分布规律，可以很容易的利用确定的分布特性获得控制置信限。对非高斯分布的变量，不能按某种特定的分布规律推导出控制线。因此，对PCA算法的改进一直以来都被广为探索。
    下面简单介绍PCAS算法在TE过程故障诊断的应用，代码如下：

PCA主程序

函数调用

function [contribution,contribution1,T,SPE,T1,SPE1,T2beta,SPEbeta,cout_T,cout_SPE,feature,feature1,feature2,feature3,PC,PC1]=PCA(train,test)

对训练数据的标准化处理（均值为0，方差为1）

tr=zscore(train);

测试数据的标准化处理

meantr=repmat(mean(train),size(test,1),1) ;%repmat(A,m,n)是将矩阵A扩展，行数*m倍，列数*n倍
stdtr=repmat(std(train),size(test,1),1);
test=(test-meantr)./stdtr;

协方差矩阵的特征值分解

[U,S,V]=svd(tr'*tr./(size(tr,1)-1));
 D=S;
% D= S^2;%X协方差矩阵的特征值为X本身矩阵的奇异值的平方，故取D=S^2
lamda=diag(D);

确定主元个数

numpc=1;
while sum(lamda(1:numpc))/sum(lamda)<0.85   
    numpc=numpc+1;
end

得分空间与残差空间

D=lamda(1:numpc);%主特征值
D=diag(D);%主特征值对应的对角矩阵
INVD=inv(D);
PC=V(:,1:numpc);%PC是主元特征向量；num_pc是主元个数
PC1=V(:,numpc+1:end);%PC1是残差特征向量；num_pc是主元个数
feature=tr*PC1; 
feature1=test*PC1;
feature2=tr*PC;%TR是tr对应的得分矩阵
feature3=test*PC;

计算训练数据的T(i)和SPE(i)

[r,p]=size(PC*PC');
I=eye(r,p);
[a,b]=size(tr);
[a1,b1]=size(test);
e=tr*(I-PC*PC');
for i=1:a;
    T(i)=tr(i,:)*PC*INVD*PC'*tr(i,:)';  
    SPE(i)=e(i,:)*e(i,:)';
end
% SPEE=tr*(I-PC*PC')*tr';

计算控制限T2和控制限SPE

beta=0.99;
T2beta=numpc*(a-1)*(a+1)*finv(beta,numpc,a-numpc)/(a*(a-numpc)); 
theta=zeros(3,1); 
for i=1:3     
    for j=numpc+1:b 
        theta(i)=theta(i)+lamda(j)^(i); 
    end 
end  
h0=1-2*theta(1)*theta(3)/3/(theta(2)^2); 
SPEbeta=theta(1)*(norminv(beta)*sqrt(2*theta(2)*h0^2)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0); 
m = mean(SPE); %%m是建模数据集中所有测量数据的SPE统计量的均值
v = var(SPE);%%v是对应的方差
%SPEbeta = v/(2*m)*chi2inv(0.99,2*m.^2/v);%SPE的控制限

计算测试数据的T（i）和SPE(i)

e1=test*(I-PC*PC');
cout_T=0;
cout_SPE=0;
for i=1:a1;
  T1(i)=test(i,:)*PC*INVD*PC'*test(i,:)';
  if i>=160 && T1(i)>T2beta
     cout_T=cout_T+1;
  end
  SPE1(i)=e1(i,:)*e1(i,:)';
  if i>=160 && SPE1(i)>SPEbeta
     cout_SPE=cout_SPE+1;
  end
end
cout_T=cout_T/(a1-159);
cout_SPE=cout_SPE/(a1-159);

SPE+T2的综合指标

zonghe=1/T2beta*T+1/SPEbeta*SPE;
zh= ksdensity(zonghe,0.99,'function','icdf');
zonghe1=1/T2beta*T1+1/SPEbeta*SPE1;
cout_zh=0;
cout_zhF=0;
for i=1:a1;
  if i<160 && zonghe1(i)>zh
     cout_zh=cout_zh+1;
  end 
  if i>=160 && zonghe1(i)>zh
     cout_zhF=cout_zhF+1;
  end
end
% count_m_zh=1-length(find(zonghe>zh))/960
cout_hh=cout_zhF/(a1-159)*100;
cout_hh_F=cout_zh/159*100;

求测试数据的贡献率

contribution=zeros(1,b);
contribution1=zeros(1,b1);
 for i=1:b;
    contribution(i)=(norm(e(:,i)))^2;% 即：对e1的每一列求平方和，即为贡献值；
 end
  for i=1:b1;
    contribution1(i)=(norm(e1(:,i)))^2;% 即：对e1的每一列求平方和，即为贡献值；
  end