因子分析MATLAB程序源代码

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=tztq(DATA,TESTDATA);function contribution,T2,

2、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= SA2;lamda=diag(D);num_pc=1;while su

3、m(lamda(1:num_pc)/sum(lamda)sqrt(gama)%ni=ni+1;%F=F ax(ib,:);%end% AX=F%训练集基的提取结束N=size(ax,1);for i=1:Nfor j=1:NK(i,j)=exp(-norm(ax(i,:)-ax(j,:)A2/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(lamda1e-10); % 求非零的特征值个数%求主元个数np

4、c=1;while sum(lamda(1:npc)/sum(lamda(1:B)0.9npc=npc+1;endnpc%求特征空间有效维数DimFS=1;while sum(lamda(1:DimFS)/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

5、_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:Nyfor j=1:NKy(i,j)=exp(-norm(ay(i,:)-ax(j,:)A2/c);endendt1=ones(1,N);t11=1/N*

6、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(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+

7、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:3theta (i)=trace(lamda(npc+1:DimFS,npc+1:DimFS)Ai);endh0=1-2*theta(1)*theta(3)/(3*theta(2)A2);ca=icdf(norm,0.99,0,1);s_cl=theta(1)*(ca*sqrt(2*theta(2)*h0A2)/theta (1)+1+theta(2)*h0*(h0-1)/theta ( 1)a2)a(1/h0);figureplot(1:Ny,T2,k);ho