%FUNCTION: find the linear regression estimates by the least squares method
function out=lrest(x,y,order,show);

% INPUT: x==source values
%        y==target values
%        order == the linear function order
%        show == plot the linear line 
% OUTPUT: out==[b0 b1 de], for the linear mapping function y=b0+b1*x
%         sse==standard deviation for erroe 
%         2.0*db0 95% confidence interval for b0 ( +-t standard error estimate)
%         2.0*db1 95% confidence interval for b1
%%%

if nargin<3
   order=1;
   show=0;
elseif nargin<4
   show=0;
end

n=length(x);

mx=mean(x);
my=mean(y);

sxx=sum((x-mx).^2);
sxy=sum((x-mx).*(y-my));

if order==0

   b1=1;
   b0=my-mx;

elseif order==1
   if sxx==0
      b1=1;
   else
      b1=sxy/sxx;   %estimate for b1
   end

   b0=my-b1*mx;  %estimate for b0
end

sse=sum((y-b0-b1*x).^2);             %sum of squares for error
de=sqrt(sse/(n-2));         %standard deviation for error
%db0=de*sqrt( sum(x.^2)/n/sxx );  %standard deviation for b0
%db1=de/sqrt(sxx);             %standard deviation for b1

out=[b0 b1 de];


if show==1
    plot(x,y,'y*',x,b0+b1*x,'r-')
    title(['Linear Regression Estimates (Std= ' num2str(de) ')']);
    xlabel('Soucre');
    ylabel('target');
    drawnow;
end