close;clc;clear all;
%--------------------------------------------------------------------------
% 载入数据
d1=xlsread('zc.xls');
d2=xlsread('gz.xls');
%--------------------------------------------------------------------------
% 数据的标准化
stantData=d1(:,:); %引入标准数据
c=d2; %引入待检测数据
c1=ones(960,1);
Xmean=mean(stantData);%标准数据的均值
Xstd=std(stantData); %标准数据的方差
X=(c-c1*Xmean)./(c1*Xstd); %用待检测数据减去标准数据的均值并除以方差后，得到标准化的待检测数据
%X=zscore(c); %标准化 %在此之所以用标准数据，是为了知道检测数据的偏离程度，因而减去标准数据的均值和方差
%为什么在此用标准数据的均值和方差
[P, T, LATENT, TSQUARED] = princomp(X); % LATENT是X的特征值；T=T得分向量；P=P负荷矩阵；（即T=XP）
%--------------------------------------------------------------------------
%确定主元素个数k
percent = input('确定方差贡献率限(0～1）：'); % alpha0 是预设的贡献率,通常为85%
k=0;
for i=1:size(LATENT,1); %应用size函数返回LATENT的行维数（语法r=size(A,n)返回n指定的A的方向的维数：n=1为行方向的维数；n=2为列方向的维数）
alpha(i)=sum(LATENT(1:i))/sum(LATENT); %求解贡献率
if alpha(i)>=percent ; %判断是否大于给定贡献率
k=i; %把主元个数提供给K
break;
end;
end;
%--------------------------------------------------------------------------
%存储新的变量
Xp=zeros(size(X));
for j=1:k;
Xp=Xp+T(:,j)*P(:,j)'; %Xp是新的变量
end
%--------------------------------------------------------------------------
beta = input('确定统计阈值置信水平(0～1）：');%0.9 --0.99
theta=zeros(4,1);
for i=1:3 ;
for j=k+1:size(X,2) ; %j=k+1:52
theta(i)=theta(i)+LATENT(j)^(i);
end
end
h0=1-2*theta(1)*theta(3)/[3*(theta(2)^2)];%h0=1-2*theta(1)*theta(3)/3/(theta(2)^2);
SPEbeta=theta(1)*(norminv(beta)*(2*theta(2)*h0^2)^0.5/theta(1)+1+theta(2)*h0*(h0-1)/(theta(1)^2))^(1/h0);%求取SPE的阈值
T2knbeta=finv(beta,k,(size(X,1)-k))*k*(size(X,1)^2-1)/(size(X,1)*(size(X,1)-k)); %利用F分布来求T2统计的阈值
%--------------------------------------------------------------------------
%求取T2的公式
T2=zeros(960,1);
for i=1:960;
T2(i)=X(i,:)*P(:,1:k)*inv(diag(LATENT(1:k)))*P(:,1:k)'*X(i,:)'; % 书上的公式
end
%--------------------------------------------------------------------------
% 求取Q的公式
Q=zeros(960,1);
for i=1:960;
Q(i)=X(i,1:k)*(eye(k)-P(1:k,1:k)*P(1:k,1:k)')'*(eye(k)-P(1:k,1:k)*P(1:k,1:k)')*X(i,1:k)'; % 书上的公式
end
%--------------------------------------------------------------------------
figure(1);
plot(T2(1:960));
hold on;
TS=T2knbeta*ones(960,1); %产生全1矩阵{size(X,1) x 1}维，目的是画一条直线
plot(TS,'r--');
xlabel('采样次数');
ylabel('T2统计图');
hold off;
print -djpeg 144T2; %保存t2统计图像
%--------------------------------------------------------------------------
figure(2);
plot(Q(1:960));
hold on;
SPES=SPEbeta*ones(960,1);%产生全1矩阵{size(X,1) x 1}维，目的是画一条直线
plot(SPES,'r--');
xlabel('采样次数');
ylabel('SPE统计图');
hold off;
print -djpeg 144SPE; %保存SPE统计图像
评论0