* -------------------------------------------------------------------- * PROGRAM MOWINREG * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Performs a robust regression (least absolute deviation method) or * * a least sum-of-squared deviation regression in a moving window ov- * * er a 2-dim. surface. The least-square regression is fully implemen- * * ted although normally no s.d. for single data points is submitted. * * * * written: (08/20/97) by: Nick Zimmermann * * Utah State University, Dept of Forest Resources * * * * modified: (11/17/97) by: Nick Zimmermann * * * * Version: 1.2 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 2D-Arrays (real) REAL xcov(600,600) ! X-coverage data of actual window REAL ycov(600,600) ! Y-coverage data of actual window REAL icept(600,600) ! Icept-output of actual regression REAL lpsrt(600,600) ! Lsprt-output of actual regression * 1D-Arrays (real) REAL x(90000) ! 1D-array with x-data REAL y(90000) ! 1D-array with y-data REAL sig(90000) ! * Real variables REAL abdev ! REAL a ! Returned Icept of local window position REAL b ! Returned Lpsrt of local window position REAL siga ! REAL sigb ! REAL chi2 ! REAL q ! REAL*4 missin ! Real expression of missing value, as indicated in the header * Integer variables INTEGER ncell(600,600) ! test the # of cells per regression INTEGER mwt ! INTEGER winr ! Row-Size of moving window INTEGER winc ! Col-Size of moving window INTEGER ndata ! Length of 1D-array containing local X/Y-data INTEGER*2 ierr1 ! Input-error (missing file #1) INTEGER*2 ierr2 ! Input-error (missing file #2) INTEGER*2 ncol ! Number of Cols INTEGER*2 nrow ! Number of Rows INTEGER*2 sttc ! Startcell-Pos.(col) for reading actual window into 1D-array INTEGER*2 sttr ! Startcell-Pos.(row) for reading actual window into 1D-array INTEGER*2 endc ! Endcell-Pos. in Cols INTEGER*2 endr ! Endcell-Pos. in Rows INTEGER*2 cend ! End of NCOL-word INTEGER*2 cf ! Right-Pos. in Cols INTEGER*2 cl ! Left-Pos. in Cols INTEGER*2 rf ! Right-Pos. in Rows INTEGER*2 rl ! Left-Pos. in Rows INTEGER*2 r ! Regression method INTEGER*2 p ! % performed INTEGER*2 count ! Count of arguments INTEGER*2 winmin ! Minimum # of cells per regression-window * Character variables CHARACTER*1 equal ! Regressions based on equal # of points CHARACTER*20 miss ! String containing missing value expression CHARACTER*40 cline ! Header-Line #1: NCOL-data CHARACTER*40 rline ! Header-Line #2: NROW-data CHARACTER*40 fline1 ! creates a format to write the icept-output CHARACTER*40 fline2 ! creates a format to write the lpsrt-output CHARACTER*40 fline3 ! creates a format to write the ncell-output CHARACTER*80 xllc ! Header-Line #3: LL-X CHARACTER*80 yllc ! Header-Line #4: LL-Y CHARACTER*80 cell ! Header-Line #5: cellsize CHARACTER*80 nodata ! Header-Line #6: NODATA CHARACTER*32 argmnt(7) ! Splitted arguments passed from the terminal (see GETARGS) * ----------- Read arguments and check for completness --------------- * * OPEN (UNIT=13,FILE='icept.out',STATUS='unknown') OPEN (UNIT=14,FILE='lpsrt.out',STATUS='unknown') OPEN (UNIT=15,FILE='ncell.out',STATUS='unknown') * mwt = 0 count=0 * * DO 1 i=1,7 * CALL GETARG(i,argmnt(i)) * count = count+1 * 1 CONTINUE * CALL GETARG(1,argmnt(1)) IF (argmnt(1)(1:1) .eq. ' ') THEN PRINT 992 STOP ELSE IF (argmnt(1)(1:4) .eq. 'help') THEN PRINT 993 STOP ENDIF * CALL GETARG(2,argmnt(2)) IF (argmnt(2)(1:1) .eq. ' ') THEN PRINT 994 STOP ENDIF * CALL GETARG(3,argmnt(3)) IF (argmnt(3)(1:1) .eq. ' ') THEN PRINT 994 STOP ENDIF READ (argmnt(3),*) winr * DO 20 i=32,1,-1 IF (argmnt(1)(i:i) .NE. ' ') THEN OPEN (UNIT=11,FILE=argmnt(1)(1:i),STATUS='old',IOSTAT=ierr1) IF (ierr1 .NE. 0) THEN PRINT 991, argmnt(1) STOP ENDIF GOTO 30 ENDIF 20 CONTINUE 30 CONTINUE DO 40 i=32,1,-1 IF (argmnt(2)(i:i) .NE. ' ') THEN OPEN (UNIT=12,FILE=argmnt(2)(1:i),STATUS='old',IOSTAT=ierr2) IF (ierr2 .NE. 0) THEN PRINT 991, argmnt(2) STOP ENDIF GOTO 50 ENDIF 40 CONTINUE 50 CONTINUE * CALL GETARG(4,argmnt(4)) IF (argmnt(4)(1:1) .eq. ' ') THEN argmnt(4)(1:32) = argmnt(3)(1:32) ELSE IF (argmnt(4)(1:1) .eq. '#') THEN argmnt(4)(1:32) = argmnt(3)(1:32) ENDIF READ (argmnt(4),*) winc * CALL GETARG(5,argmnt(5)) IF (argmnt(5)(1:1) .eq. ' ') THEN winmin = INT(0.8*winr*winc) if (winmin .lt. 3) winmin=3 ELSE IF (argmnt(5)(1:1) .eq. '#') THEN winmin = INT(0.8*winr*winc) if (winmin .lt. 3) winmin=3 ENDIF IF (winmin .eq. 0) THEN READ (argmnt(5),*) winmin ENDIF * CALL GETARG(6,argmnt(6)) IF (argmnt(6)(1:1) .eq. ' ') THEN r=2 ELSE IF (argmnt(6)(1:1) .eq. '#') THEN r=2 ENDIF IF (r .eq. 0) THEN READ (argmnt(6),*) r ENDIF * CALL GETARG(7,argmnt(7)) IF (argmnt(7)(1:1) .EQ. 'y') argmnt(7)(1:1) = 'Y' IF (argmnt(7)(1:1) .EQ. 'n') argmnt(7)(1:1) = 'N' IF (argmnt(7)(1:1) .eq. ' ') THEN argmnt(7)(1:1) = 'Y' ELSE IF (argmnt(7)(1:1) .eq. '#') THEN argmnt(7)(1:1) = 'Y' ENDIF READ (argmnt(7),*) equal * ----------- Read header and check for errors ---------------------- * * READ (11,995) cline ! check x-file READ (cline,996) ncol READ (11,995) rline READ (rline,996) nrow READ (11,995) xllc READ (11,995) yllc READ (11,995) cell READ (11,995) nodata * nc1 = ncol ! check file-type nr1 = nrow IF (nodata(1:6) .NE. 'NODATA' .OR. cline(1:6).NE.'ncols ') THEN PRINT 51 51 FORMAT(/'One of the INPUT-files is not in GRIDASCII', + '-format!'///) STOP ENDIF * READ (12,995) cline ! check y-file READ (cline,996) ncol READ (12,995) rline READ (rline,996) nrow READ (12,995) xllc READ (12,995) yllc READ (12,995) cell READ (12,995) nodata * nc2 = ncol ! check file-type nr2 = nrow IF (nodata(1:6) .NE. 'NODATA' .OR. cline(1:6).NE.'ncols ') THEN PRINT 51 STOP ENDIF * IF (nc1.NE.nc2 .OR. nr1.NE.nr2) THEN ! check NROW/NCOL PRINT*,'' PRINT*,' The two headers do not match: NROW/NCOL differ!' PRINT*,'' PRINT*,'' STOP ENDIF * IF (winr.GT.nrow .OR. winc.GT.ncol) THEN ! check window-size PRINT*,'' PRINT*,' The selected window is too big for this lattice:' PRINT*,' Window cannot bee bigger than either nrows/ncols!' PRINT*,'' PRINT*,'' STOP ENDIF * IF (r.NE.1 .AND. r.NE.2) THEN PRINT*,'' PRINT*,' Invalid number to select a regression model:' PRINT*,' -type 1 to minimize absolute deviation ' PRINT*,' -type 2 to minimize sum of square of deviations' PRINT*,'' PRINT*,'' STOP ENDIF * IF (equal(1:1).NE.'Y' .AND. equal(1:1).NE.'N') THEN PRINT*,'' PRINT*,' Invalid selection of window-mode:' PRINT*,' -type Y to force equal window size' PRINT*,' -type N to allow smaller peripheral windows' PRINT*,' ' PRINT*,'' PRINT*,'' STOP ENDIF * IF ((winc*winr).LT.winmin .OR. (winc*winr).LT.3) THEN PRINT*,'' PRINT*,' The selected window is too small:' PRINT*,' The window should contain 3 cells at least,' PRINT*,' and must contain more cells than you decla-' PRINT*,' red as the minimum (option: winmin)' PRINT*,'' PRINT*,'' STOP ENDIF * IF (winmin .LT. 3) THEN PRINT*,'' PRINT*,' The selected minimum # of cells/window is too small:' PRINT*,' The window should contain 3 cells at least,' PRINT*,' and must contain more cells than you decla-' PRINT*,' red as the minimum (option: winmin)' PRINT*,'' PRINT*,'' STOP ENDIF * ----------- Read NODATA-expression --------------------------------- * * READ (nodata,990) miss DO 60 i=20,1,-1 IF (miss(i:i) .NE. ' ') THEN cend=i GOTO 61 ENDIF 60 CONTINUE 61 CONTINUE READ (miss,*) missin NODATA(15:25) = '-9999 ' * ----------- Read data into 2D-arrays ------------------------------- * * DO 70 i=1,nrow READ (11,*) (xcov(i,j), j=1,ncol) 70 CONTINUE * DO 80 i=1,nrow READ (12,*) (ycov(i,j), j=1,ncol) 80 CONTINUE * ----------- Initialize the to Output-arrays ------------------------ * * DO 100, i=1,nrow DO 90, j=1,ncol icept(i,j) = -9999. lpsrt(i,j) = -9999. 90 CONTINUE 100 CONTINUE * ----------- Determine Start- and Endcell --------------------------- * * IF (Equal(1:1) .EQ. 'Y') THEN !** If "equally-sized" windows required sttr = winr/2 + 1 ! Start-row sttc = winc/2 + 1 ! Start-col IF (sttr.EQ.0) sttr=1 IF (sttc.EQ.0) sttc=1 * IF (winr .EQ. 1) THEN ! Row=1 endr=nrow ELSEIF (INT(winr/2) - REAL(winr)/2. .EQ. 0.) THEN ! Even number endr = nrow - (winr/2 - 1) ELSE ! Odd number endr = nrow - winr/2 ENDIF * IF (winc .EQ. 1) THEN ! Col=1 endc=ncol ELSEIF (INT(winc/2) - REAL(winc)/2. .EQ. 0.) THEN ! Even numbers endc = ncol - (winc/2 - 1) ELSE ! Odd numbers endc = ncol - winc/2 ENDIF * ELSE !** If "unequally-sized" windows allowed sttr=1 sttc=1 endr=nrow endc=ncol ENDIF * ----------- Read the window-values in a 1D-array ------------------- * * p = 0.0 DO 140, k=sttr,endr DO 130, l=sttc,endc * rf = k-winr/2 cf = l-winc/2 IF (rf.LT.1) rf=1 IF (cf.LT.1) cf=1 * IF (INT(winr/2) - REAL(winr)/2. .EQ. 0.) THEN ! Even numbers rl = k+(winr/2-1) ! End-col IF (rl.GT.nrow) rl=nrow ELSE ! Odd numbers rl = k+winr/2 ! End-col IF (rl.GT.nrow) rl=nrow ENDIF * IF (INT(winc/2) - REAL(winc)/2. .EQ. 0.) THEN ! Even numbers cl = l+(winc/2-1) ! End-row IF (cl.GT.ncol) cl=ncol ELSE ! Odd numbers cl = l+winc/2 ! End-row IF (cl.GT.ncol) cl=ncol ENDIF * DO 104, i=1,ndata ! Use old ndata- x(i) = 0. ! value to clean y(i) = 0. ! regression-queue 104 CONTINUE ! .... and ..... ndata = 0 ! reset "ndata" * DO 120, i=rf,rl DO 110, j=cf,cl IF (xcov(i,j).NE.missin.AND.ycov(i,j).ne.missin) THEN ndata = ndata+1 x(ndata) = xcov(i,j) y(ndata) = ycov(i,j) ENDIF 110 CONTINUE 120 CONTINUE * ------ Perform local regression IF (r.EQ.1 .AND. ndata.GE.winmin) THEN CALL MEDFIT(x,y,ndata,a,b,abdev) ELSEIF (ndata .GE. winmin) THEN CALL FIT(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q) ELSE GOTO 130 ENDIF icept(k,l) = a lpsrt(k,l) = b ncell(k,l) = ndata 130 CONTINUE p = INT(100 * (REAL(k)/REAL(endr))) ! Calculate progress count=p-(Mod(p,10)) IF (p-count .eq. 0.0) THEN PRINT 998, (p-Mod(p,10)) ! Show progress ENDIF 140 CONTINUE * ----------- Write results to 2D-array (OUTPUT) --------------------- * * DO 160 i=17,15,-1 IF (cline(i:i) .NE. ' ') THEN cend=i GOTO 170 ENDIF 160 CONTINUE 170 CONTINUE fline1 = '('//cline(15:cend)//'(F13.6,1X))' fline2 = '('//cline(15:cend)//'(F16.9,1X))' fline3 = '('//cline(15:cend)//'(I3,1X))' * WRITE (13,999) cline,rline,xllc,yllc,cell,nodata WRITE (14,999) cline,rline,xllc,yllc,cell,nodata WRITE (15,999) cline,rline,xllc,yllc,cell,nodata DO 180, i=1,nrow WRITE (13,fline1) (icept(I,J), j=1,ncol) WRITE (14,fline2) (lpsrt(i,j), j=1,ncol) WRITE (15,fline3) (ncell(i,j), j=1,ncol) 180 CONTINUE * ----------- FORMAT-statements in the MAIN-program ------------------ * * 990 FORMAT(T15,A20) 991 FORMAT(/1X,'Cannot find: ',A32/1X,'Check for proper file-name!'//) 992 FORMAT(//1X,' MoWinReg: Regressions in a Moving Window:'/1X, + ' ---------'/1X, + ' This program reads two 2D-arrays in GRIDASCII-', + 'format and calculates multiple'/1X,' linear regressions ', + 'in a moving window of variable size (max=300x300 pixels)', + '.'/' You can choose from two regression methods, addres', + 's the minimum number of '/, + ' pixels for the regression, and set the procedure to ful', + 'l or partial window '/, + ' mode. The output (Lpsrt & Icept) is written into two fi', + 'les in the same for-'/, + ' mat. Import the files into ArcInfo using the ASCIIG', + 'RID-command. The result'/, + ' of each regression is written to the center cell of the', + ' actual window. Using'/, + ' the partial window mode, up to 75% of the actual window', + ' is empty (corners;'/, + ' 50% along the sidelines). You can force the program to ', + 'calculate regressions'/, + ' in full windows only. Alternatively you can set the min', + 'mum number of cells,'/, + ' necessary to run a regresion. If the number of cells is', + ' too small, a missing'/, + ' value is written to the windows'' center cell.', + //1X, + ' mowinreg help for more information'///, + T36,'written by:'/, + T41,'Niklaus E. Zimmermann,'/1X, + T41,'USU, Logan, UT 84322-5215'/1X, + T41,''/1X, + T36,'now: Swiss Federal Research Institute WSL'/1X, + T41,'CH-8903 Birmensdorf, Switzerland'/1X, + T41,'(niklaus.zimmermann@wsl.ch)'//) 993 FORMAT(//'Usage: mowinreg []', + ' [] [] []'//, + ' Xcover: a GRIDASCII-file containing X-values'/, + ' Ycover: a GRIDASCII-file containing Y-values'/, + ' Wnrows: Number of (horizontal) rows of the moving ', + 'window'/, + ' Wncols: Number of (vertical) cols of the moving ', + 'window'/, + ' maximum size for wncol*wnrow is 90,000'/, + ' DEFAULT: wncols = wnrows'/, + ' Winmin: Minimum number of pixels per window used to ', + 'calculate a regression.'/, + ' If the number in a window is lower (Peripher', + 'y/Nodata), a nodata-va-'/, + ' lue is written to the center cell, and the', + 'window moves to the next'/, + ' position. ', + 'DEFAULT: winmin = (0.8*(wncols*winrows))'/, + ' Method: Select a method to fit the linear regression', + ' model:'/, + ' 1 = minimize absolute deviations (robust)'/, + ' 2 = minimize sum of squares of deviations'/, + ' DEFAULT: method = 2 (least-squares regress', + 'ion)'/, + ' Equal: Type Y to allow only equal-sized windows'/, + ' Type N to allow smaller sized peripheral ', + 'windows',/, + ' With you force MoWinReg to skip calcu', + 'lations at peripheral cells'/, + ' (=1/2 window), since MoWinReg writes the ', + 'output of the regression'/, + ' into the center cell of the actual window.'/, + ' DEFAULT: equal = Y'//, + ' Defaults can be selected by typing ''#'' instead ', + 'of a value. If all optional argu-'/, + ' ments shall be default, type the required (3)', + ' arguments in the command line only.'//, + ' For a general overview, type: movinreg '///) 994 FORMAT(//,' The argument-list is incomplete, type:'/, + ' mowinreg help , for more information'///) 995 FORMAT(A) 996 FORMAT(T15,I6) 997 FORMAT(T15,F7.1) 998 FORMAT(' Running... ',I3,'% done') 999 FORMAT(2(A40/),3(A80/),A80) * END * * ==================================================================== * SUBROUTINE MEDFIT(X,Y,NDATA,A,B,ABDEV) * Fits y=a+bx by the criterion of least absolute deviations. The ar- * * rays X and Y, of length NDATA, are the input experimental points. * * The fitted parameters A & B are output, along with ABDEV which is * * the mean absolute deviation (in Y) of the experimental points from * * the fitted line. This routine uses the routine ROFUNC, with commu- * * nication via a common block. [~ltsreg in S-Plus] * * Press, W.H., et al., 1989. Numerical Recipes: p544. * PARAMETER (NMAX=1000) EXTERNAL ROFUNC COMMON /ARRAYS/ NDATAT,XT(NMAX),YT(NMAX),ARR(NMAX),AA,ABDEVT * DIMENSION X(NDATA),Y(NDATA) SX = 0. SY = 0. SXY= 0. SXX= 0. * DO 11, J=1,NDATA ! As a first guess for A and B, we XT(J) = X(J) ! will find the least-squares fit- YT(J) = Y(J) ! ting line. SX = SX + X(J) SY = SY + Y(J) SXY = SXY + X(J)*Y(J) SXX = SXX + X(J)**2 11 CONTINUE * NDATAT = NDATA DEL = (NDATA*SXX)-(SX**2) AA = ((SXX*SY)-(SX*SXY)) / DEL ! Least-squares solution BB = ((NDATA*SXY)-(SX*SY)) / DEL CHISQ = 0. * DO 12, J=1,NDATA CHISQ = CHISQ + (Y(J)-(AA+BB*Y(J)))**2 12 CONTINUE * SIGB = SQRT(CHISQ/DEL) ! The sd will give some idea of how B1 = BB ! big an iteration step to take. F1 = ROFUNC(B1) B2 = BB+SIGN(3.*SIGB,F1) ! Guess bracket as 3 sd away, in the F2 = ROFUNC(B2) ! downhill direction know from F1. * 1 IF (F1*F2 .GT. 0.) THEN ! Bracketing. BB = 2. * B2-B1 B1 = B2 F1 = F2 B2 = BB F2 = ROFUNC(B2) GO TO 1 ENDIF * SIGB = 0.01*SIGB ! Refine until error a negl. nb. of sd * 2 IF (ABS(B2-B1) .GT. SIGB) THEN ! Bisection. BB = 0.5 * (B1+B2) IF (BB .EQ. B1 .OR. BB .EQ. B2) GO TO 3 F = ROFUNC(BB) IF (F*F1 .GE. 0.) THEN F1 = F B1 = BB ELSE F2 = F B2 = BB ENDIF GO TO 2 ENDIF * 3 A = AA B = BB ABDEV = ABDEVT/NDATA * RETURN END * ==================================================================== * FUNCTION ROFUNC(B) * Evaluates the right-hand side of equation 14.6.16 (parameter b-es- * * timation) for a given value B. Communication with the subroutine * * MEDFIT is through a common block statement. * * Press, W.H., et al., 1989. Numerical Recipes: p545. * PARAMETER (NMAX=1000) COMMON /ARRAYS/ NDATA,X(NMAX),Y(NMAX),ARR(NMAX),AA,ABDEV * N1 = NDATA + 1 NML = N1/2 NMH = N1-NML * DO 11, J=1,NDATA ARR(J) = Y(J) - B*X(J) 11 CONTINUE * CALL SORT(NDATA,ARR) * AA = 0.5 * (ARR(NML)+ARR(NMH)) SUM = 0. ABDEV = 0. * DO 12, J=1,NDATA D = Y(J) - (B*X(J) + AA) ABDEV = ABDEV + ABS(D) SUM = SUM + X(J)*SIGN(1.0,D) 12 CONTINUE * ROFUNC = SUM * RETURN END * ==================================================================== * SUBROUTINE SORT(N,ARRAY) * Sorts an array ARRAY of length N into ascending numerical order * * using the Heapsort algorithm. N is input; ARRAY is replaced on * * output by its sorted rearrangement. * * Press, W.H., et al., 1989. Numerical Recipes: p231. * DIMENSION ARRAY(N) L = N/2 + 1 IR = N * * The index L will be decremented from its initial value down to * * 1 during the "hiring" (heap creation) phase. Once it reaches 1, * * the index IR will be decremented from its initial value down to * * 1 during the "retirement-and-promotion" (heap selection) phase. * * 10 CONTINUE IF (L .GT. 1) THEN ! Still in hiring phase. L = L-1 RRA = ARRAY(l) ELSE ! In retirement-&-promotion phase. RRA = ARRAY(IR) ! Clear a space at end of array. ARRAY(IR) = ARRAY(1) ! Retire the top of the heap into it. IR = IR-1 ! Decrease the size of the corpor'n. IF (IR .EQ. 1) THEN ! Done with the last promotion. ARRAY(1) = RRA ! The least component workerof all! * RETURN * ENDIF ENDIF * I = L ! Whether we are in "hiring" or "promotion", we he- J = L+L ! re set up to shift down rra to its proper level. * 20 IF (J .LE. IR) THEN ! Do while L. LE. IR IF (J .LT. IR) THEN IF (ARRAY(J) .LT. ARRAY(J+1)) J=J+1 ENDIF IF (RRA .LT. ARRAY(J)) THEN ! Demote RRA. ARRAY(I) = ARRAY(J) I = J J = J+J ELSE ! This is RRA's level. Set J to J = IR+1 ! terminate the shift-down. ENDIF * GO TO 20 ENDIF * ARRAY(I) = RRA ! Put RRA into its slot. * GO TO 10 * END * ==================================================================== * SUBROUTINE FIT(X,Y,NDATA,SIG,MWT,A,B,SIGA,SIGB,CHI2,Q) * Given a set of NDATA points X(I), Y(I) with standard deviations * * SIG(I), fit them to a straight line Y=A+BX by minimizing CHI**2. * * Returned are A,B and their respective probable uncertainties SIGA * * and SIGB, the chi-square CHI**2, and the goodness-of-fit probabi- * * lity Q (that the fit would have CHI**2 this large or larger). If * * MWT=0 on input, then the standard deviations are assumed to be un- * * available: Q is returned as 1.0 and the normalization of CHI**2 is * * to unit standard deviation on all points. (Press et al, p508) * * Press, W.H., et al., 1989. Numerical Recipes: p508. * DIMENSION X(NDATA),Y(NDATA),SIG(NDATA) SX = 0. ! Initialize sums to zero. SY = 0. ST2 = 0. * IF (MWT .NE. 0) THEN ! Accumulate sums SS = 0. DO 11, I=1,NDATA ! With weights .... WT = 1./(SIG(I)**2) SS = SS + WT SX = SX + X(I)*WT SY = SY + Y(I)*WT 11 CONTINUE ELSE ! ... or without weights. DO 12, I=1,NDATA SX = SX + X(I) SY = SY + Y(I) 12 CONTINUE SS = FLOAT(NDATA) ENDIF * SXOSS = SX/SS * IF (MWT .NE. 0) THEN DO 13,I=1,NDATA T = (X(I)-SXOSS) / SIG(I) ST2 = ST2 + T*T B = B + T*Y(I)/SIG(I) 13 CONTINUE ELSE DO 14,I=1,NDATA T = X(I) - SXOSS ST2 = ST2 + T*T B = B + T*Y(I) 14 CONTINUE ENDIF * B = B/ST2 ! Solve for A, B, sd(A) and sd(B) A = (SY - SX*B) / SS SIGA = SQRT((1. + SX*SY / (SS*ST2)) / SS) SIGB = SQRT(1. / ST2) CHI2 = 0. ! Calculate chi**2 * IF (MWT .EQ. 0) THEN DO 15, I=1,NDATA CHI2 = CHI2 + (Y(I) - A - B*X(I))**2 15 CONTINUE * Q = 1. SIGDAT = SQRT(CHI2 / (NDATA-2)) ! For unweighted data evaluate typi- SIGA = SIGA*SIGDAT ! cal SIG using chi**2, and adjust SIGB = SIGB*SIGDAT ! the standard deviations. ELSE DO 16, I=1,NDATA CHI2 = CHI2 + ((Y(I) - A - B*X(I)) / SIG(I))**2 16 CONTINUE Q = GAMMQ(0.5 * (NDATA-2), 0.5*CHI2) ENDIF * RETURN END * ==================================================================== * FUNCTION GAMMQ(A,X) * Returns the incomplete gamma-function P(a,x) -> 1-P(a,x). * * Press, W.H., et al., 1989. Numerical Recipes: p162. * IF (X .LT. 0. .OR. A .LE. 0.) PAUSE IF (X .LT. A+1.) THEN ! Use the series repersentation CALL GSER(GAMSER,A,X,GLN) GAMMQ = 1. - GAMSER ! and take its complement. ELSE ! Use the continued fraction representation CALL GCF(GAMMCF,A,X,GLN) GAMMQ = GAMMCF ENDIF * RETURN END * ==================================================================== * SUBROUTINE GSER(GAMSER,A,X,GLN) * Returns the incomplete gamma function P(a,x) evaluated by its se- * * ries representation as GAMSER. Also returns ln(..a) as GLN. * * Press, W.H., et al., 1989. Numerical Recipes: p162. * PARAMETER (ITMAX=100, EPS=3.E-7) * GLN = GAMMLN(A) IF (X .LE. 0.) THEN IF (X .LT. 0.) PAUSE GAMSER = 0. * RETURN ENDIF * AP = A SUM = 1./A DEL = SUM DO 11, N=1,ITMAX AP = AP + 1. DEL = DEL*X/AP SUM = SUM+DEL IF (ABS(DEL) .LT. ABS(SUM)*EPS) GO TO 1 11 CONTINUE * PAUSE 'A too large, ITMAX too small' * 1 GAMSER = SUM*EXP(-X + A*LOG(X) - GLN) * RETURN END * ==================================================================== * SUBROUTINE GCF(GAMMCF,A,X,GLN) * Returns the incomplete gamma function Q(a,x) evaluated by its con- * * tinued fraction representation as GAMMCF. Also returns ln(..a) * * as GLN. Press, W.H., et al., 1989. Numerical Recipes: p162. * PARAMETER (ITMAX=100, EPS=3.E-7) * GLN = GAMMLN(A) GOLD = 0. ! This is the previous value, tested A0 = 1. ! against for convergence. A1 = X ! We are setting up A's & B's of equation(5.2.4) B0 = 0. ! for evaluating the continued fraction. B1 = 1. FAC = 1. ! FAC is the renomaliz. factor for prevent. overflow DO 11, N=1,ITMAX ! of the partial numerators and denominators. AN = FLOAT(N) ANA = AN-A A0 = (A1 + A0*ANA)*FAC ! One step of the recurrence (5.2.5) B0 = (B1 + B0*ANA)*FAC ANF = AN*FAC A1 = X*A0 + ANF*A1 ! The next step of the recurr. (5.2.5) B1 = X*B0 + ANF*B1 * IF (A1 .NE. 0.) THEN ! Shall we normalize? FAC = 1./A1 ! Yes. Set FAC so it happens. G = B1*FAC ! New value of answer. IF (ABS((G-GOLD)/G) .LT. EPS) GO TO 1 ! Converged? If so, exit. GOLD = G ! If not. Save value. ENDIF * 11 CONTINUE * PAUSE 'A too large, ITMAX too small' 1 GAMMCF = EXP(-X + A*ALOG(X) - GLN)*G ! Put factors in front. * RETURN END * ==================================================================== * FUNCTION GAMMLN(XX) * See Press, W.H., et al., 1989. Numerical Recipes: p157. for full * * explanation... * REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER * * Internal arithmetic could be done in double precision, a nicety that * has been omitted here; five-figure accuracy is considered as good enough * DATA COF,STP / 76.18009173,-86.50532033,24.01409822, + -1.231739516,.120858003E-2,-.536382E-5,2.50662827465 / DATA HALF,ONE,FPF / 0.5,1.0,5.5 / * X = XX-ONE TMP = X+FPF TMP = (X+HALF) * LOG(TMP) - TMP SER = ONE * DO 11, J=1,6 X = X+ONE SER = SER + COF(J)/X 11 CONTINUE * GAMMLN = TMP+LOG(STP*SER) * RETURN END * ==================================================================== *