%  This script fits a 20th order polynomial to noisy observations of 
%     f(x) = sin(x) + 2*exp( - 30 * x^2)
%
%  Figure caption:  Data simulated from function f(x) = sin(x) +
%  2exp(-30*x*x) together with twentieth order polynomial fit
%  (shown as line).  Note that the polynomial is over-fitting
%  (under-smoothing) in the relatively smooth regions of f(x), and
%  under-fitting (over-smoothing) in the peak.  (In the data shown
%  here, the noise standard deviation is 1/50 times the standard
%  deviation of the function values.)

nObservations = 101;
polynomialOrder = 20;

%  Set up observations and function values
xAtObservations = -2:0.04:2;
fAtObservations = sin(xAtObservations) + 2*exp(-30*(xAtObservations.^2));

%  Obtain a reasonable noise level and add noise.  
fsd = std(fAtObservations);
addedNoiseSD = fsd/50;
addedNoise = normrnd(0, addedNoiseSD, 1, 101);
observations = fAtObservations + addedNoise;

%  Fit polynomial model.
polynomialModel = polyfit(xAtObservations, observations, 20);
modelsPredictions = polyval(polynomialModel, xAtObservations);

%  Plot the model's predictions.
plot(xAtObservations, modelsPredictions, '-k', 'LineWidth', 2)

hold on;
%  Add the actual observations.
plot(xAtObservations, observations, 'ok', ...
     'MarkerSize', 6, 'MarkerFaceColor', 'k')

%  Add axes and legends.
set(gca, 'FontSize', 22, ...
         'XLim', [-2.1, 2.1], 'YLim', [-1.1, 2.1], ...
         'XTick', -2:1:2, 'YTick', -1:0.5:2)

xlabel({'';'x'}, 'FontSize', 22)
ylabel({'y';''}, 'FontSize', 22)
set(gcf, 'Position', [200, 100, 1000, 800])