www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregtwostage/diagnosticStats.m

    function Diags= diagnosticStats(TS,Xg,Yrf,Sigma)
%DIAGNOSTICSTATS diagnostic stats for Global two-stage
%
% Diags= diagnosticStats(TS,Xg,Yrf,Sigma)
%  Diags is a structure with fields
%    'Observed','Yhat''Residuals','SResiduals'

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

[Xgc,Yrf,W0]=checkdata(TS,Xg,Yrf,Sigma);

[TS,J,y,W0,S,SF]= mle_scale(TS,Xgc,Yrf,W0);

Nf= size(Yrf,2);
Ns= size(Yrf,1);


Wci= choltinv(TS.covmodel,W0);
Wf= cov(TS.covmodel,W0);
D= unstruct(TS.covmodel);

% degrees of freedom for MLE model
% df= numel(Yrf)-numParams(TS)-length(TS.covmodel) - length(covmodel(TS.Local))-1; 
OK= true;

% Predicted RF Values
P= zeros(size(Yrf));
for i= 1:Nf
   P(:,i)= eval(TS.Global{i},Xgc);
end

% Error
e= Yrf-P;	%Equation 13



% Calculate stats
% Correlation of G

Zc= Wci*J;
% make sure we use full qr algorithm - sparse is nasty

ri= var(TS);
if isempty(ri)
	R=qr(full(Zc));
	R= triu(R(1:size(R,2),:));
	% Predicted Variances
	rd= abs(diag(R));
	tol= length(Zc)*eps*max(rd);
    if tol==0
        tol=eps;
    end
    OK= min(abs(rd))<tol;
    if OK
        ri= inv(R);
    end

end

if OK
 
    Jw= J*ri;
    VarP = sum((Jw).^2,2);	%Equation 18
    % Error Variances
    H=  sqrt(diag(Wf)-VarP);

    % Parameter Variance
    % VarB = ri*ri';						%Equation 16

    % Reshape Data into Response Feature Order
    Zind= reshape(1:Nf*Ns,Nf,Ns)';

    Wc= diag(sqrt(diag(D)));
    H= reshape(H(Zind),Ns,Nf)*Wc;

    VarE= e./H*SF;
else
    % Make studentised residuals NaN
	VarE= NaN(size(e));
end
% Output Structure
Diags= struct('Observed',Yrf,'Yhat',P,'Residuals',e,'SResiduals',VarE);