PCA完整程序（MATLAB）

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,:)']; % 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=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;