% Least squares polynomial fitting, L. Braile, September 2002
% Calculate a 5th order polynomial
% 
at = [-6.4917 11.9745 -5.7484 1.2424 -0.1218 0.0045];
% Calculate data
xd = [0.5:0.25:9.5]';
% set up the function of the form y = a0 + a1*x +a2*x^2 ...
% the ones are used for the constant (a0) term
Xd = [ones(size(xd)) xd xd.^2 xd.^3 xd.^4 xd.^5];
yd = Xd*at';
% add random noise to y data
randn('state',0);  %set state so random numbers will be same each run
noise = 0.5*randn(max(size(xd)),1);
yd = yd + noise;
% calculate theoretical data (correct 5th order polynomial without noise)
xt = [0.5:0.1:9.5]';
Xt = [ones(size(xt)) xt xt.^2 xt.^3 xt.^4 xt.^5];
yt = Xt*at';
% use Matlab function polyfit to determine least squares fit to data
[p,s] = polyfit(xd,yd,5);  % last input is order of the polynomial
% evaluate polynomial at the xt locations
[yt1,del2] = polyval(p,xt,s);
% plot the data, theoretical line, least squares fit and 95% error bound on y
plot(xd,yd,'or',xt,yt,'--',xt,yt1,'-',xt,yt1+2*del2,':',xt,yt1-2*del2,':')
xlabel('x-value')
ylabel('y-value')
title('Least squares fit to data with noise') 
text(1,-2,'o = data with noise (std = 0.5); long dash = data with no noise') 
text(1,-3,'line = least squares fit; short dash lines = 95% confidence limits on y')
% calculate statistics, sum of squared residuals
yd1 = polyval(p,xd);
sumsq = sum((yd-yd1).^2)