function [T2]=partiald(ax,ay) %ax训练数据集，ay测试数据集
[Nx]=size(ax,1);
mean_X = mean(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; %此参数可调
[N]=size(ax,1);
for i=1:N
for j=1:N
K(i,j)=exp(-(norm(ax(i,:)-ax(j,:)))^2/c);%求核矩阵，采用径向基核函数，参数c
end
end
Summationkx=sum(K);
Summationkx1=Summationkx';
Summationkx2=sum(Summationkx1);
n1=ones(N,N);
N1=(1/N)*n1;
Kp=K-N1*K-K*N1+N1*K*N1;
[u,s,v]=svd(Kp/N);
v=-1*v;
lamda=s;
P=v;
lamda=diag(lamda);
%求主元个数
npc=1;
while sum(lamda(1:npc))/sum(lamda)<0.9
npc=npc+1;
end
lamda=diag(lamda);
for i=1:npc
P(:,i)=P(:,i)/(norm(P(:,i))*sqrt(N*s(i,i)));
end
[Ny]=size(ay,1);
[num_sample] = size(ay,1);
ay = ay-mean_X(ones(num_sample,1),:);
ay = ay./std_X(ones(num_sample,1),:);
for i=1:Ny
for j=1:N
Ky(i,j)=exp(-(norm(ay(i,:)-ax(j,:)))^2/c);
end
end
Summationky=sum(Ky')';
etest=ones(N,1);
etest=(1/N)*etest;
for i=1:Ny
kp1(:,i)= Ky(i,:)'-K*etest-N1*Ky(i,:)'+N1*K*etest;
end
for i=1:Ny
for j=1:npc
tnew(j)=P(:,j)'*kp1(:,i);
end
tnew1(i,:)=tnew;
end
for i=1:Ny
T2(i)=tnew1(i,:)*inv(lamda(1:npc,1:npc))*tnew1(i,:)';
SPE(i)=1-2/N*Summationky(i)+1/(N^2)*Summationkx2-tnew1(i,:)*tnew1(i,:)';
end
clT2=npc*(N-1)/(N-npc)*finv(0.99,npc,N-npc); %T2控制线
S=var(SPE);
miu=mean(SPE);
p=round(2*miu^2/S);
cl=(S/(2*miu)*chi2inv(0.99,p)); %Q控制线
figure;
plot(T2,'r');
title('TSQUARE');
hold on;
plot(1:Ny,clT2,'g');
figure;
plot(SPE,'r');
hold on;
title('SPE');
plot(1:Ny,cl,'g');
