因子分析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 [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;
end
D=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,:)';
end
for 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,:)'];
% en

d
%
%
% 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;
end
npc
%求特征空间有效维数
DimFS=1;
while sum(lamda(1:DimFS))/sum(lamda(1:B))<=0.99
DimFS=DimFS+1;
end
lamda=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);
end
end
t1=ones(1,N);
t11=1/N*t1;
for i=1:Ny
kp1(i,:)= Ky(i,:)-t11*K- Ky(i,:)*N1+t11*K*N1;
end
for i=1:Ny
for k=1:B
t(i,k)=P(:,k)'*kp1(i,:)';
end
end
% 求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,SPE,'k')
hold on;
plot(1:Ny, s_cl,'r');
title('SPE');
hold off;