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

    function [data,factors,specialPlots,olIndex]= diagnosticStats(m,X,Y,NaturalResiduals)
%DIAGNOSTICSTATS diagnostic statistics for linear models
%
%  [data,factors, specialPlots]= DIAGNOSTICSTATS(m)
%
%  This is an overloaded function to return stats and factors for the
%  diagnostic plots (scatter plots in browser)

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

% This gets the standard statistics
if nargin<4
    NaturalResiduals = true;
end

[Xm,y] = checkdata(m,X,Y);

if isempty(m.Store)
   m = InitModel(m,Xm,y); 
end
yhat = eval(m,Xm);

H = m.Store.H;
[~,mse,df]= var(m);
nobs = length(y);

% leverage values
hneg= 1-H;
% make sure this is a small positive number
hneg(hneg<=0)=eps;

residuals = y - yhat;

if isfinite(mse)
    sse = sum(residuals.^2);
    % Studentised residuals
    if all(m.lambda==0) && df>1
        % externally scaled residuals
        SSExt= (sse*hneg- residuals.^2)/(df-1);
        SSExt = max(SSExt,eps);
        studext=residuals./sqrt(SSExt);
    else
        studext = zeros(size(y));
        studext(:)= NaN;
    end
    % internally scaled residuals
    studint= residuals./sqrt(sse/df*hneg);
else
    studext = zeros(size(y));
    studext(:)= NaN;
    studint= studext;
end
p= nobs-df;
if p>0
    % Cooks Distance
    cookd = studint.^2 .*(H./(hneg))./p;
else
    cookd = zeros(size(residuals));
    cookd(:)= NaN;
end

% Plot indexes.
obs = (1:nobs)';
% Get the global variables.
Xm= invcode(m,Xm);

if NaturalResiduals && (isTransOneSide(m) || isTBS(m))
    % use natural units
    y = yinv(m,y);
    yhat = yinv(m,yhat);
    if ~isreal(yhat)
        yhat(abs(imag(y))>eps)= NaN;
        yhat =real(yhat);
    end
    
    residuals = y - yhat;
end


data= [yhat,studext,studint,residuals,cookd,H,obs,y,Xm];

% This returns the names of the statistics
sname = InputLabels(m);
% make it a row
sname= sname(:)';

output= getOutput(m);
if NaturalResiduals && ~isempty(output.Units)
    units = sprintf(' [%s]',output.Units);
    snameY= ResponseLabel(m);
else
    units = '';
    snameY = varname(m);
end
factors={['Predicted ',snameY],...
    'Studentized residuals (External)',...
    'Studentized residuals (Internal)',...
    ['Residuals',units],...
    'Cook''s distance',...
    'Leverage',...
    'Obs. number',...
    snameY};

if any(m.lambda~=0)
    pos= strcmp(factors,'Studentized residuals (External)');
    % delete external studentised residuals when there is some ridge
    % regresion
    factors(pos)=[];
    data(:,pos)= [];
end

factors= [factors,sname];

% This is a list of special plots. This can be empty
specialPlots= {'Predicted/Observed','Normal plot'};
% outlier index.
olIndex= outliers(m,data,factors);