simca.rar_SIMCA聚类分析_matlab simca_simca matlab_simca怎么聚类_单类SIMCA

function simca_model=simca(x,classes) % % Model for Soft Independent Modeling of Class Analogy (SIMCA) % % simca_model=simca(x,classes) % % input: % x (samples x descriptors) for calibration % classes (samples x 1) classes numbers (must be >0) % % output: % simca_model struct with: % models (1 x No. classes) struct with: (one for each class) % x (samples x descriptors) input for this class % clas (samples x 1) numbers for this class % S (samples x pc) scores for this model % L (descriptors x pc) loadings for this model % ev (pc x 1) eigenvalues for this model % T2 (samples x 1) T?for samples of this model % Q2 (samples x 1) Q?for samples of this model % T2lim (1 x 1) T?limit for this model % Q2lim (1 x 1) Q?limit for this model % stats (1 x 1) struct with: % classes (No. classes x 1) class types % T2 (samples x 1) T?for calibration samples % Q2 (samples x 1) Q?for calibration samples % class_c (samples x 1) Predicted classes for calibration samples % RMSEP (1 x 1) Root Mean Square Error for Calibration % R2 (1 x 1) Correlation Coefficient for calibration % suc (1 x 1) Success (%) of classification for calibration samples % % See also simcapred % % By Cleiton A. Nunes % UFLA,MG,Brazil % %---classes--- y=classes; yu=unique(y); for i=1:length(yu) ncl=yu(i); ncli=find(y==ncl); simca_model.models(i).x=x(ncli,:); simca_model.models(i).clas=y(ncli,:); end %---CV--- for inc=1:length(yu) wb=waitbar(0,'Please wait...','Name','Performing cross-validation'); x0=simca_model.models(inc).x; y0=simca_model.models(inc).clas; if min(size(x0))<=20; ncp=(min(size(x0)))-1; else ncp=20; end for incp=1:ncp waitbar(incp/ncp); for iam=1:size(x0,1) xv=x0(iam,:); xc=x0; xc(iam,:)=[]; yv=y(iam,:); yc=y; yc(iam,:)=[]; [E L varpc eive]=pcar(xc); ns=xv*L(:,1:incp); r=sum((ns*(L(:,1:incp))'-xv).^2); rCP(iam,incp)=r; vCP(incp)=sum(varpc(1:incp)); end end close(wb); press=sum(rCP); figure; subplot(1,2,1); plot(press,'-o'); t = sprintf('PRESS for class %g',yu(inc)); title(t) xlabel('No. of PC on model'); ylabel('PRESS'); subplot(1,2,2); plot(vCP,'-o'); t = sprintf('Explained variance for class %g',yu(inc)); title(t) xlabel('No. of PC on model'); ylabel('Explained variance (%)'); cpret = input('How many PC do you want to retain? '); close; clear rCP; clear vCP; %---PCA--- [E L varpc eive]=pcar(x0); L=L(:,1:cpret); E=x0*L; %---T--- if cpret>1 t = E*inv(diag(eive(1:cpret))); T = sum((t.^2)')'; else T = E.^2/eive(1,1); end %---Q--- res = (x0-E*L')'; Q = (sum(res.^2))'; %---Tlim--- [m,n]=size(x0); if m > 300 Tlim = (cpret*(m-1)/(m-cpret))*finv(0.95,cpret,300); elseif m == cpret; Tlim = 0; else Tlim = (cpret*(m-1)/(m-cpret))*finv(0.95,cpret,m-cpret); end %---Qlim--- me=size(eive,1); if cpret < min([m n]) cl=2*95-100; t1 = sum(eive(cpret+1:me,1)); t2 = sum(eive(cpret+1:me,1).^2); t3 = sum(eive(cpret+1:me,1).^3); h0 = 1-2*t1*t3/3/(t2^2); if h0<0.001 h0=0.001; end ca = sqrt(2)*erfinv(cl/100); h1 = ca*sqrt(2*t2*h0^2)/t1; h2 = t2*h0*(h0-1)/(t1^2); Qlim = t1*(1+h1+h2)^(1/h0); else Qlim = 0; end simca_model.models(inc).S=E; simca_model.models(inc).L=L; simca_model.models(inc).ev=eive(1:cpret); simca_model.models(inc).T2=T; simca_model.models(inc).Q2=Q; simca_model.models(inc).T2lim=Tlim; simca_model.models(inc).Q2lim=Qlim; end simca_model.stats.classes=yu; %---pred calib--- pred=simcapred(x,simca_model,classes); simca_model.stats.T2=pred.T2; simca_model.stats.Q2=pred.Q2; simca_model.stats.class_c=pred.class_p; simca_model.stats.RMSEC=pred.RMSEP; simca_model.stats.R2=pred.R2; simca_model.stats.suc=pred.suc; %--- PCA function --- function [E L varpc eive]=pcar(x) m=size(x,1); [u s L]=svd(x,'econ'); % [U,S,V] = SVD(X,'econ') also produces the "economy size" % decomposition. If X is m-by-n with m >= n, then it is % equivalent to SVD(X,0). For m < n, only the first m columns % of V are computed and S is m-by-m. E=u*s; eive=diag(s); %取对角元素 var=diag(s).^2; varpc=var/sum(var)*100;