% Listing 2.2 of the book
% For Example 2.4, Page 34

% Arun K. Tangirala
% March 21, 2015
% Freely distributed under no warranty

%% Generation of training data
xvec = (0:0.02:4)';
yvec = 1.2 + 0.4*xvec + 0.3*xvec.^2 + 0.2*xvec.^3;

% Add measurement noise
ymeas = yvec + 0.12*randn(length(yvec),1);

%% Polynomial fits

% 1st order fit
[pvechat1,S1] = polyfit(xvec,ymeas,1);
% 2nd order fit
[pvechat2,S2] = polyfit(xvec,ymeas,2);
% 3rd order fit
[pvechat3,S3] = polyfit(xvec,ymeas,3);
% 4th order fit
[pvechat4,S4] = polyfit(xvec,ymeas,4);
% 5th order fit
[pvechat5,S5] = polyfit(xvec,ymeas,5);
% 6th order fit
[pvechat6,S6] = polyfit(xvec,ymeas,6);
% 7th order fit
[pvechat7,S7] = polyfit(xvec,ymeas,7);

% Calculation of standard error in estimates
Rinv = inv(S3.R);
covmat3 = (Rinv*Rinv')*S3.normr^2/S3.df;
errvec3 = sqrt(diag(covmat3))';

Rinv4 = inv(S4.R);
covmat4 = (Rinv4*Rinv4')*S4.normr^2/S4.df;
errvec4 = sqrt(diag(covmat4))';

% Test data for model validation
xtest = (4:0.02:8)';
ytest = 1.2 + 0.4*xtest + 0.3*xtest.^2 + 0.2*xtest.^3;
ytmeas = ytest + 0.12*randn(length(ytest),1);

% Predictions of the model on test data
[yhat3,delta3] = polyval(pvechat3,xtest,S3);
[yhat4,delta4] = polyval(pvechat4,xtest,S4);
[yhat5,delta5] = polyval(pvechat5,xtest,S5);

% Plot training data
plot(xvec,ymeas)
xlabel('Input'); ylabel('Output')
%prepare_figure

% Plot the residual norm vs. order
errnorm = [S1.normr S2.normr S3.normr S4.normr S5.normr S6.normr S7.normr]';
figure
plot((1:7),errnorm,'-o','Markerface','blue')
xlabel('Order')
ylabel('Residual norm')
%prepare_figure

% Plot test results
figure
plot(xtest,[ytmeas yhat3 yhat4 yhat5])
xlabel('Input')
ylabel('Output')
title('Comparison of predictions on test data')
legend({'Actual' ; '3^{rd} order' ; '4^{th} order' ; '5^{th} order'})
%prepare_figure
