因子分析MATLAB程序源代码[1]

1、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=tz

2、tq(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);endU,S,V=svd(Xtrain./sqrt(size(Xtrain

3、,1)-1),0); D= S2;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,

4、:)*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; % d

5、f=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)*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

6、;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,&

7、#39;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(

8、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

9、,:)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=on

10、es(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=d

11、iag(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 = mean(ay);% std_y=std(ay);% ay = ay-mean_y(o

12、nes(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 % c

13、ovtyb=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)*h02)/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)2