subroutine york(nptsf, xfit, yfit, xfite, yfite, + WK2, A,SIGMAA,B,SIGMAB) include 'typecom' include 'robcom' double precision xfit(*), yfit(*), xfite(*), yfite(*) double precision wk2(*) integer iter c print*,'give eps' c read *,eps if(nptsf .le. 1)then call xtext('NOT ENOUGH POINTS TO CALCULATE!') write(tstring, *)'Number of points = ', nptsf call xtext(tstring) return end if eps = 1.0d-5 istat = 0 call regrwt(nptsf, xfit, yfit, xfite, yfite, wk2, iter, + b, a, sigmab, sigmaa, xbar, ybar, istat, eps) C Print Results call xtext('Results Of York Straight Line Fit Are:') write(tstring, *)' no. of iterations = ', iter call totext(tstring) if(istat .eq. 1)then call xtext('YORK failed') call xtext('No convergence') else write(tstring, *)' no. of iterations = ', iter write(tstring,*)'GRADIENT= ',B,' +/- ',SIGMAB call totext(tstring) write(tstring,*)'INTERCEPT= ',A,' +/- ',SIGMAA call totext(tstring) WRITE(ifil8,*)'RESULTS OF York STRAIGHT LINE FIT ARE:' WRITE(ifil8,*)'GRADIENT= ',B,' +/- ',SIGMAB WRITE(ifil8,*)'INTERCEPT= ',A,' +/- ',SIGMAA end if end C.IDENTIFICATION C C Program REGRWT.FOR C C.AUTHOR C C F. Murtagh, ST-ECF, Garching. Version 1.0 9-July-1986 C (Following discussions with J. Melnick and L. Lucy.) C C.KEYWORDS C C Least squares fitting, linear regression. C C.INPUT/OUTPUT C C.ALGORITHM C C Error weights in both dep. and indep. vbes. C Following initial guess at slope (not using weights), iterative C improvement of slope is carried out until convergence. C C Reference: D. York, "Least squares fitting of a straight line", C Canadian Journal of Physics, 44, 1079-1086, C 1966. C See also: M. Lybanon, "A better least squares method when both C variables have uncertainties", Americal C Journal of Physics, 52, 22-26, 1984. C and: F.S. Acton, "Analysis of Straight Line Data", Dover, C New York, 1959, Chapter 5. C C.INPUTS C C N ... dimension of X, Y, WTX, WTY, W. C X, Y ... vectors of independent and dependent variables. C WTX, WTY ... vectors of weights associated with X and Y. C W ... work vector. C C.OUTPUTS C C ITER ... number of iterations taked by the optimisation algorithm. C SL1 ... slope of least squares solution. C XIN ... intercept of least squares solution. C SIGMA ... uncertainty in slope (sigma). C SIGMAI ... uncertainty in intercept (sigma). C XBAR, YBAR ... centroid of points. C ISTAT ... error indicator (= 1 if no convergence; normally 0). C C--------------------------------------------------------------------- subroutine regrwt(n,x,y,wtx,wty,w,iter,sl1,xin, x sigma,sigmai,xbar,ybar,istat, eps) include 'typecom' double precision x(*),y(*),wtx(*),wty(*),w(*) double precision sl0,sl1,xin,sigma,sigmai,xbar,ybar,a1,a2,b1, x b2,c1,sum,phi,rphi,u,v,sumi,eps c data eps /1.0d-10/ C C Initial guess for the slope, SL0 c 40 CALL SOLVE(X,Y,WTX,WTY,N,SL0) call fitlin(x, y, wty, n, 1, xin, sigmai, SL0, sigma, r) C ITER= 0 C C Reference for the following: York, Can. J. Physics, 1966. C In the following, the "least squares cubic" is solved iteratively. C 55 do 1000 i = 1, n w(i) = wtx(i)*wty(i)/(wtx(i)+sl0*sl0*wty(i)) 1000 continue c xbar = 0.0d0 ybar = 0.0d0 sum = 0.0d0 do 1010 i = 1, n xbar = xbar + w(i)*x(i) ybar = ybar + w(i)*y(i) sum = sum + w(i) 1010 continue xbar = xbar/sum ybar = ybar/sum C C In York's notation, A1, B1 and C1 here become, respectively, C alpha, beta and gamma. C a1 = 0.0d0 a2 = 0.0d0 b1 = 0.0d0 b2 = 0.0d0 c1 = 0.0d0 do 1020 i = 1, n u = x(i) - xbar v = y(i) - ybar a1 = a1 + w(i)*w(i)*u*v/wtx(i) a2 = a2 + w(i)*w(i)*u*u/wtx(i) b1 = b1 + w(i)*w(i)*v*v/wtx(i) b2 = b2 + w(i)*u*u c1 = c1 + w(i)*u*v 1020 continue a1 = 2.0d0*a1/(3.0d0*a2) b1 = (b1 - b2)/(3.0d0*a2) c1 = - c1/a2 c if ((a1*a1 - b1).le.eps) then istat = 2 return endif phi = (a1*a1*a1 - 1.5d0*a1*b1 + 0.5d0*c1)/ x ((a1*a1 - b1)**1.5d0) c rphi is angle given by arccos of the above expr. (in radians). rphi = acos(phi) rphi = rphi/3.0d0 + 4.0d0*3.1415926d0/3.0d0 c sl1 = a1 + 2.0d0*sqrt(a1*a1 - b1)*cos(rphi) c iter = iter+1 if(abs(sl1-sl0) .gt. eps) then sl0 = sl1 if (iter .gt. 100) then istat = 1 return endif goto 55 endif C C Intercept. C XIN = YBAR - SL1*XBAR C C Errors. C sigma = 0.0d0 sigmai = 0.0d0 sum = 0.0d0 sumi = 0.0d0 do 1030 i = 1, n u = x(i) - xbar v = y(i) - ybar sigma = sigma + w(i)*(sl1*u-v)*(sl1*u-v) sum = sum + w(i)*u*u sigmai = sigmai + w(i)*x(i)*x(i) sumi = sumi + w(i) 1030 continue sigma = sigma/(dble(n-2)*sum) sigmai = sigma*sigmai/sumi sigma = sqrt(sigma) sigmai = sqrt(sigmai) c c return end C---------------------------------------------------------------- SUBROUTINE SOLVE(X,Y,ERRX,ERRY,N,SL1) include 'typecom' C C Get estimate of slope, by ignoring weights; C Equation (10), p. 135, of Acton, F.S., "Analysis of Straight- C Line Data", Dover, 1959, is used. C double precision X(*),Y(*),ERRX(*),ERRY(*) double precision SIGMA(5),AVE(4) c REAL*8 W(30),RMS C EQUIVALENCE (SIGMA(1),SXX),(SIGMA(2),SYY),(SIGMA(3),SXY), X (SIGMA(4),SU) ,(SIGMA(5),SV) C C CALL W_MEAN(X,Y,N,ERRX,ERRY,AVE,SIGMA) RHO = 0.0d0 R = RHO*SQRT(SV/SU) R2 = SV/SU A = SXY - R*SXX B = SYY - R2*SXX C = -R2*SXY+R*SYY c slope SL1 = 0.5d0*(B+SQRT(B*B-4.0d0*A*C))/A c Intercept C This was commented out in the original!! C XIN = AVE(4) - SL1*AVE(3) C RETURN END C----------------------------------------------------------------- SUBROUTINE W_MEAN(X,Y,N,ERRX,ERRY,AVE,SIGMA) include 'typecom' double precision X(*),Y(*),ERRX(*),ERRY(*) double precision AVE(*),SIGMA(*) IF(N.LT.2.) RETURN C C AVEX = 0.0d0 AVEY = 0.0d0 WAVX = 0.0d0 WAVY = 0.0d0 WX = 0.0d0 WY = 0.0d0 SU = 0.0d0 SV = 0.0d0 DO I=1,N WTX = 1.0d0/ERRX(I)**2 WTY = 1.0d0/ERRY(I)**2 WX = WX + WTX WY = WY + WTY WAVX = WAVX + X(I)*WTX WAVY = WAVY + Y(I)*WTY AVEX = AVEX + X(I) AVEY = AVEY + Y(I) SU = SU + ERRX(I)*ERRX(I) SV = SV + ERRY(I)*ERRY(I) ENDDO AVEX = AVEX/N AVEY = AVEY/N WAVX = WAVX/WX WAVY = WAVY/WY C SXX = 0.0d0 SYY = 0.0d0 SXY = 0.0d0 DO I=1,N SXX = SXX + (X(I)-AVEX)*(X(I)-AVEX) SYY = SYY + (Y(I)-AVEY)*(Y(I)-AVEY) SXY = SXY + (X(I)-AVEX)*(Y(I)-AVEY) ENDDO C C AVE(1) = AVEX AVE(2) = AVEY AVE(3) = WAVX AVE(4) = WAVY SIGMA(1) = SXX SIGMA(2) = SYY SIGMA(3) = SXY SIGMA(4) = SU/(N-1.0d0) SIGMA(5) = SV/(N-1.0d0) RETURN END