function [a,b,cr,cr2]=linreg(x,y)
% LINREG - computes linear regression of y on x
% 
% LINREG(x,y) computes linear regression of y on x, i.e.,  
%                    y = a + b*x
% and returns coefficients a and b with 95% uncertainity. 
% Assumes x,y vectors have no NaN elements, and that x includes
% no measurement error. Uses TINV to estimate student-t value 
% for confidence limits.  See NLINREG when both x and y are 
% considered as random variabiles.
%
% Ref: Draper and Smith, Applied Regression Analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ver 1.0: 11/25/96 (RB)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%	$Revision: 1.6 $
%	$Date: 2008/01/19 18:45:58 $
%	$Id: linreg.m,v 1.6 2008/01/19 18:45:58 duncombe Exp duncombe $
%

% convert to row vectors
[N,M]=size(x);
if M==1
  x=x';y=y';
end

% 95% confidence limit student-t distribution
L=length(x);
P=.95;
p=.5.*(1 + P);
st=tinv(p,N-2);

mx=mean(x);
my=mean(y);
dx=x-mx;
dy=y-my;
reg=sum(dx.*dy)./sum(dx.^2);
b=sum(dx.*y)./sum(dx.^2);
a=my-b.*mx;
syx=(sum(dy.^2)-(sum(dx.*dy)^2)./sum(dx.^2))./(N-2);
syx=sqrt(syx);
aerr=st.*syx.*sqrt(((mx.^2)./sum(dx.^2)+1./N));
berr=st.*syx./sqrt(sum(dx.^2));
corr=corrcoef(x,y);
stdy=std(y-a-b*x);
%Standard error of the estimate CAF
sse=(y-(a+b*x)).^2;
se=sqrt(sum(sse)/(N-2));
%
A=[a aerr b berr corr(1,2)];
disp(['mean x = ',num2str(mx),';   mean y = ',num2str(my)])
disp(['std(y-a-bx) = ',num2str(stdy),';   N = ',int2str(N)])
disp(['a   aerr  = ',num2str(a),'  +-  ',num2str(aerr)])
disp(['b   berr  = ',num2str(b),'  +-  ',num2str(berr)])
disp(['corr coef = ',num2str(corr(1,2))])
disp(['(corr coef)^2 = ',num2str(corr(1,2)^2)])
%
disp(['se = ',num2str(se)])
%
cr=corr(1,2);
cr2=corr(1,2)^2;
%
xr=[min(x) max(x)];
yp=a+b*xr;
ypu=a+aerr + (b+berr)*xr;
ypl=a-aerr + (b-berr).*xr;
plot(x,y,'ro',xr,yp,xr,ypu,'r:',xr,ypl,'r:')
xlabel('X')
ylabel('Y')
title('Least-squares fit of Y on X with maximum 95% CIs')
gtext({'offset= ',num2str(a),'slope= ',num2str(b)})
gtext({'corr coef = ',num2str(corr(1,2)),'(corr coef)^2 = ', ...
		num2str(corr(1,2)^2)})
gtext({'se = ', num2str(se)})


