close all     % closes all windows opened during the previous runs
clear all     % clears the memory from the results of previous runs
clc           % clears screen
disp('Finding the linear best fit y=ax+b using the least square')
disp('approximation (this assumes that x values are exact while y')
disp('values have some errors)')
disp(' ')
disp('Scattered data points (x,y):')
x=[0,.5,1]      
y=[0,0,1]
plot(x,y,'o')   %scattered data plot
xlabel('x')
ylabel('y')
axis('equal')
axis([-.5 1.5 -.5 1.5])
hold on;
pause
disp('Least square fit. Need to compute the following quantities:')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% you may wish to memorise these useful formulae as you
%%%%%%% can use them to compute best fit line by hand
xm=mean(x)       %these four coefficients are necessary to compute
ym=mean(y)       %coefficients for the best straight line fit
xxm=mean(x.^2)
xym=mean(x.*y)
pause
disp('Then the slope of the least square fit is')
a=(xym-xm*ym)/(xxm-xm^2)
pause
disp('and the intercept of the straight line is')
b=ym-a*xm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1=linspace(-0.5,1.5);
y1=a*x1+b;
plot(x1,y1)              % plots the least squre fit; it could
pause                    % be obtained by calling the polyfit function

disp('Now assume that both x and y values have errors. The best fit')
disp('line is different now. To find its coefficients:')
disp('compute the data covariance matrix:')
k=cov(x,y)           % Note that this command will automatically 
                     % subtract the mean from x and y vectors.
                     % The normalisation is 1/(m-1)
%k=cov(x,y,1)        % This will produce the covariance matrix 
                     % normalised with 1/m. See that this would change
                     % the eigenvalues but not the eigenvectors
pause
disp('Find its eigenvectors v and eigenvalues d:')
[v,d]=eig(k)         
t=linspace(-2,2);
pause
disp('Now plot two straight lines passing through the point')
disp('(xm,ym) and containing the eigenvectors:')
disp(' ')
plot(mean(x)+v(1,2)*t,mean(y)+v(2,2)*t,'g') 
disp('This one corresponds to the larger eigenvalue and is')
disp('the best fit.')
disp(' ')
pause 
plot(mean(x)+v(1,1)*t,mean(y)+v(2,1)*t,'r') 
disp('This one corresponds to the smaller eigenvalue and is')
disp('perpendicular to the best fit. Data scattering in this')
disp('direction is minimal')
pause
disp(' ')
disp('The slopes of these lines are:')
slp1=v(2,2)/v(1,2)
slp2=v(2,1)/v(1,1)
pause
disp('Recollect that we derived using geometrical considerations')  
disp('that the slope a of the best fit line is given by the one of')
disp('the roots of the quadratic equation A*a^2+B*a-A=0, where')
yym=mean(y.*y)
A=xym-xm*ym
B=xxm-yym+ym^2-xm^2
pause
disp('i.e the slopes are:')
aslp1=-B/2/A+sqrt(1+B^2/4/A^2)
aslp2=-B/2/A-sqrt(1+B^2/4/A^2)
pause
disp('Observe that they are identical to')
slp1 %and
slp2
pause
disp('Thus when both x and y have errors, the covariance matrix') 
disp('approach leads to the correct best fit line.')