% This generates a least squares fit from n data points
% David Bloore, July 2011

clc;
clear all;

% fixed user inputs

n=input('Number of Data Points?');
xb=zeros(1,n);
b=zeros(1,n);

% data point input

counter=1;

while counter<n+1
   display(sprintf('For data point %d',counter));
   xb(counter)=input('Value for x?');
   b(counter)=input('Value for y?');
   counter=counter+1;
end

% degree of polynomial fit input

boolean=0;

while boolean==0;
    d=input('Desired degree of polynomial fit?');
    if d<n;
        boolean=1;
    else
        clc;
        display(sprintf('Please enter a degree less than %d',n));
    end
end

% construction of A

A=zeros(n,d+1);
A(1:n,1)=ones(1,n);

counter=1;

while counter<d+1;
    xtemp=xb.^counter;
    A(1:n,counter+1)=xtemp;
    counter=counter+1;
end

% solution of A'*A*u=A'*b

AtA=A'*A;
Atb=A'*b';
u=AtA\Atb;

% construction of fit funtion domain

x1p=(11*xb(1)-xb(n))/10;
xnp=(11*xb(n)-xb(1))/10;
x=linspace(x1p,xnp,10000);
y=u(1)*ones(1,10000);

% evaluation of fit function

counter=1;

while counter<d+1;
    y = y + u(counter+1)*(x.^counter);
    counter=counter+1;
end

plot(xb,b,'o',x,y,'--');
title(sprintf('Degree %d Least Squares Fit of %d Data Points',d,n));
legend('Data Points',sprintf('Fit of Degree %d',d));
xlim([x1p xnp]);
xlabel('x');
ylabel('y');
grid on;