www.gusucode.com > mbctools 工具箱 matlab 源码程序 > mbctools/@mdev_local/mle_validate.m

    function [Table,Models,mdev]= mle_validate(mdev,hFig)
%MLE_VALIDATE

%  Copyright 2000-2004 The MathWorks, Inc. and Ford Global Technologies, Inc.



if ~isfield(mdev.MLE,'Validate') | ~mdev.MLE.Validate
   ModelNo= BMIndex(mdev);
else
   ModelNo= 1;
end   


TSmle= mdev.MLE.Model;
[Xg,Yrf,Sigma] = mledata(mdev,1,mdev.MLE.Modes(2));
[f,Gamma,Li,OK]= loglikFcn(TSmle,Xg,Yrf,Sigma);
sigma2= mdev.MLE.Pooled_MSE;
Ng= size(Gamma,1);

bd= ~OK;

rf = mdev.ResponseFeatures(1,:);

% get global data

presp= Parent(mdev);
L= model(mdev);
G= children(mdev,'BestModel');

[Xts,Y]= getdata(mdev,'FIT');
X = Xts{1};
XG= Xts{2};
Y(:,:,~mdev.FitOK)= NaN;



Yf_pred= zeros(size(Y,3),1);


% Get Response Feature data from Local Fit
Yf = mdev.RFData.double;
if RFstart(L)
   Yf(:,1)=[];
end

nl= nfactors(L);

RMSE= zeros(size(XG,1),1);
SSr1= Y;

[Y_pred,Lfpred]= EvalModel(TSmle,Xts);

SSr1(:,1)= double(Y)-Y_pred;

N= size(Y,3);
for k=1:N											
   
   if mdev.FitOK(k)
      % Response Feature residuals 
      % need to account for ytrans here !!!!!!!!!!!????????????
      r= Lfpred(k,:)-Yf(k,rf);
      % Response Feature Covariance Matrix for sweep k
      S= squeeze( mdev.MLE.Sigma(rf,rf,k) )*sigma2 + Gamma;
      
      if all(isfinite(r)) & all(isfinite(S(:))) & rank(S)== size(S,1)
         Yf_pred(k,1) = r/S*r';
      else
         Yf_pred(k,1) = NaN;
      end
      
      r= double(SSr1(:,:,k));
      rok= isfinite(r);
      if any(rok)
         RMSE(k) = sqrt(sum(r(rok).^2)/sum(rok));
      end
   else
      RMSE(k)= NaN;
      Yf_pred(k,1)= NaN;
   end
   
end

ResAll= double(Y)- Y_pred; %(i)
rOK= isfinite(ResAll);
ResAll=ResAll(rOK);
if ~isempty(ResAll)
	VarR= sum(ResAll.^2,1)/length(ResAll);
	s= sqrt(VarR);
else
	s= NaN;
end

% Look at Response Features
Yf_pred(bd)=NaN;
resids= Yf_pred;
residsOK= all(isfinite(resids),2);
resids= resids(residsOK,:);
dfFeat= length(resids);
VarRF(1)= ( sum(resids,1)./dfFeat );

% Sweep Contributions to log L
LogLi= zeros(N,1);
LogLi(~bd,1)=Li;
LogLi(bd,1)= NaN;

if isfield(mdev.TSstatistics,'Summary') & ~isempty(mdev.TSstatistics.Summary)
	s_un= mdev.TSstatistics.Summary(1,:);
	switch length(s_un)
	case 4
		if islinear(mdev.TwoStage{1})
			% need press r,mse
			s_un= [s_un(1:2) NaN s_un(3:end)];
		end
	case 5
	otherwise
		% old summary table of wrong size
		s_un= statistics(mdev);
	end
else
	s_un= statistics(mdev);
end

sLp= yinv( model(mdev) ,sqrt(sigma2));

% display [Local RMSE , Pred RMSE , Pred T^2 , -log L]
% log L is last element in table
chead= colhead(mdev);
lind= length(chead);
if length(s_un)==5
	% use NaN for MLE PRESS RMSE
	Table= [s_un;s_un(1) s NaN VarRF(1) f];
else
	Table= [s_un;s_un(1) s VarRF(1) f];
end

% store the results
mdev.ResponseFeatures= mdev.ResponseFeatures([ModelNo,ModelNo],:);
mdev.TSstatistics.Summary= Table;
mdev.TSstatistics.LogL = [mdev.TSstatistics.LogL(:,ModelNo) LogLi];
mdev.TSstatistics.RespFeat= [mdev.TSstatistics.RespFeat(:,ModelNo) Yf_pred];
mdev.TSstatistics.RMSE= [mdev.TSstatistics.RMSE(:,ModelNo) RMSE];

% twostage info setup
%[mdev.MLE.Model        ,mdev] = tsinfo(mdev,mdev.MLE.Model);
%[mdev.TwoStage{1},mdev] = tsinfo(mdev,mdev.TwoStage{1});

% initialise PEV calcs
[Xg,Yrf,Sigma] = mledata(mdev,1);
if ~pevcheck(mdev.MLE.Model)
   mdev.MLE.Model= pevinit(mdev.MLE.Model,Xg,Yrf,Sigma);
end
if ~pevcheck(mdev.TwoStage{1})
   mdev.TwoStage{1}= pevinit(mdev.TwoStage{1},Xg,Yrf,Sigma);
end
pointer(mdev);

Models= {mdev.TwoStage{1} mdev.MLE.Model};

mdev.MLE.Validate=1;

mdev.TwoStage= Models;

pointer(mdev);