% model discrimination
% 多项式模型 Model: y=beta1+beta2*t+beta3*t^2+...+betap*t^(p-1)
% To test whether certain subsets of parameters are needed in the model
% The test is based on the relative reduction in residual sum of squares.
% 假设：beta_new=0，如引入该参数后预测残差平方和下降较大，则假设无效
clear all; format short g
n=17; %Measured temperature in Celsius.  Time step is 96 sec.
y=[137.55 129.37 121.96 115.17 108.99 102.91 97.7 92.98 88.69 84.69...
   80.91 77.45 73.91 71.05 68.44 65.99 63.71]';
R=y'*y; % 初始残差平方和，起始时尚未引入参数，认为所有参数为0，残差即为y值
t=0:1:n-1;t=t';%Only for compactness is more than one statment on line
for m=1:6 % 待分析参数的个数设为 p=5
   X=ones(n,1); % 针对截距参数 beta1
   for jj=1:m-1
       tj=t.^jj; X=[X tj];
   end
   b=Xy; % OLS estimate; same as b=(X'*X)X'*y;
   e=y-X*b; % 计算残差
   sum=e'*e; % 残差平方和，自由度为n-p,p=m
   DR=R-sum; % 残差平方和降幅(服从卡方分布），因为只有一个值，自由度为1
   s2=sum/(n-m); % 自由度：n-m
   Con=cond(X); % condition number
   mm=n-m;
   F=DR/s2; % F(q,n-p)分布，q=1为因变量个数，p为参数个数
   % Fcr=finv(0.95,n-p),只有一尾所以为0.95；对于t检验，Tcr=tinv(0.975,n-p)
   K=m+1;
   AICc=n*log(sum/n)+2*K+2*K*(K+1)/(n-K-1); % Akaike criterion
   % The model with lower AICc will be more likely to be correct
   B(m,:)=[m mm sum s2 DR F Con AICc]; % 随着增加参数，逐行记录统计结果
   % 参数个数，自由度，残差平方和,自由度平均残差平方和，残差降幅，F值，Cond num. AICc
   R=sum; % 逐个参数评价，残差平方和依次更新
   ypred=X*b; % 预测值
end
disp('       p    deg. of fr.     R       Mean sq.      Delta R        F      Cond(X)   AICc');
for i=1:6 % fprintf不能输出矩阵，以loop实现
fprintf('%8i %8i %14.2f %10.2f %12.2f %10.2f %10.1f %10.2fn',B(i,:))
% %d or %i for integer format; %f for decimal format
end
% Fcrit=finv(.95,1,n-3)
% disp(B); aa    

% Comment:

% F(.95,1,17-p) is 4.75 and 4.84 for p=5 and p=6, respectively

% Then at the 5% level of significance, the p=5 model is selected.    

       p    deg. of fr.     R       Mean sq.      Delta R         F           Cond(X)       AICc

      1       16        8526.11     532.88    147120.43     276.08             1.0       110.56
      2       15         276.07       18.40        8250.03     448.25            18.1        55.23
      3       14            4.98        0 .36          271.09     761.65          331.3        -9.53
      4       13            0.34         0.03             4.65      178.71        6032.2      -51.16
      5       12            0.22         0.02             0.12         6.20     114976.0      -53.29
      6       11            0.18         0.02               0.4         2.59   2388793.7      -50.84