function linregctd(x,y,variable,unit,sensor,data)
% 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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 
% VERSION: 
%	$Revision: 1.5 $
%	$Date: 2008/01/19 18:44:55 $
%	$Id: linregctd.m,v 1.5 2008/01/19 18:44:55 duncombe Exp duncombe $
%
% CHANGELOG:
% 2007/12/17 CMDR 
% * excessive use of the gtext function is driving me insane. test
% and skip these lines

CMDR_insanity_check=0;

% 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=t_inv(p,N-2);

mx=mean(x);
my=mean(y);
%mx=median(x);
%my=median(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(' ')
disp([sensor,' ' ,variable,' (',data,')'])
disp(['mean x = ',num2str(mx),';   mean y = ',num2str(my)])
disp(['std(y-a-bx) = ',num2str(stdy),';   N = ',int2str(N)])
disp(['a   a-err  = ',num2str(a),'  +-  ',num2str(aerr)])
disp(['b   b-err  = ',num2str(b),'  +-  ',num2str(berr)])
disp(['corr coef = ',num2str(corr(1,2))])
disp(['(corr coef)^2 = ',num2str(corr(1,2)^2)])
%
disp(['se = ',num2str(se)])
%
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(['Bottle ',variable,' (',unit,')'])
ylabel(['CTD ',variable,' (',unit,')'])
title(['Least-squares fit of Y on X with maximum 95% CIs'])
if CMDR_insanity_check, gtext({['offset= ',num2str(a),'  +-  ',num2str(aerr)],[' slope= ',num2str(b),'  +-  ',num2str(berr)]}) ; end;
if CMDR_insanity_check, gtext({['corr coef = ',num2str(corr(1,2))],['(corr coef)^2 = ',num2str(corr(1,2)^2)]}) ; end;
if CMDR_insanity_check, gtext({['se = ', num2str(se)]}) ; end;
grid on

