www.gusucode.com > mbctools 工具箱 matlab 源码程序 > mbctools/@mdev_local/validate.m
function [BMIndex,Table]= validate(mdev,hFig) %VALIDATE % % [BMIndex,Table]= validate(mdev,hFig) % Copyright 2000-2015 The MathWorks, Inc. and Ford Global Technologies, Inc. TSModels = mdev.TwoStage; selrf = mdev.ResponseFeatures; [Xag,Y]= getdata(mdev,'FIT'); X = Xag{1}; % Local Input Data XG = Xag{end}; % Global Input Data Y(:,:,~mdev.FitOK)= NaN; L= model(mdev); G= children(mdev,'BestModel'); % Get Response Feature data from Local Fit Yf= mdev.RFData.info; Yf= double(Yf); Nf= length(get(L,'Values')); if size(Yf,2)>Nf rf1=RFstart(L); Yf(:,rf1)=[]; end % You can only calculate PRESS for linearmods CalcPRESS=1; for i=1:length(G) CalcPRESS= CalcPRESS & islinear(G{i}); end switch DatumType(L) case {1,2} dm= children(mdev,1,'model'); CalcPRESS = CalcPRESS & islinear(dm{1}); case 3 pdatum= datumlink(mdev); CalcPRESS = CalcPRESS & islinear(pdatum.model); end N= size(X,3); sigma2= mdev.MLE.Pooled_MSE; LogLi= NaN(N,length(TSModels)); Yf_pred= zeros(size(Y,3),length(TSModels)); RMSE= Yf_pred; XGts= [zeros(size(XG,1),size(X,2)) double(XG)]; Lfpred = cell(1,length(TSModels)); SSr1= Y; Yr= double(Y); if CalcPRESS Table= zeros(length(TSModels),4+CalcPRESS); Yf_PRESS= Yf_pred; RPRESS= Yf_pred; SSr2= Y; else Table= zeros(length(TSModels),4); end if nargin>1 && isgraphics(hFig) mbH = MBrowser; DO_STATUS_UPDATE = true; setProgressBar(mbH, 'value', 0); msgID = mbH.addStatusMsg('Building Two-stage models...'); else DO_STATUS_UPDATE = false; end Y_pred = cell(1, length(TSModels)); Y_PRESS = cell(1, length(TSModels)); for i= 1:length(TSModels) % log likelihood TS= TSModels{i}; doMLE = canMLE(TS); if doMLE [Xg,Yrf,Sigma] = mledata(mdev,i); % initialise pev for univariate case TS= pevinit(TS,Xg,Yrf,Sigma,1); TS = mle_Algorithm(TS); [f,G{i},Li,OK]= loglikFcn(TS,Xg,Yrf,Sigma); LogLi(OK,i)=Li; Table(i,end)= f; end % Ordinary Predictions Y_pred{i}= TS(Xag); SSr1(:,1)= Yr-Y_pred{i}; [~,Lfpred{i}]= EvalModel(TS,XGts); if CalcPRESS % PRESS Predictions [Y_PRESS{i},LfPRESS] = presspred(TS,{X,double(XG)}); SSr2(:,1)= Yr-Y_PRESS{i}; end rf= selrf(i,:); for k=1:N % for i=1:length(TSModels) % Response Feature residuals % need to account for ytrans here !!!!!!!!!!!???????????? if mdev.FitOK(k) if doMLE % Response Feature Covariance Matrix for sweep k S= squeeze( mdev.MLE.Sigma(rf,rf,k) )*sigma2 + G{i}; else S=NaN; end [Yf_pred(k,i),RMSE(k,i)]= i_CalcRFStats(Yf(k,rf),Lfpred{i}(k,:),S,SSr1(:,:,k)); if CalcPRESS [Yf_PRESS(k,i),RPRESS(k,i)]= i_CalcRFStats(Yf(k,rf),LfPRESS(k,:),S,SSr2(:,:,k)); end else Yf_pred(k,i)= NaN; RMSE(k,i) = NaN; Yf_PRESS(k,i)= NaN; RPRESS(k,i)= NaN; end end if DO_STATUS_UPDATE setProgressBar(mbH,'value',i./length(TSModels)); end end if DO_STATUS_UPDATE setProgressBar(mbH,'value', 0); mbH.removeStatusMsg(msgID); end sn= statistics(mdev); for i= 1:length(TSModels) ResAll= double(Y)- Y_pred{i}; rOK= isfinite(ResAll); if any(rOK) ResAll=ResAll(rOK); VarR= sum(ResAll.^2,1)/(length(ResAll)); else VarR=0; end s= sqrt(VarR); % Look at Response Features resids= Yf_pred(:,i); residsOK= all(isfinite(resids),2)'; resids= resids(residsOK,:); dfFeat= length(resids); if dfFeat VarRF(1)= ( sum(resids,1)./dfFeat ); else VarRF(1)= 0; end if CalcPRESS % PRESS RMSE ResAll= double(Y)- Y_PRESS{i}; rOK= isfinite(ResAll); if any(rOK) ResAll=ResAll(rOK); VarPRESS= sum(ResAll.^2,1)/(length(ResAll)); else VarPRESS= 0; end sp= sqrt(VarPRESS); Table(i,1:end-1)= [sn(1) s sp VarRF(1)]; else Table(i,1:end-1)= [sn(1) s VarRF(1)]; end end mdev.TSstatistics.Summary= Table; mdev.TSstatistics.RespFeat= Yf_pred; mdev.TSstatistics.LogL = LogLi; mdev.TSstatistics.RMSE = RMSE; if CalcPRESS mdev.TSstatistics.RPRESS = RPRESS; [dum,ind]= min(Table(:,3)); else [dum,ind]= min(Table(:,2)); end if ~isempty(Table) BMIndex= ind(1); mdev.MLE.Validate=0; if isfield(mdev.MLE,'BestModel') mdev.MLE= rmfield(mdev.MLE,'BestModel'); end if isfield(mdev.MLE,'Model') mdev.MLE= rmfield(mdev.MLE,'Model'); end pointer(mdev); else BMIndex= []; end return % sweep residual stats function [Yf_pred,RMSE]= i_CalcRFStats(Yf,Lf,S,SSr) r= Yf - Lf; % T^2 stats if all(isfinite(r)) && all(isfinite(S(:))) && rank(S) == size(S,1) Yf_pred = r/S*r'; else Yf_pred = NaN; end % Predicted RMSE res= double(SSr); rok= isfinite(res); if any(rok) RMSE = sqrt(sum(res(rok).^2)/sum(rok)); else RMSE = NaN; end