% ch13_polynomial_models.m
% Companion code for Chapter 13: Optimization and Data Fitting.
%
% Compare polynomial models of different degrees using training and
% independent validation data. A smaller training error does not necessarily
% mean better predictions. Change the data and degrees below.
% Uses MATLAB's polyfit with centering and scaling for improved conditioning.
% No fitting toolbox is required. All fits here are unweighted.

clear; close all; clc;

%% Example data and settings
seed = 7;
rng(seed, 'twister');
% Alternative: comment out the line above and uncomment the following.
% rng('shuffle');

true_function = @(x) sin(2*pi*x) + 0.5*x;
noise_sigma = 0.18;
x = linspace(0,1,18)';
y = true_function(x) + noise_sigma*randn(size(x));

% Independent measurements, not used to fit any coefficients.
x_validation = linspace(0,1,101)';
y_validation = true_function(x_validation) ...
    + noise_sigma*randn(size(x_validation));

degrees = [1, 5, 14];
make_plot = true;
show_true_function = true; % Set false for measured data without a known model
plot_x = [0, 1];
plot_y = [-1.5, 2];        % Use [] for automatic vertical limits

% For measured data, replace x,y and x_validation,y_validation above with
% separate training and validation sets. Keep validation data out of the fit.
% Do not repeatedly tune models against a final test set: reserve fresh test
% data if an unbiased assessment after model selection is needed.

%% Fit each degree and compare predictions
validateattributes(x,{'double'},{'real','vector','finite','nonempty'});
validateattributes(y,{'double'},{'real','vector','finite','nonempty'});
validateattributes(x_validation,{'double'},{'real','vector','finite','nonempty'});
validateattributes(y_validation,{'double'},{'real','vector','finite','nonempty'});
x=x(:); y=y(:); x_validation=x_validation(:); y_validation=y_validation(:);
if numel(x)~=numel(y) || numel(x_validation)~=numel(y_validation)
    error('Each data set must contain matching x and y vectors.');
end
validateattributes(degrees,{'double'},{'real','vector','integer','nonnegative','finite','nonempty'});
if any(degrees >= numel(unique(x)))
    error('Each degree must be smaller than the number of distinct training x values.');
end

n_models = numel(degrees);
models = cell(n_models,1);
training_rmse = zeros(n_models,1);
validation_rmse = zeros(n_models,1);
xx = linspace(plot_x(1),plot_x(2),600)';
curves = zeros(numel(xx),n_models);
residuals = zeros(numel(x),n_models);

fprintf('\nPolynomial model comparison\n');
fprintf('%8s %18s %18s\n','Degree','Training RMSE','Validation RMSE');
for k=1:n_models
    [p,S,mu] = polyfit(x,y,degrees(k));
    fitted = polyval(p,x,[],mu);
    predicted = polyval(p,x_validation,[],mu);
    residuals(:,k) = y-fitted;
    training_rmse(k) = sqrt(mean(residuals(:,k).^2));
    validation_rmse(k) = sqrt(mean((y_validation-predicted).^2));
    curves(:,k) = polyval(p,xx,[],mu);
    % Coefficients are descending powers of z = (x-mu(1))/mu(2).
    models{k} = struct('degree',degrees(k),'coefficients',p,'mu',mu,'fit_info',S);
    fprintf('%8d %18.8g %18.8g\n',degrees(k),training_rmse(k),validation_rmse(k));
end
fprintf('RMSE uses N in the denominator; it is a prediction-error measure.\n');
fprintf('A single validation set does not establish the best model in general.\n');
% Predict with model k using:
% y_new = polyval(models{k}.coefficients, x_new, [], models{k}.mu);

%% Optional fitted curves and training residuals
if make_plot
    styles = {'k:','k-','k-.','k--'};
    labels = arrayfun(@(d) sprintf('Degree %d',d),degrees,'UniformOutput',false);
    figure('Color','w');
    subplot(2,1,1); hold on;
    handles = gobjects(0);
    names = {};
    handles(end+1) = plot(x,y,'ko','MarkerFaceColor','w','MarkerSize',6);
    names{end+1} = 'Training data';
    if show_true_function
        handles(end+1) = plot(xx,true_function(xx),'k--','LineWidth',1.2);
        names{end+1} = 'True function';
    end
    for k=1:n_models
        handles(end+1) = plot(xx,curves(:,k),styles{mod(k-1,numel(styles))+1},'LineWidth',1.5);
        names{end+1} = labels{k};
    end
    xlabel('x','FontSize',14); ylabel('y','FontSize',14);
    title('Polynomial model comparison','FontSize',14,'FontWeight','normal');
    legend(handles,names,'Location','best');
    xlim(plot_x);
    if ~isempty(plot_y), ylim(plot_y); end
    set(gca,'FontSize',14); grid on; box on;

    subplot(2,1,2); hold on;
    h = gobjects(n_models,1);
    for k=1:n_models
        h(k) = plot(x,residuals(:,k),styles{mod(k-1,numel(styles))+1},'LineWidth',1.2);
    end
    plot(plot_x,[0 0],'k-','LineWidth',0.5);
    xlabel('x','FontSize',14); ylabel('Training residual','FontSize',14);
    legend(h,labels,'Location','best');
    xlim(plot_x); set(gca,'FontSize',14); grid on; box on;
end