linreg

 LINREG computes linear regression of y on x, i.e.

 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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001    function linreg(x,y)
0002 
0003 % LINREG computes linear regression of y on x, i.e.
0004 %
0005 % LINREG(x,y) computes linear regression of y on x, i.e.,
0006 %                    y = a + b*x
0007 % and returns coefficients a and b with 95% uncertainity.
0008 % Assumes x,y vectors have no NaN elements, and that x includes
0009 % no measurement error. Uses TINV to estimate student-t value
0010 % for confidence limits.  See NLINREG when both x and y are
0011 % considered as random variabiles.
0012 %
0013 % Ref: Draper and Smith, Applied Regression Analysis
0014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0015 % ver 1.0: 11/25/96 (RB)
0016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0017 
0018 % convert to row vectors
0019 [N,M]=size(x);
0020 if M==1
0021   x=x';y=y';
0022 end
0023 
0024 % 95% confidence limit student-t distribution
0025 L=length(x);
0026 P=.95;
0027 p=.5.*(1 + P);
0028 st=tinv(p,N-2);
0029 
0030 mx=mean(x);
0031 my=mean(y);
0032 dx=x-mx;
0033 dy=y-my;
0034 reg=sum(dx.*dy)./sum(dx.^2);
0035 b=sum(dx.*y)./sum(dx.^2);
0036 a=my-b.*mx;
0037 syx=(sum(dy.^2)-(sum(dx.*dy)^2)./sum(dx.^2))./(N-2);
0038 syx=sqrt(syx);
0039 aerr=st.*syx.*sqrt(((mx.^2)./sum(dx.^2)+1./N));
0040 berr=st.*syx./sqrt(sum(dx.^2));
0041 corr=corrcoef(x,y);
0042 stdy=std(y-a-b*x);
0043 A=[a aerr b berr corr(1,2)];
0044 disp(['mean x = ',num2str(mx),';   mean y = ',num2str(my)])
0045 disp(['std(y-a-bx) = ',num2str(stdy),';   N = ',int2str(N)])
0046 disp(['a   aerr  = ',num2str(a),'  +-  ',num2str(aerr)])
0047 disp(['b   berr  = ',num2str(b),'  +-  ',num2str(berr)])
0048 disp(['corr coef = ',num2str(corr(1,2))])
0049 disp(['(corr coef)^2 = ',num2str(corr(1,2)^2)])
0050 
0051 xr=[min(x) max(x)];
0052 yp=a+b*xr;
0053 ypu=a+aerr + (b+berr)*xr;
0054 ypl=a-aerr + (b-berr).*xr;
0055 plot(x,y,'o',xr,yp,xr,ypu,'r:',xr,ypl,'r:')
0056 xlabel('X')
0057 ylabel('Y')
0058 title('Least-squares fit of Y on X with maximum 95% CIs')
0059 end
0060