% lsExample: An example script for least squares regression.
%
% This function computes best fits for two model functions to a small set
% of data.
%
% Variable names are not the standard used in the rest of the course, to
% give practice working with other variable names.
%
% Ian Mitchell for CPSC 303, 2/9/06

% Data.  Abscissae are not ordered, but that doesn't matter.
ti = [ 0.9; -0.5; 0.2; 0.0; -1.0 ];
zi = [ 5.4;  4.0; 2.2; 1.9;  5.2 ];

% Amount of data.
m = length(ti);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Plot the data;
figure;
hi = plot(ti, zi, 'ko');
hold on;
set(hi, 'MarkerSize', 12);
legend_handles = hi(1);
legend_strings = { 'Data Values' };

% Points at which to plot the regression curves.
ts = linspace(1.25 * min(ti), 1.25 * max(ti), 101);

disp('Press any key to fit polynomial.');
pause;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% First, we'll fit to a quadratic polynomial.  Three coefficients.  Use
% "subscript" q to denote the quadratic fit.
nq = 3;
Aq = zeros(m,nq);
% Basis zero is a constant.
Aq(:,1) = ones(m,1);
% Basis one is linear.
Aq(:,2) = ti;
% Basis two is quadratic.
Aq(:,3) = ti.^2;

% Right hand side is just the data values.
zq = zi;

% Solve by normal equations A' * A x = A' * z.  We'll store the coefficients
% in vector x.
xq = (Aq' * Aq) \ (Aq' * zq);

% Plot the fit of the polynomial.  Note that Matlab uses indices starting at one.
zq = xq(1) + xq(2) * ts + xq(3) * ts.^2;
hq = plot(ts, zq, 'b-');
set(hq, 'LineWidth', 4);
legend_handles = [ legend_handles, hq(1) ];
legend_strings = { legend_strings{:}, 'Least squares quadratic fit' };

disp('Press any key to fit hyperbolic cosine.');
pause;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Now, let's try a fitting a nonlinear function.  We noticed that the
% linear component of the fit above was pretty small, so we'll leave it
% out.  For the quadratic component, let's substitute another symmetric
% convex function: cosh(pi * t).  Use "subscript" h to denote a
% hyperbolic trig fit.
nh = 2;
Ah = zeros(m, nh);
% Basis zero is still a constant.
Ah(:,1) = ones(m,1);
% Basis one is the cosh.
Ah(:,2) = cosh(pi * ti);
% There is no basis two.

% Right hand side is still the data values.
zh = zi;

% Now let's solve with Matlab's builtin.  In exact arithmetic, answer would
% be the same as with normal equations.  Computationally, Matlab's method is
% more stable, although for such a small matrix it will not matter.
xh = Ah \ zh;

% Plot the fit of this hyperbolic function.
zh = xh(1) + xh(2) * cosh(pi * ts);
hh = plot(ts, zh, 'r--');
set(hh, 'LineWidth', 4);
legend_handles = [ legend_handles, hh(1) ];
legend_strings = { legend_strings{:}, 'Least squares offset cosh fit' };

disp('Press any key to fit a quartic interpolant.');
pause;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% For comparison, let's look at the quartic interpolant.  Easiest way to
% find it is to use Matlab's buildin polyfit.
pp = polyfit(ti, zi, 4);
zp = polyval(pp, ts);
hp = plot(ts, zp, 'm-.');
set(hp, 'LineWidth', 4);
legend_handles = [ legend_handles, hp(1) ];
legend_strings = { legend_strings{:}, 'Quartic interpolant' };

disp('Press any key to fit a cubic spline.');
pause;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Finally, let's also compare a cubic spline.  We'll use Matlab's routine
% and the default not-a-knot condition.
ss = spline(ti, zi);
zs = ppval(ss, ts);
hs = plot(ts, zs, 'k:');
set(hs, 'LineWidth', 4);
legend_handles = [ legend_handles, hs(1) ];
legend_strings = { legend_strings{:}, 'Cubic spline interpolant' };


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Finally, pretty up the figure.
title({ 'Two least squares regression fits for a data set'; ...
        'Ian Mitchell for CPSC 303' });
xlabel('t');
ylabel('z');
legend(legend_handles, legend_strings);