% Demonstrates use of polyfit to perform least squares 
% polynomial fit to data --- in this case a linear 
% fit.

format long;

% First, define some test data which is precisely linear in x.
%
%   y = a x + b

a  = 2.31;
b = -4.25;
xmin = -1.0;
xmax = 1.0;
nx = 1001;

x = linspace(xmin, xmax, nx);
y = a*x + b;

% Perform the fit; coefficients are returned in pfit. 
% Note that supplying 1 as the 3rd argument means request linear fit.

pfit = polyfit(x, y, 1);

afit = pfit(1);
bfit = pfit(2);

% Should have afit == a, bfit == b (to within machine precision) in this case.

fprintf('tpolyfit: Exactly linear data:\n a=%g\n b=%g\n afit=%g\n bfit=%g\n', ...
   a, b, afit, bfit);

% Now add random values distributed uniformly on 
% interval [-noiseamp .. noiseamp] to the data.

noiseamp = 1.0;
y = y + (-noiseamp + 2.0 * noiseamp * rand(1,nx));

% Fit to the "noisy" data.

pfit = polyfit(x, y, 1);

afit = pfit(1);
bfit = pfit(2);

fprintf('\ntpolyfit: Linear data with noise amplitude %g:\n', noiseamp);  
fprintf(' a=%g\n b=%g\n afit=%g\n bfit=%g\n', a, b, afit, bfit);

% Make plot of noisy data and fit---note how the fit data is generated 
% in place (using vector arithmetic) in the call to plot.

plot(x, y, x, afit*x + bfit, 'ro');
