%%%%%%%%%%%%%%5
%  Matlab Script to demonstrate curve fitting by  Least squares solutions of
% A'*A=A'*b  (and A' is A transpose in matlabese)
%%%%%%%%%%%%%%
close all
clear all

%%%%%
%   A simple problem from class to fit 5 points with a straight line
%    of form  y=c(1)+c(2)*x.
%       c(1), c(2)  are the coefficients of the line and are
%    the  Least Squares solution of  
%             Ac=y_p  
%    where  A is the Vandermonde matrix [ 1 x_p ] 
%     here x_p is the array of x coordinates for all the points
%     y_p is  the array of y coordinates for the points
%%%%%
pnts=[ -5 -4 ;  -2 -1;   0   0;  2  2; 3 4.5]  % an array of 5 X-Y points
xaxis_label='X';
yaxis_label='Y';
%%
%  for more fun, uncomment the next 3 lines line to get
%  X-midterm grades vs Y average homework grades
%
% pnts=load('Midterm_vs_Homework.txt')
% xaxis_label='Midterm Grades';
% yaxis_label='Average Homework Grade';
%
%%


Xp=pnts(:,1);  % The x array
Yp=pnts(:,2);  % The y array
%%%%
%  plot the data points
%%%%
figure
plot(Xp,Yp,'o','MarkerSize',[10],'MarkerFaceColor','b'); 
xlabel(xaxis_label);
ylabel(yaxis_label); 
grid

%%%%%
%  now create the vandermonde matrix A= [ 1 x ]
%%%%%
A=[ ones(size(Xp)) Xp ]

%%%%
% and solve the Least Squares Problem
%%%%

ATA=A'*A  % form A transpose A
ATb=A'*Yp  % form A transpose b

c=ATA\ATb  % and solve for the coefficients c by elimination
                   % (matlab style)

%%%%
%  now form the best fit line and plot it
%%%%      
%
%  first set up an array for x from min(Xp) to max(Xp) with npnts pnts
%
%
xmin=min(Xp);  
xmax=max(Xp);
npnts=20;
dx=(xmax-xmin)/npnts; % x array spacing
x=xmin:dx:xmax;  % the x array

%
% calculate the line y=c(1)+c(2)*x
%
y=c(1)+c(2)*x;

hold on
hl=plot(x,y,'r','linewidth',[2]); %plot the line

%%%% 
%  Now calculate the best fit quadratic by using the VanderMonde Matrix
%  A=[ 1 x x^2]
%%%%

A=[ ones(size(Xp)) Xp Xp.^2 ]

%%%%
% and solve the Least Squares Problem
%%%%

ATA=A'*A  % form A transpose A
ATb=A'*Yp  % form A transpose b

c=ATA\ATb  % and solve for the coefficients c

%%%%
%  now form the best fit quadratic y=c(1)+c(2)*x+c(3)*x^2
%       and plot it;
%%%%

y=c(1)+x.*(c(2)+x*c(3));  % this is actually the correct way to  evaluate a polynomial;

hold on
hq=plot(x,y,'g','linewidth',[2]);


%%%
%  and for good luck, let's get the least squares cubic with a QR decomposition
%%%

A=[ ones(size(Xp)) Xp Xp.^2 Xp.^3 ]; % VanderMonde matrix of order 3
[Q,R]=qr(A,0);  % matlab's economy QR decompostion
c=R\(Q'*Yp)    % solve Rc=Q'*yb by Elimination (R is upper
                      % triangular so it's a cheap backsubstition)

%
%   calculate the curve and plot
%

y=c(1)+x.*(c(2)+x.*(c(3)+x*c(4))); 
hc=plot(x,y,'k','linewidth',[2]);

legend([hl hq hc],'Least squares line','Least Squares quadratic','Least Squares cubic',0);