%  demo_ls.m
%
%  Demonstration of the use of (discrete) linear least squares
%  in polynomial regression.

%  Set up data vectors.

  tvec = [1 2 2 3 3 4 5]';
  yvec = [1 2 3 5 3 7 10]';
  N = length(tvec);
  
%  Set up design matrix corresponding to degree 1 polynomial approximation.

  onevec = ones(N,1);
  A = [onevec tvec];
  
%  Display design matrix and data vector.

  fprintf(' Design matrix and data vector.\n')
  A
  yvec
  
%  Solve for the coefficients of the best least squares polynomial 
%  using the pseudoinverse.

  coef = pinv(A) * yvec;

%  Display the coefficients.

  fprintf(' Coefficients of best least squares polynomial fit.\n')
  coef
  
%  Plot the data and the polynomial fit.

  tmin = min(tvec);
  tmax = max(tvec);
  delta_t = (tmax - tmin) / 1000;
  t = [tmin:delta_t:tmax]';
  poly = coef(1)*t.^0 + coef(2)*t.^1;
  
  plot(tvec,yvec,'*', t,poly,'-')
%  Put legend in upper left corner.
  legend('Data Points','Polynomial Fit','Location','NorthWest');   
  xlabel('t axis'), ylabel('y axis')
  title('Least squares polynomial fit to data')