function linearfit
   t = [ 1 2 3 4 5 6 7 8 9 ]
   y = [ 1 2 1 3 4 3 4 6 6 ]

   [alpha,beta] = lsq_linear_fit(t,y)

   f = @(t) alpha + beta*t
   hold off
   plot(t,y,'ro')
   range = max(t)-min(t);
   hold on
   myfplot(f,[min(t)-range/5,max(t)+range/5])
end

function [alpha,beta] = lsq_linear_fit(t,y) 
                         % Least squares straight line fit
                         % by solving the Normal Equations
   [q,m] = size(t);
   a = [ones(m,1) t']
   b = y'
   [alpha,beta] = solve_normal_eqns(a,b)
end

function [alpha,beta]=solve_normal_eqns(a,b)
   xbar = ((a'*a)\(a'*b))';
   alpha = xbar(1);
   beta  = xbar(2);
end

function myfplot(f,interval)
	x = linspace(interval(1),interval(2),200);
	y = f(x);
	plot(x,y,'b');
end