KPCAA.zip_pca

%正常模型建立 Xtrain=load('d05_te.txt'); %取训练数据阵的行和列数 [XX_row,XX_col]=size(Xtrain); %平均值 Xt_mean = mean(Xtrain); %方差 Xt_std = std(Xtrain); %训练集矩阵每一列减去平均值除以方差 for i = 1:XX_col Xtrain(:,i) = (Xtrain(:,i)-Xt_mean(i))./Xt_std(i); end %计算训练数据的核矩阵K for i=1:XX_row for j=1:XX_row kX(i,j)=exp(-(norm(Xtrain(i,:)-Xtrain(j,:))).^2/480); end end lmX=1/XX_row*ones(XX_row);%定义训练数据的IN %对训练数据的核矩阵K均值中心化处理 klX=kX-lmX*kX-kX*lmX+lmX*kX*lmX; % trac=trace(klX)/XX_row; % klX=klX/trac; %此函数是利用协方差矩阵对其进行主元分析。pk是主成分（特征向量构成的矩阵），lamd是核矩阵的协方差矩阵的特征值（特征值构成的向量。列向量），baifen是每个特征向量表征在观测量总方差中所占的百分数（向量。列向量）。 [pkX,lamdX,baifenX]=pcacov(cov(klX)); %对训练数据的特征向量进行归一化 for i=1:XX_row psX(:,i)=pkX(:,i)/sqrt(lamdX(i)); end %选择主元个数 num_pc=1;%主元个数的初值为1 while sum(lamdX(1:num_pc))/sum(lamdX)<0.9 num_pc=num_pc+1; end sigmaXtrain = cov(Xtrain); [T,lamda] = eig(sigmaXtrain); D = flipud(diag(lamda)); P = T(:,XX_col-num_pc+1:XX_col); TT=Xtrain*T; TT1=Xtrain*P; %求置信度分别在99%的T2统计量的T2控制限 t2a99=num_pc*(XX_row^2-1)*finv(0.99,num_pc,XX_row-num_pc)/(XX_row*(XX_row-num_pc)); tX=klX*psX;%计算SPE时要用到 for i=1:XX_row spe(i)=sum(tX(i,num_pc+1:XX_col).^2); end a=mean(spe); b=var(spe); Qa99=b/(2*a)*chi2inv(0.999,2*a^2/b); %% 在线监测 Xtest =load('d06_te.txt');%载入检测数据 [XC_row,XC_col]=size(Xtest);%取测试数据阵的行和列数 % 用原始数据的均值和标准差对检测数据进行标准化 Xte_mean=repmat(Xt_mean,XC_row,1); Xte_std=repmat(Xt_std,XC_row,1); for i=1:XC_row for j=1:XX_col Xtest(i,j)=(Xtest(i,j)-Xte_mean(i,j))/Xte_std(i,j); end end %计算测试数据的核矩阵K for i=1:XC_row for j=1:XC_row kC(i,j)=exp(-(norm(Xtest(i,:)-Xtrain(j,:))).^2/480); end end lmC=1/XC_row*ones(XC_row);%定义测试数据的IN klC=kC-lmC*kX-kC*lmC+lmC*kX*lmX;%对测试数据的核矩阵K均值中心化处理 % klC=klC/trac; tC=klC*psX; %T2统计量 lamdX=lamdX(1:num_pc,:); D=diag(lamdX); for i=1:XC_row % T2(i)=tC(i,1:num_pc)*inv(D(1:num_pc,1:num_pc))*tC(i,1:num_pc)'; T2(i)=Xtest(i,:)*P*pinv(lamda(XX_col-num_pc+1:XX_col,XX_col-num_pc+1:XX_col))*P'*Xtest(i,:)'; end for i=1:XC_row Q(i)=sum(tC(i,num_pc+1:XC_col).^2); end %% 出图 %绘图 figure subplot(2,1,1); plot(1:XC_row,T2,'k'); title('对于故障1，用于检测的KPCA多变量统计量'); xlabel('采样数'); ylabel('T^2'); hold on; line([0,XC_row],[t2a99,t2a99],'LineStyle','--','Color','r'); subplot(2,1,2); plot(1:XC_row,Q,'k'); xlabel('采样数'); ylabel('SPE'); hold on; line([0,XC_row],[Qa99,Qa99],'LineStyle','--','Color','r'); %% 贡献图 for i = 1:3 theta(i) = sum((D(num_pc+1:XX_col)).^i); end h0 = 1 - 2*theta(1)*theta(3)/(3*theta(1)^2); ca = norminv(0.95,0,1); QUCL = theta(1)*(h0*ca*sqrt(2*theta(2))/theta(1) + 1 + theta(2)*h0*(h0 - 1)/theta(1)^2)^(1/h0); n = size(Xtest,1); [r,y] = size(P*P'); I = eye(r,y); Q = zeros(n,1); %1.确定造成失控状态的得分 S = Xtest(160,:)*P(:,1:num_pc); r = [ ]; for i = 1:num_pc if S(i)^2/lamda(i) > t2a99/num_pc r = cat(2,r,i); end end %2.计算每个变量相对于上述失控得分的贡献 cont = zeros(length(r),XC_col); for i = length(r) for j = 1:XC_col cont(i,j) = abs(S(i)/D(i)*P(j,i)*Xtest(160,j)); end end %3.计算每个变量的总贡献 CONTJ = zeros(XC_col,1); for j = 1:XC_col CONTJ(j) = sum(cont(:,j)); end %4.计算每个变量对Q的贡献 e = Xtest(160,:)*(I - P*P'); contq = e.^2; %置信度为95%的Q统计控制限 for i = 1:3 theta(i) = sum((D(num_pc+1:XX_col)).^i); end %5. 绘制贡献图 figure; subplot(2,1,1); bar(CONTJ,'k') xlabel('变量号'); ylabel('T^2贡献率 %'); subplot(2,1,2); bar(contq,'k'); xlabel('变量号'); ylabel('Q贡献率 %'); %计算控制限 alpha=0.9; S=lamda(XC_col-num_pc+1:XC_col,XC_col-num_pc+1:XC_col); FAI=P*pinv(S)*P'/t2a99+(eye(XC_col)-P*P')/QUCL; S=cov(Xtrain); g=trace((S*FAI)^2)/trace(S*FAI); h=(trace(S*FAI))^2/trace((S*FAI)^2); kesi =g*chi2inv(alpha,h); %% 综合指标 figure; fai=(Q/QUCL)+(T2/t2a99); plot(fai) hold on; line([0,n],[kesi,kesi],'LineStyle','--','Color','r'); title('混合指标'); %% 可视化 % 原始得分矩阵可视化 figure; subplot(2,3,1) plot(TT(:,1),TT(:,2),'*'); xlabel('PC1'); ylabel('PC2'); title('SCORE PLOT FOR PC1 AND PC2'); subplot(2,3,2) plot(TT(:,1),TT(:,3),'*'); xlabel('PC1'); ylabel('PC3'); title('SCORE PLOT FOR PC1 AND PC3'); subplot(2,3,3) plot(TT(:,1),TT(:,4),'*'); xlabel('PC1'); ylabel('PC4'); title('SCORE PLOT FOR PC1 AND PC4'); subplot(2,3,4) plot(TT(:,2),TT(:,3),'*'); xlabel('PC2'); ylabel('PC3'); title('SCORE PLOT FOR PC2 AND PC3'); subplot(2,3,5) plot(TT(:,2),TT(:,4),'*'); xlabel('PC2'); ylabel('PC4'); title('SCORE PLOT FOR PC2 AND PC4'); subplot(2,3,6) plot(TT(:,3),TT(:,4),'*'); xlabel('PC3'); ylabel('PC4'); title('SCORE PLOT FOR PC3 AND PC4'); %% 主元得分可视化 figure; subplot(2,2,1) plot(TT1(:,1),TT1(:,2),'*'); xlabel('PC1'); ylabel('PC2'); title('SCORE PLOT FOR PC1 AND PC2'); subplot(2,2,2) plot(TT1(:,1),TT1(:,3),'*'); xlabel('PC1'); ylabel('PC3'); title('SCORE PLOT FOR PC1 AND PC3'); subplot(2,2,3) plot(TT1(:,2),TT1(:,3),'*'); xlabel('PC2'); ylabel('PC3'); title('SCORE PLOT FOR PC2 AND PC3'); subplot(2,2,4) plot3(TT1(:,1),TT1(:,2),TT1(:,3),'*'); xlabel('PC1'); ylabel('PC2'); zlabel('PC3'); grid on; title('得分矩阵');