《因子分析MATLAB程序源代码.docx》由会员分享，可在线阅读，更多相关《因子分析MATLAB程序源代码.docx（11页珍藏版）》请在三一办公上搜索。
1、因子分析MATLAB程序源代码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 co
2、ntribution,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_colXtrain(:,i) = (Xtrain(:,i)-X_mean(i)./X_std(i);Xtest(:,i) = (Xtest(:,i)-X_mean(i)./X_std(i);endU,S,V=svd(Xtrain./sqrt(size(Xtrain,1)-1),0); D= S2;lamda=diag(D);
3、num_pc=1;while sum(lamda(1:num_pc)/sum(lamda)<0.9num_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:aSPE(l)=e(l,:)*e(l,:);endfor j=1:bcontribution(j)=(norm(e(:,j)2;e
4、ndt2cl=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:3theta(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)
5、2);ca=icdf(norm,0.99,0,1);s_cl=theta(1)*(ca*sqrt(2*theta(2)*h02)/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.tx
6、t);% Xtest=X2(:,3:62);% Xtest=double(Xtest);contribution,T2,SPE,t2cl,s_cl=PCA_model(Xtrain,Xtest);figurex=size(Xtest,1);plot(1:x,SPE,k);hold on;plot(1:x,s_cl,r-);title(SPE);hold off;figureplot(1:x,T2,K);title(T2);hold on;plot(1:x,t2cl,r-);hold off;figurebar(contribution,group)title(贡献图);function Kp,
7、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
8、(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:Nfor j
9、=1:N K(i,j)=exp(-norm(ax(i,:)-ax(j,:)2/c);%求核矩阵 endendn1=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.9npc=npc+1;endnpc%求特征空间有效维数DimFS=1;while sum(lamda(1:DimFS)/
10、sum(lamda(1:B)<=0.99DimFS=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);endNy=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
11、= 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:Nyfor j=1:N Ky(i,j)=exp(-norm(ay(i,:)-ax(j,:)2/c);endendt1=ones(1,N);t11=1/N*t1;for i=1:Nykp1(i,:)= Ky(i,:)-t11*K- Ky(i,:)*N1+t11*K*N1;endfor i=1:Nyfor k=1:Bt(i,k)=P(:,k)*kp1
12、(i,:);endend% 求T2,SPE% covtyb=inv(t*t);for i=1:NyT2(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
13、i=1:3theta(i)=trace(lamda(npc+1:DimFS,npc+1:DimFS)i); endh0=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)*h02)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)2)(1/h0);figureplot(1:Ny,T2,k);hold on;plot(1:Ny,t2cl,r);title(T2);hold off;figureplot(1:Ny,SPE,k)hold on;plot(1:Ny, s_cl,r);title(SPE);hold off;
