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

    function mdev= Change_RespFeat(mdev,SweepNos)
%CHANGE_RESPFEAT response feature data update
%
%  mdev= Change_RespFeat(mdev);                % change all sweeps/rfs

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


if nargin<2
    SweepNos = 1:size(mdev.AllModels,2);
end

L= model(mdev);

% this is now a way to get clean data
[X,Yresp]= getdata(mdev);
[X,Y,DataOK,bd]= checkdata(L,X,Yresp);
Yresp(bd)= NaN;
% bad sweeps
Yresp(:,:,~DataOK)= NaN;
if any(DataOK)
	X(:,:)= invcode(L,double(X));
end
if ~isempty(mdev.GLSWeights)
	Wc= mdev.GLSWeights;
else
	Wc= cell(1,size(X,3));   
end

if isnull(mdev.RFData)
    % make response feature pointer if required
    XF =getdata(mdev,'fit');
    XG = XF{2};
    S= sweepset('Variable', 1 , '%g' , {}, 'Response Feature'  , {} , 'none' , zeros(size(XG,1),0) );
    D= [XG(:,~1) S];
    pD=xregpointer(D);
    mdev.RFData = pD;
end


Nf= numfeats(L);
Ns= length(SweepNos);
NewRF=zeros(Ns,Nf);
nl= nfactors(L);
D= zeros(Ns,nl);


if isempty(mdev.MLE) || size(mdev.AllModels,2) ~= size(mdev.MLE.Sigma,3)
	% size of data may have changed
	Ns= size(mdev.AllModels,2);
	Sigma   = zeros(Nf,Nf,Ns);
	SSE_nat = zeros(Ns,1);
	sse     = zeros(Ns,1);
	df      = ones(Ns,1);
    ValRMSE    = NaN(Ns,1);
    ValDataSize= zeros(Ns,1);
else
	Sigma   = mdev.MLE.Sigma;
	SSE_nat = mdev.MLE.SSE_nat;
	sse     = mdev.MLE.SSE;
	df      = mdev.MLE.df;
    if ~isfield(mdev.MLE,'ValRMSE')
        % initialise validation data storage
        mdev.MLE.ValRMSE = NaN(1,size(mdev.AllModels,2));
        mdev.MLE.ValDataSize = NaN(1,size(mdev.AllModels,2));
    end
    ValRMSE    = mdev.MLE.ValRMSE;
    ValDataSize= mdev.MLE.ValDataSize;
end


if Nf~=size(Sigma,1)
	Sigma   = zeros(Nf,Nf,Ns);
end

mdev.IsLinearised = false;
m=1;

[Xval,Yval]= getLocalValidationData(mdev,SweepNos);
for testCounter= 1:length(SweepNos)
    testNumber = SweepNos(testCounter);
	if mdev.FitOK(testNumber)
		j = sum(DataOK(1:testNumber));
		L= LocalModel(mdev,testNumber);
		Xc= code(L,X{j});
		D(m,:)= datum(L);
		NewRF(m,:)= evalfeatures(L);
		% need to re-evaluate covariance matrix as well
		if ~isempty(Wc{testNumber}) && length(Y{j})~=size(Wc{testNumber},1)
            Wc{testNumber} = CalcWeights(L,Xc);
			mdev.GLSWeights{testNumber}= Wc{testNumber};
		end
		
		[SummStats,sig]= localstats(L,Xc,Y{j},Wc{testNumber});
		
		SSE_nat(testNumber) = SummStats.SSE_natural;
		sse(testNumber)     = SummStats.sse;
		df(testNumber)      = SummStats.df;
        
        if ~isempty(Xval) && ~isempty(Xval{testCounter})
            % only calculate validate RMSE if available
            [ValRMSE(testNumber),ValDataSize(testNumber)] = LocalValidationRMSE(mdev,testNumber,L,Xval(testCounter),Yval(testCounter));
        end
        
        if ~isempty(sig)
			Sigma(:,:,testNumber)= sig;
        end
	else
		SSE_nat(testNumber) = Inf;
		sse(testNumber)     = Inf;
		df(testNumber)      = 0;
		
		NewRF(m,:)  = NaN;
		D(m,:)        = NaN;
	end
	m=m+1;
end

mdev.MLE.Sigma  = Sigma  ;
mdev.MLE.SSE_nat= SSE_nat;
mdev.MLE.SSE    = sse;
mdev.MLE.df     = df; 
mdev.MLE.ValRMSE = ValRMSE;
mdev.MLE.ValDataSize = ValDataSize;

if RFstart(L)
	NewRF= [D NewRF];
end

% set up guids for response feature data and remove any global outliers which change
sguids = getSweepGuids(Yresp,'goodonly');


RF = mdev.RFData.info;
RF  = setGuids(RF,sguids);
% update rf data
if ~isempty(NewRF)
    RF(SweepNos,:)= NewRF;
end
mdev.RFData.info= RF;
mdev   = updateOutlierIndices(mdev,Yresp,RF);


% Update statistics
NotFitted = ~mdev.FitOK;
mdev.MLE.SSE(NotFitted)= 0;
mdev.MLE.SSE_nat(NotFitted)= 0;
mdev.MLE.df(NotFitted)= 0;
mdev.MLE.ValRMSE(NotFitted) = NaN;
mdev.MLE.ValDataSize(NotFitted) = 0;

SSE= mdev.MLE.SSE;
df = mdev.MLE.df;

% Summary Statistics for local model
DFT= sum(df(mdev.FitOK));
if DFT>0
    Pooled_MSE= sum(SSE(mdev.FitOK))/DFT;
else
    Pooled_MSE= NaN;
end
mdev.MLE.Pooled_MSE= Pooled_MSE;

SSE= mdev.MLE.SSE_nat;
ind= isfinite(SSE(:)) & mdev.FitOK(:);
DFT= sum(df(ind));
if DFT>0
    Pooled_MSE= sum(SSE(ind))/DFT;
else
    Pooled_MSE= NaN;
end

% mdev_local summary stats
ch  = colhead(mdev);
stats = [sqrt(Pooled_MSE) NaN(1,length(ch)-1)];
if isPointByPoint(mdev)
    %  point-by-point validation statistics
    sse = mdev.MLE.ValRMSE.^2.*mdev.MLE.ValDataSize;
    HasVal = isfinite(sse);
    if any(HasVal)
        ValRMSE = sqrt(sum(sse(HasVal))/sum(mdev.MLE.ValDataSize(HasVal)));
        stats(2) = ValRMSE;
    end
end
mdev= statistics(mdev,stats);


pointer(mdev);