Home > atmlab > math > wtls_line.m

wtls_line

PURPOSE ^

weighted total least squares (wtls) fit of a straigth line

SYNOPSIS ^

function [a,b,alpha,p,chiopt,Cab,Calphap]=wtls_line(xin,yin,uxin,uyin,guck);

DESCRIPTION ^

 weighted total least squares (wtls) fit of a straigth line 
 to a set of points with uncertainties in both coordinates

 FORMAT    [a,b,alpha,p,chiopt,Cab,Calphap]=...
            wtls_line(xin,yin,uxin,uyin,guck)

 OUT:          a, b    usual straight line parameters
                       y=a*x+b
               alpha,p more stable parametrisation
                       y*cos(alpha)-x*sin(alpha)-p=0
                       alpha: slope angle in radians
                       p: distance of straight line from (0,0)
                       conversion: a=tan(alpha),b=p/cos(alpha)
               chiopt  minimum chisquare found
               Cab     covariances, [var(a),var(b),cov(a,b)]
               Calphap covariances, [var(alpha),var(p),cov(alpha,p)]

  IN:          xin     abscissa vector
               yin     ordinate vector
               uxin    (standard) uncertainties of xin, same size as xin
               uyin    (standard) uncertainties of yin, same size as yin
               guck    flag, if >0, chisquare is plotted over whole
                       alpha-interval, optimum alpha is shown as red +

 algorithm:    M. Krystek & M. Anton
               Physikalisch-Technische Bundesanstalt Braunschweig, Germany
               Meas. Sci. Tech. 18 (2007), pp3438-3442

 tested for Matlab 6 and Matlab 7
 testdata: script file pearson_york_testdata.m

 2007-03-08

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

wtls_line.m

SOURCE CODE ^

0001 % weighted total least squares (wtls) fit of a straigth line
0002 % to a set of points with uncertainties in both coordinates
0003 %
0004 % FORMAT    [a,b,alpha,p,chiopt,Cab,Calphap]=...
0005 %            wtls_line(xin,yin,uxin,uyin,guck)
0006 %
0007 % OUT:          a, b    usual straight line parameters
0008 %                       y=a*x+b
0009 %               alpha,p more stable parametrisation
0010 %                       y*cos(alpha)-x*sin(alpha)-p=0
0011 %                       alpha: slope angle in radians
0012 %                       p: distance of straight line from (0,0)
0013 %                       conversion: a=tan(alpha),b=p/cos(alpha)
0014 %               chiopt  minimum chisquare found
0015 %               Cab     covariances, [var(a),var(b),cov(a,b)]
0016 %               Calphap covariances, [var(alpha),var(p),cov(alpha,p)]
0017 %
0018 %  IN:          xin     abscissa vector
0019 %               yin     ordinate vector
0020 %               uxin    (standard) uncertainties of xin, same size as xin
0021 %               uyin    (standard) uncertainties of yin, same size as yin
0022 %               guck    flag, if >0, chisquare is plotted over whole
0023 %                       alpha-interval, optimum alpha is shown as red +
0024 %
0025 % algorithm:    M. Krystek & M. Anton
0026 %               Physikalisch-Technische Bundesanstalt Braunschweig, Germany
0027 %               Meas. Sci. Tech. 18 (2007), pp3438-3442
0028 %
0029 % tested for Matlab 6 and Matlab 7
0030 % testdata: script file pearson_york_testdata.m
0031 %
0032 % 2007-03-08
0033 
0034 function [a,b,alpha,p,chiopt,Cab,Calphap]=...
0035     wtls_line(xin,yin,uxin,uyin,guck);
0036 
0037 tol=1e-6; %"tolerance" parameter of fnimbnd, see there
0038 
0039 global ux uy x y
0040 % inputs
0041 if nargin<5,guck=0;end
0042 % force column vectors
0043 x=xin(:);
0044 y=yin(:);
0045 ux=uxin(:);
0046 uy=uyin(:);
0047 % "initial guess":
0048 p0=polyfit(x,y,1);
0049 alpha0=atan(p0(1));
0050 % one-dimensional search, use p=p^
0051 [alphaopt,chiopt,exitflag] = ...
0052     fminbnd(@chialpha,alpha0-pi/2,alpha0+pi/2,optimset('TolX',tol));
0053 % get optimum p from alphaopt
0054 alpha=alphaopt;
0055 uk2=ux.^2*sin(alpha)^2+uy.^2*cos(alpha)^2;
0056 u2=1./mean(1./uk2);
0057 w=u2./uk2;
0058 xbar=mean(w.*x);
0059 ybar=mean(w.*y);
0060 p=ybar*cos(alpha)-xbar*sin(alpha);
0061 % convert to a, b parameters y=a*x+b
0062 a=sin(alpha)/cos(alpha);
0063 b=p/cos(alpha);
0064 % --- uncertainty calculation, covariance matrix = 2*inv(Hessian(chi2)) ---
0065 n=length(x);
0066 vk=y.*cos(alpha)-x*sin(alpha)-p;
0067 vka=-y*sin(alpha)-x*cos(alpha);
0068 vkaa=-vk-p;
0069 fk=vk.*vk;
0070 fka=2*vk.*vka;
0071 fkaa=2*(vka.^2+vk.*vkaa);
0072 gk=uk2;
0073 gka=2*sin(alpha)*cos(alpha)*(ux.^2-uy.^2);
0074 gkaa=2*(ux.^2-uy.^2)*(cos(alpha)^2-sin(alpha)^2);
0075 Hpp=2*n/u2;
0076 Halphap=-2*sum((vka.*gk-gka.*vk)./gk.^2);
0077 Halphaalpha=...
0078     sum(fkaa./gk-2*fka.*gka./gk.^2+2*gka.^2.*fk./gk.^3-gkaa.*fk./gk.^2);
0079 NN=2/(Hpp*Halphaalpha-Halphap^2);
0080 var_p=NN*Halphaalpha;
0081 var_alpha=NN*Hpp;
0082 cov_alphap=-NN*Halphap;
0083 Calphap=[var_alpha,var_p,cov_alphap];
0084 % ------ convert to a & b covariance matrix, following DIN 1319 (4)
0085 var_a=var_alpha/cos(alpha)^4;
0086 var_b=(var_alpha*p*p*sin(alpha)^2+var_p*cos(alpha)^2+...
0087     2*cov_alphap*p*sin(alpha)*cos(alpha))/cos(alpha)^4;
0088 cov_ab=(var_alpha*p*sin(alpha)+cov_alphap*cos(alpha))/cos(alpha)^4;
0089 Cab=[var_a,var_b,cov_ab];
0090 % ------ end of uncertainty calculation -----------------------------------
0091 %
0092 %-------------------- plotting section ------------------------------------
0093 if guck~=0
0094     alphatest=linspace(alpha-pi/2,alpha+pi/2,1000);
0095     for k=1:length(alphatest),chitest(k)=chialpha(alphatest(k));end
0096     plot(alphatest,chitest,alphaopt,chiopt,'+r')
0097     xlabel('alpha')
0098     ylabel('\chi^2')
0099     grid on
0100     disp([mfilename,' :: paused, hit any key to continue'])
0101     pause
0102 end
0103 % ------------------ end of plotting section ------------------------------
0104 return
0105 % ------------------ subfunction ------------------------------------------
0106 function chi=chialpha(alpha)
0107     global ux uy x y
0108     uk2=ux.^2*sin(alpha)^2+uy.^2*cos(alpha)^2;
0109     u2=1./mean(1./uk2);
0110     w=u2./uk2;
0111     xbar=mean(w.*x);
0112     ybar=mean(w.*y);
0113     p=ybar*cos(alpha)-xbar*sin(alpha);
0114     chi=sum((y*cos(alpha)-x*sin(alpha)-p).^2./uk2);
0115 return
0116 % ----------- end of subfunction ------------------------------------------

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005