% Example 2.9
% using F test for polynomial model building
%
clear all
sigma = 0.02;  % noise level
beta = [ 1 0.1 -0.005 ]';   % true values for the parameters
t = [ 0:0.4:10]';
n = length(t);
X = [ ones(n,1) t t.^2 ];  % for the true model
% generate random noise for measurements
rng(129,'v5normal');  % use same seed value
eps = sigma*randn( n,1 );
y = X*beta + eps;  % the simulated measurements
figure(1); clf;
plot(t, y, 'linewidth', 1.0 )
hold on; 
axis( [ 0 10 -0.2 1.6 ] )
% now test the models polynomial 1 through 5
clear X
R_last = y'*y;  
fprintf('  i   dof   R_LS      s-sqr.      delta-R        F      F_test     cond(X)\n');
for p = 1:5
    clear X
    dof = n - p;
    X = zeros( n, p );  % just set it up
    X(:,1) = ones(n,1); % always ones in first column
    for ip = 2:p
        X(:,ip) = t.^(ip-1);
    end
    b_LS = (X'*X) \ X'*y;
    e = y - X*b_LS;
    R_LS = e'*e;          % Eq. (2.4.13)
    s_squared = R_LS / (n-p);
    delta_R = R_last - R_LS;
    F = delta_R/s_squared;
    F_test = finv( 0.95, 1, n-p );
    if( F/F_test > 1.0 )    % save largest value of p for F/F_test > 1
        b_opt = b_LS;
        p_opt = p;
        e_opt = e;
    end
    fprintf(' %2d   %2d  %8.5f   %8.5f   %9.5f   %9.3f   %6.2f   %8.2f\n', ...
              p, dof, R_LS, s_squared, delta_R, F, F_test, cond(X) );
    R_last = R_LS;
end
plot( t, 4*e_opt, 'o', 'MarkerEdgeColor', [ 0.1 0.5 0.1], 'linewidth', 1.0)
plot( [ 0 10], [0 0], 'r-', 'linewidth', 1.0 )
xlabel('t')
ylabel('e                                           y                   ')
text(5,0.2,'4 times residuals')
text(5,1.2,'simulated measurements')
% print the 600 dpi figure to tiff format
fig=figure(1);
print(fig,fullfile(pwd, 'Figure2_3.tif'),'-dtiff','-r600')


