因子分析MATLAB程序源代码.doc
11页
clear all;DATA=load('D:0.m');DATA=double(DATA);DATA=DATA';TESTDATA=load('D:14f.m');TESTDATA=double(TESTDATA);% DATA=load('D:正常.txt');% DATA=double(DATA);% DATA=DATA(:,3:12);% TESTDATA=load('D:异常.txt');% TESTDATA=double(TESTDATA);% TESTDATA=TESTDATA(:,3:12);[Kp,T2]=tztq(DATA,TESTDATA);function [contribution,T2,SPE,t2cl,s_cl] = PCA_model(Xtrain,Xtest)X_mean = mean(Xtrain); X_std = std(Xtrain); [X_row ,X_col]= size(Xtrain); for i = 1:X_col Xtrain(:,i) = (Xtrain(:,i)-X_mean(i))./X_std(i); Xtest(:,i) = (Xtest(:,i)-X_mean(i))./X_std(i);end[U,S,V]=svd(Xtrain./sqrt(size(Xtrain,1)-1),0); D= S^2;lamda=diag(D);num_pc=1;while sum(lamda(1:num_pc))/sum(lamda)<0.9 num_pc=num_pc+1;endD=diag(lamda);P=V(:,1:num_pc);[a,b]=size(Xtest);[r,y]=size(P*P');I=eye(r,y);e=Xtest*(I-P*P');for i=1:a T2(i)=Xtest(i,:)*P*inv(D(1:num_pc,1:num_pc))*P'*Xtest(i,:)';endfor l=1:a SPE(l)=e(l,:)*e(l,:)';end for j=1:b contribution(j)=(norm(e(:,j)))^2; end t2cl=num_pc*(X_row-1)*(X_row+1)*icdf('f',0.99,num_pc,X_row-num_pc)/(X_row*(X_row-num_pc));for i=1:3 theta(i)=trace((D(num_pc+1:X_col,num_pc+1:X_col))^i); end% 另一种SPE控制线算法% h=(theta(1)^2)/theta(2);% g=theta(2)/theta(1);% conf=0.95; % df=round(h); % delta2a1=g*pinv(df,2);h0=1-2*theta(1)*theta(3)/(3*theta(2)^2);ca=icdf('norm',0.99,0,1);s_cl=theta(1)*(ca*sqrt(2*theta(2)*h0^2)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0);clear all;X1=load('D:0.m');Xtrain=X1';Xtrain=double(Xtrain);X2=load('D:14f.m');Xtest=X2;Xtest=double(Xtest);% X1=load('D:正常br.txt');% Xtrain=X1(:,3:62);% Xtrain=double(Xtrain);% X2=load('D:异常br.txt');% Xtest=X2(:,3:62);% Xtest=double(Xtest);[contribution,T2,SPE,t2cl,s_cl]=PCA_model(Xtrain,Xtest); figure x=size(Xtest,1); plot(1:x,SPE,'k'); hold on; plot(1:x,s_cl,'r-'); title('SPE'); hold off; figure plot(1:x,T2,'K'); title('T2'); hold on; plot(1:x,t2cl,'r-'); hold off; figure bar(contribution,'group') title('奉献图');function [Kp,T2]=tztq(ax,ay)[Nx]=size(ax);mean_X = mean(ax);axb=ax;std_X=std(ax);ax=ax-mean_X(ones(Nx,1),:);std_X(find(std_X==0))=1;%数据预处理ax=ax./std_X(ones(Nx,1),:);c=10000;% gama=0.05;% ni=1;% F1=ax(1,:);% F=F1';% for ib=2:Nx% for i=1:ni% for j=1:ni% % batar1(ib).block(i,j)=exp(-norm(ax(i,:)-ax(j,:))^2/c);% end% batar2(i,1)=exp(-norm(ax(i,:)-ax(ib,:))^2/c);% batar3(1,i)=exp(-norm(ax(ib,:)-ax(i,:))^2/c);% end% s1=exp(-norm(ax(ib,:)-ax(ib,:))^2/c);% batar= batar3(1,i)*inv(batar1(ib).block)* batar2(i,1);% s=(s1- batar)/s1;% if s>sqrt(gama)% ni=ni+1;% F=[F ax(ib,:)'];% end% % % end% AX=F'%训练集基旳提取结束 [N]=size(ax,1); for i=1:N for j=1:N K(i,j)=exp(-norm(ax(i,:)-ax(j,:))^2/c);%求核矩阵 end end n1=ones(N,N); N1=1/N*n1; Kp=K-N1*K-K*N1+N1*K*N1; [u,s,v]=svd(Kp/N); lamda=s; P=v; lamda=diag(lamda); B=length(find(lamda>1e-10)); %求非零旳特性值个数 %求主元个数npc=1;while sum(lamda(1:npc))/sum(lamda(1:B))<0.9 npc=npc+1;endnpc%求特性空间有效维数DimFS=1;while sum(lamda(1:DimFS))/sum(lamda(1:B))<=0.99 DimFS=DimFS+1;endlamda=diag(lamda);for i=1:B % P(:,i)=P(:,i)/norm(P(:,i)*s(i,i));P(:,i)=P(:,i)/(norm(P(:,i))*sqrt(N*lamda(i,i)));end[Ny]=size(ay,1);mean_X =mean(axb);std_X = std(axb);[num_sample] = Ny; ay = ay-mean_X(ones(num_sample,1),:); ay = ay./std_X(ones(num_sample,1),:); % mean_y = mean(ay);% std_y=std(ay);% ay = ay-mean_y(ones(Ny,1),:);% std_y(find(std_y==0))=1;%数据处理% ay = ay./std_y(ones(Ny,1),:);for i=1:Ny for j=1:N Ky(i,j)=exp(-norm(ay(i,:)-ax(j,:))^2/c); endend t1=ones(1,N); t11=1/N*t1; for i=1:Ny kp1(i,:)= Ky(i,:)-t11*K- Ky(i,:)*N1+t11*K*N1; endfor i=1:Ny for k=1:B t(i,k)=P(:,k)'*kp1(i,:)'; endend % 求T2,SPE % covtyb=inv(t'*t); for i=1:Ny T2(i)=t(i,1:npc)*inv(lamda(1:npc,1:npc))*t(i,1:npc)'; %也可以% SPE(i)=t(i,1:npc)*t(i,1:npc)';% T2(1,i)=t(i,1:npc)*(covtyb(1:npc,1:npc))*t(i,1:npc)'; SPE(i)=t(i,(npc+1):B)*t(i,(npc+1):B)'; end %T2,SPE控制线 t2cl=npc*(N-1)*(N+1)*icdf('f',0.99,npc,N-npc)/(N*(N-npc)); for i=1:3 theta(i)=trace((lamda(npc+1:DimFS,npc+1:DimFS))^i); end h0=1-2*theta(1)*theta(3)/(3*theta(2)^2); ca=icdf('norm',0.99,0,1); s_cl=theta(1)*(ca*sqrt(2*theta(2)*h0^2)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0); figure plot(1:Ny,T2,'k'); hold on; plot(1:Ny,t2cl,'r'); title('T2'); hold off; figure plot(1:Ny。
