linreg1

 LINREG1 solves directly for the coefficients in a linear model

 LINREG1(x,y,stdx,stdy,b0) solves directly for the coefficients in
 the linear model 
          y = a + b*x
 for the case when both x and y are considered independent random 
 variables with zero-mean noise characterized by stdx and stdy, 
 which are independent of specific observations. Inputs are:
 the data xi,yi; the assumed known stdx, stdy; and the initial 
 guess of the slope b0. Code solves nonlinear problem exactly for 
 the slope b which minimizes the merit function
         S = sum((yi-a-b*xi).^2/(stdy^2 + b^2*stdx^2))
 which represents the sum of the distances normal to the least
 squares line. Then a is determined from a=ym-b*xm. Approximate 
 95% confidence intervals given for a,b.  NOTE 1: solution reduces 
 to linear regression result (from LINREG) when the ratio of 
 stdx/stdy gets very small or very large.  NOTE 2: for the case 
 of very small stdx, then substitution of stdydata for stdy 
 gives ci's for a in close agreement with LINREG.  NOTE 3: code 
 uses FZERO to find value of b near b0. If b0 is poor quess, 
 then b may be a poor solution. Check using plot of results 
 and experimentation with diffeent b0 that b given is robust.
 (See LINREG for usual case of regressing y on x, when x has 
 no measurement error.)

 Solution approach follows Reed (1992), Am. J. Phys., 60, 59-62, 
 and Reed (1989), Am. J. Phys., 57, 642-646, who corrected earlier 
 formulas given by York (1966), Can. J. Phys., 44, 1079-1086.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ver 1.0: 11/25/96 (RB)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001   function linreg1(x,y,stdx,stdy,b0)
0002 
0003 % LINREG1 solves directly for the coefficients in a linear model
0004 %
0005 % LINREG1(x,y,stdx,stdy,b0) solves directly for the coefficients in
0006 % the linear model
0007 %          y = a + b*x
0008 % for the case when both x and y are considered independent random
0009 % variables with zero-mean noise characterized by stdx and stdy,
0010 % which are independent of specific observations. Inputs are:
0011 % the data xi,yi; the assumed known stdx, stdy; and the initial
0012 % guess of the slope b0. Code solves nonlinear problem exactly for
0013 % the slope b which minimizes the merit function
0014 %         S = sum((yi-a-b*xi).^2/(stdy^2 + b^2*stdx^2))
0015 % which represents the sum of the distances normal to the least
0016 % squares line. Then a is determined from a=ym-b*xm. Approximate
0017 % 95% confidence intervals given for a,b.  NOTE 1: solution reduces
0018 % to linear regression result (from LINREG) when the ratio of
0019 % stdx/stdy gets very small or very large.  NOTE 2: for the case
0020 % of very small stdx, then substitution of stdydata for stdy
0021 % gives ci's for a in close agreement with LINREG.  NOTE 3: code
0022 % uses FZERO to find value of b near b0. If b0 is poor quess,
0023 % then b may be a poor solution. Check using plot of results
0024 % and experimentation with diffeent b0 that b given is robust.
0025 % (See LINREG for usual case of regressing y on x, when x has
0026 % no measurement error.)
0027 %
0028 % Solution approach follows Reed (1992), Am. J. Phys., 60, 59-62,
0029 % and Reed (1989), Am. J. Phys., 57, 642-646, who corrected earlier
0030 % formulas given by York (1966), Can. J. Phys., 44, 1079-1086.
0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 % ver 1.0: 11/25/96 (RB)
0033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0034 
0035 N=length(x);
0036 x=x(:);
0037 y=y(:);
0038 
0039 % compute coefficients
0040 vary=stdy^2;
0041 varx=stdx^2;
0042 r2=vary/varx;
0043 xm=mean(x);
0044 ym=mean(y);
0045 U=x-xm;
0046 V=y-ym;
0047 UU=U'*U;
0048 VV=V'*V;
0049 UV=U'*V;
0050 
0051 % solve for b, a
0052 global B C
0053 B=(r2*UU-VV)./UV;
0054 C=-r2;
0055 b=fzero('quadrat',b0,0.00001); 
0056 a=ym-b*xm;
0057 cc=corrcoef(x,y);
0058 
0059 % compute yp, xp, Rs, S
0060 yp=a+b*x;
0061 r=yp-y;
0062 xp=(y-a)./b;
0063 b2=b^2;
0064 varydata=(std(yp-y))^2;
0065 varxdata=(std(xp-x))^2;
0066 r2data=varydata/varxdata;
0067 S=((V-b*U)'*(V-b*U))/(vary + b2*varx);
0068 
0069 % compute confidence intervals for a,b
0070 A=UV;
0071 B=r2*UU-VV;
0072 C=-r2*UV;
0073 SD=(S/(N-2))/(2*A*b+B)^2;
0074 b2r2=(b2-r2);
0075 br2=b*r2;
0076 
0077 % start with b
0078 am=(b2r2)^2;
0079 bm=(2*br2)^2;
0080 cm=4*br2*b2r2;
0081 varbx=varx*(am*VV+bm*UU+cm*UV);
0082 am=4*b2;
0083 bm=(b2r2)^2;
0084 cm=(-4*b*b2r2);
0085 varby=vary*(am*VV+bm*UU+cm*UV);
0086 varb=SD*(varbx+varby);
0087 stdb=sqrt(varb);
0088 
0089 % now do a
0090 XD=xm^2/(2*A*b+B)^2;
0091 am=(2*br2)^2;
0092 bm=b2r2^2;
0093 cm=4*br2*b2r2;
0094 varax=b2/N+XD*(am*UU+bm*VV+cm*UV);
0095 am=(2*b)^2;
0096 cm=(4*b)*b2r2;
0097 varay=1/N+XD*(bm*UU+am*VV-cm*UV);
0098 vara=varx*varax+vary*varay;
0099 stda=sqrt(vara);
0100 
0101 % now compute 95% ci's
0102 st=tinv(0.975,N-2);
0103 aerr=st*stda;
0104 berr=st*stdb;
0105 
0106 % disp results
0107 disp(['mean x = ',num2str(xm),';  mean y = ',num2str(ym)])
0108 disp(['std(y-a-bx) = ',num2str(std(y-a-b*x)),';  N = ',int2str(N)])
0109 disp(['a   aerr  = ',num2str(a),'  +-  ',num2str(aerr)])
0110 disp(['b   berr  = ',num2str(b),'  +-  ',num2str(berr)])
0111 disp(['corr coef = ',num2str(cc(1,2))])
0112 disp(['(corr coef)^2 = ',num2str(cc(1,2)^2)])
0113 disp(['stdx = ',num2str(stdx),';  stdxdata = ',num2str(sqrt(varxdata))])
0114 disp(['stdy = ',num2str(stdy),';  stdydata = ',num2str(sqrt(varydata))])
0115 disp(['vary/varx = ',num2str(r2),'; varydata/varxdata = ',num2str(r2data)])
0116 
0117 xr=[min(x) max(x)];
0118 yp=a+b*xr;
0119 ypu=a+aerr + (b+berr)*xr;
0120 ypl=a-aerr + (b-berr).*xr;
0121 plot(x,y,'o',xr,yp,'g-',xr,ypu,'r:',xr,ypl,'r:')
0122 xlabel('X')
0123 ylabel('Y')
0124 title(['Least-squares nonlinear fit of Y on X with maximum 95% CIs for stdy/stdx = ',num2str(sqrt(r2))])
0125 end
0126