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

    function [s,list,SummStatsType,Description] = summary(m,Criteria,X,y,s2full,s2mle)
%SUMMARY Calculate summary statistics
%
%  [s,list,SummStats] = SUMMARY(m, Criteria, X, Y, s2full)
%  [MINorMAX,list,SummStatsType] = SUMMARY(m) returns a list of criteria
%  and whether to use as summary:
%    -1 special calculation done differently in summary stats
%     1 use as calculated
%     2 only compare if data and transforms the same for all models
%     0 don't use for comparison in summary table (only in prune)
%
%  Criteria list = {'PRESS RMSE','RMSE','GCV','Weighted  PRESS',...
%                  '-2logL','AIC','AICc','BIC','R^2','R^2 adj','Cp'}

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


list= {'PRESS RMSE','RMSE','GCV','Weighted PRESS','-2logL','AIC','AICc','BIC','R^2','R^2 adj','PRESS R^2','DW','Cp'};
SummStatsType= [-1 -1 1 1 2 2 2 2 1 1 1 1 0];
if nargout>3
    Description= {'Predicted Standard Error'        'sqrt(PRESS/N)'
        'Root Mean Square Error'                    'sqrt(SSE/(N-p))'
        'Generalized Cross-validation Variance'     'N*SSE/(N-p)^2'
        'Weighted Predicted Standard Error'         'sqrt(PRESS/(N-p-1))'
        '-2 * log likelihood'                       'N*log(SSE/N)'
        'Akaike Information Criteria'               '-2logL + 2*(p+1)'
        'Small Sample Akaike Information Criteria'  '-2logL + 2(p+1)*N/(N-p)'
        'Bayesian Information Criteria'             '-2logL + 2*log(N)*(p+1)'
        'R^2'                                       '1 - SSE/SST'
        'Adjusted R^2'                              '1 - SSE/SST*(N-1)/(N-p)'
        'PRESS R^2'                                 '1 - PRESS/SST'
        'Durbin-Watson Statistic'                   'sum((e_i-e_{i+1})^2)/sum((e_i^2) '
        'Mallow''s Statistic'                       'SSE/(SSEmax/(N-pmax)) - N + 2*p'};
end


if nargin<2
    % true indicates minimise and false maximise
    s= true(size(list));
    s(9:11)= false;
    return
end

if ischar(Criteria)
    switch lower(Criteria)
        case 'all'
            sel= 1:length(list);
        case 'data'
            sel= find(abs(SummStatsType)~=1);
        case 'summary'
            sel= find(abs(SummStatsType)~=0);
        case 'multi'
            sel= find(SummStatsType<=0);
        otherwise
            sel= find(ismember(list,Criteria));
    end
elseif iscell(Criteria)
    sel= find(ismember(list,Criteria));
elseif islogical(Criteria)
    sel= find(Criteria);
else
    sel= Criteria;
end

nObs=size(y,1);

[ri,mse,df] = var(m);% degrees of freedom

[anova,R2,R2adj] = AnovaTable(m.mv3xspline);

sst= anova(end,1);
sse= anova(2,1);

if nargin<6
    s2mle= sse/nObs;
end
if nargin<5
    s2full = mse;
end

nk= get(m.mv3xspline,'numknots');
K= numParams(m.mv3xspline) + nk + 1;

s= zeros(1,length(sel));
for i=1:length(sel)
    s(i) = NaN;
    switch sel(i)
        case 1
            % PRESS RMSE
            mspline = set(m.mv3xspline,'ytrans',get(m,'ytrans'));
            s(i) = NaturalPressRMSE(mspline,X,y);
        case 2
            % RMSE
            s(i) = NaturalRMSE(m,X,y);
        case 3
            % GCV
            s(i) = calcGCV(m.mv3xspline);
        case 4
            % Weighted PRESS
            if df>1
                s(i) = sqrt(press(m.mv3xspline,y)/(df-1));
            end
        case 5
            % log likelihood  actually -2logL
            if s2mle>0
                s(i)= nObs*log(s2mle);
            end
        case 6
            % Akaike Information Criteria
            if s2mle>0
                s(i) = nObs*log(s2mle) + 2*K;
            end
        case 7
            % AICc Akaike Information Criteria for small samples
            if s2mle>0 && nObs>K+1
                s(i) = nObs*log(s2mle) + 2*K*nObs/(nObs-K-1);
            end
        case 8
            % Bayesian Information Criteria
            if s2mle>0
                s(i) = nObs*log(s2mle) + log(nObs)*K;
            end
        case 9
            % R^2
            s(i)= R2;
        case 10
            %Adjusted R-squared statistic
            s(i)= R2adj;
        case 11
            if sst>0
                s(i)= 1-press(m.mv3xspline,y)/sst;
            else
                s(i)= 1;
            end
        case 12
            % Durban-Watson
            yhat= eval(m,X);
            r= y - yhat;
            s(i)= sum( (r(1:end-1)- r(2:end)).^2 ) /sum(r.^2);
        case 13
            % Mallow's Statistic
            s(i)= sse/s2full - df + numParams(m);
    end
end