C	From edf@astro.psu.edu Tue Jan 16 01:57 EST 1996
C	From: "Eric D. Feigelson" <edf@astro.psu.edu>
C	To: ali@merlin.anu.edu.au
C	Subject: SLOPES.f
C
C     *                    PROGRAM SLOPES                          *
C     *                                                            *
C     *  PROGRAM TO COMPUTE THE THEORETICAL REGRESSION SLOPES      *
C     *  AND UNCERTAINTIES (VIA DELTA METHOD), AND UNCERTAINTIES   *
C     *  VIA BOOTSTRAP AND BIVARIATE NORMAL SIMULATION FOR A       *
C     *  (X(I),Y(I)) DATA SET WITH UNKNOWN POPULATION DISTRIBUTION.*
C     *                                                            *
C     *     WRITTEN BY ERIC FEIGELSON, PENN STATE.  JUNE 1991      *
C     *                                                            *
C     *    INPUT (FROM A FILE : FILE NAME FORMAT=<A14>              *
C     *         THE FILE CONTAINS X(I), INDEPENDENT VARIABLE, AND  *
C     *         Y(I) WITH I UP TO 1000.  FORMAT 2F10.3             *
C     *                                                            *
C     *    OUTPUT -- FILE `SLOPES.OUT' (ADJUSTABLE BY USER)        *
C     *                                                            *
C     *    SUBROUTINES --                                          *
C     *         DATSTT(NMAX,NTOT,X,VARX,Y,VARY,VARXY,RHO),         *
C     *              CALCULATES SIMPLE STATISTICS OF DATASET       *
C     *         SIXLIN(NMAX,NTOT,X,Y,A,SIGA,B,SIGB), CALCULATES 6  *
C     *              LINEAR FITS ANALYTICALLY                      *
C     *         BOOTSP(N,X,Y,XSIM,YSIM,ISEED), MAKES BOOTSTRAP     *
C     *              SIMULATED DATASETS. CALLS FUNCTION RAN3.      *
C     *         TULSIM (VERS. 3.2) BIVARIATE NORMAL SIMULATIONS    *
C     *         (M. T. BOSWELL, C. 1989, MS-DOS ONLY):  
C     *                BASGET                                      *
C     *                RBNORM(X4,Y4,N,RXY(IT),RMU,SD)              *
C     *                BASSAV                                      *
C     *                                                            *
C     A full description of these methods can be found in:
C       Isobe, T., Feigelson, E. D., Akritas, M. and Babu, G. J.,
C          Linear Regression in Astronomy I, Astrophys. J. 364, 104 
C          (1990)
C       Babu, G. J. and Feigelson, E. D., Analytical and Monte Carlo
C          Comparisons of Six Different Linear Least Squares Fits,
C          Communications in Statistics, Simulation & Computation, 
C          21, 533 (1992)
C       Feigelson, E. D. and Babu, G. J., Linear Regression in 
C          Astronomy II, Astrophys. J. 397, 55 (1992).
C
C
C     *  THE USER MAY CHANGE PARAMETER NMAX AS NEEDED              *
C

       subroutine slopy(OrthA,OrthB,OrthC,OrthD) 

       PARAMETER (NMAX = 203)     ! XXX modify to the no. of data points
       IMPLICIT REAL*8 (A-H,O-Z)
       CHARACTER*14 FILE
       DIMENSION X(NMAX),Y(NMAX),XSIM(NMAX),YSIM(NMAX)
       DIMENSION A(6),SIGA(6),B(6),SIGB(6)
       DIMENSION ASIM(6),SASIM(6),BSIM(6),SBSIM(6),
     +           ASUM(6),ASSUM(6),BSUM(6),BSSUM(6),
     +           AAVG(6),SDA(6),BAVG(6),SDB(6)
C        COMMENT OUT THE FOLLOWING LINE IF TULSIM IS UNAVAILABLE
C       REAL*4 AVG4(2),SD4(2),RHO4,X4(NMAX),Y4(NMAX)

       FILE='TMP.DAT'   ! Added by Alister
       OPEN (UNIT=50, FILE='TMP.DAT', STATUS='OLD', FORM='FORMATTED')
C
C     *            READ THE DATA SET AND OPEN OUTPUT FILE
C
       read (50,*) nummax 
   15  do i=1,nummax
         X(I) = 0.0
         Y(I) = 0.0
         READ(50,14) X(I),Y(I)
       enddo
   25  CONTINUE
       NTOT=nummax
       CLOSE (UNIT = 50)
       OPEN (UNIT=10, FILE='TMP.OUT', STATUS='UNKNOWN')
       WRITE(10,28)
       WRITE(10,16)
       WRITE(10,17)
       WRITE(10,28)
       WRITE(10,19) FILE,NTOT
       WRITE(10,28)
C
C     * OPTIONAL OUTPUT OF DATA
C
       WRITE(10,999) (X(J),Y(J),J=1,NTOT)
C
C     *         COMPUTE MEANS, STANDARD DEVIATIONS AND CORREL. COEFF.
C
       CALL DATSTT(NumMAX,NTOT,X,Y,XAVG,VARX,YAVG,VARY,VARXY,RHO)
       SIGX = DSQRT(VARX)
       SIGY = DSQRT(VARY)
       RN = DFLOAT(NTOT)
       SUMXX = VARX*RN
       SUMYY = VARY*RN
       SUMXY = VARXY*RN
       WRITE(10,28)
       WRITE(10,21) XAVG, SIGX, YAVG, SIGY, SUMXX, SUMYY, SUMXY, RHO
C
C     *    COMPUTE THE SIX LINEAR REGRESSION WITH ANALYTIC ERROR ANALYSIS
C
       CALL SIXLIN(NumMAX,NTOT,X,Y,A,SIGA,B,SIGB)
       WRITE(10,28)
       WRITE(10,30)
       WRITE(10,32)
       WRITE(10,28)
       WRITE(10,34)
       WRITE(10,28)
       WRITE(10,36) A(1),SIGA(1),B(1),SIGB(1)
       WRITE(10,38) A(2),SIGA(2),B(2),SIGB(2)
       WRITE(10,40) A(3),SIGA(3),B(3),SIGB(3)
       WRITE(10,42) A(4),SIGA(4),B(4),SIGB(4)
       WRITE(10,44) A(5),SIGA(5),B(5),SIGB(5)
       WRITE(10,46) A(6),SIGA(6),B(6),SIGB(6)

c       OrthA=A(5)   !XXX
c       OrthB=SIGA(5)    !XXX
c       OrthC=B(5)   !XXX
c       OrthD=SIGB(5)    !XXX
c       Use solutions from Jackknife sampling below

C
C     *     MAKE BOOTSTRAP SIMULATED DATASETS, AND COMPUTE AVERAGES
C     *     AND STANDARD DEVIATIONS OF REGRESSION COEFFICIENTS
C
c       WRITE(*,101)
c       READ(*,102) NSIM
       NSIM=4000
c       WRITE(*,103)
c       READ(*,102) ISEED
       ISEED=786
       RNSIM = DFLOAT(NSIM)
       
       DO 100 J = 1,6
          ASUM(J) = 0.0
          ASSUM(J) = 0.0
          BSUM(J) = 0.0
          BSSUM(J) = 0.0
          SDA(J) = 0.0
          SDB(J) = 0.0
  100  CONTINUE
       DO 110 I = 1, NSIM
          IF(I.EQ.1) ISEED = -1*IABS(ISEED)
          CALL BOOTSP(NumMAX,NTOT,X,Y,XSIM,YSIM,ISEED)
          CALL SIXLIN(NumMAX,NTOT,XSIM,YSIM,ASIM,SASIM,BSIM,SBSIM)
C                 OPTIONAL DETAILED PRINTOUT
C                 WRITE(10,99) I,ASIM(1),SASIM(1),BSIM(1),SBSIM(1)
          DO 120 J = 1,6
             ASUM(J) = ASUM(J) + ASIM(J)
             ASSUM(J) = ASSUM(J) + ASIM(J)*ASIM(J)
             BSUM(J) = BSUM(J) + BSIM(J)
             BSSUM(J) = BSSUM(J) + BSIM(J)*BSIM(J)
  120     CONTINUE
  110  CONTINUE
       DO 130 J = 1,6
          AAVG(J) = ASUM(J)/RNSIM
          SDTEST = ASSUM(J) - RNSIM*AAVG(J)*AAVG(J)
          IF(SDTEST.GT.0.0) SDA(J) = DSQRT(SDTEST/(RNSIM-1.0))
          BAVG(J) = BSUM(J)/RNSIM
          SDTEST = BSSUM(J) - RNSIM*BAVG(J)*BAVG(J)
          IF(SDTEST.GT.0.0) SDB(J) = DSQRT(SDTEST/(RNSIM-1.0))
  130  CONTINUE
       WRITE(10,28)
       WRITE(10,28)
       WRITE(10,28)
       WRITE(10,28)
       WRITE(10,105) NSIM,ISEED
       WRITE(10,28)
       WRITE(10,34)
       WRITE(10,28)
       WRITE(10,36) AAVG(1),SDA(1),BAVG(1),SDB(1)
       WRITE(10,38) AAVG(2),SDA(2),BAVG(2),SDB(2)
       WRITE(10,40) AAVG(3),SDA(3),BAVG(3),SDB(3)
       WRITE(10,42) AAVG(4),SDA(4),BAVG(4),SDB(4)
       WRITE(10,44) AAVG(5),SDA(5),BAVG(5),SDB(5)
       WRITE(10,46) AAVG(6),SDA(6),BAVG(6),SDB(6)

C
C     *     MAKE JACKKNIFE SIMULATED DATASETS AND COMPUTE AVERAGES AND
C     *     STANDARD DEVIATIONS OF REGRESSION COEFFICIENTS.  NOTE THAT
C     *     THE S.D. FORMULAE ARE DIFFERENT THAN FOR THE OTHER ESTIMATORS.
C
       RNTOT = DFLOAT(NTOT)
       DO 150 J = 1,6
          ASUM(J) = 0.0
          ASSUM(J) = 0.0
          BSUM(J) = 0.0
          BSSUM(J) = 0.0
          SDA(J) = 0.0
          SDB(J) = 0.0
  150  CONTINUE
       NJACK = NTOT - 1
       DO 160 I = 1, NTOT
          DO 170 K = 1, NTOT
             IF(K.EQ.I) GO TO 170
             IF(K.LT.I) XSIM(K) = X(K)
             IF(K.LT.I) YSIM(K) = Y(K)
             IF(K.GT.I) XSIM(K-1) = X(K)
             IF(K.GT.I) YSIM(K-1) = Y(K)
  170     CONTINUE
          CALL SIXLIN(NumMAX,NJACK,XSIM,YSIM,ASIM,SASIM,BSIM,SBSIM)
          DO 175 J = 1, 6
             ASUM(J) = ASUM(J) + ASIM(J)
             BSUM(J) = BSUM(J) + BSIM(J)
             ASSUM(J) = ASSUM(J) + ASIM(J)*ASIM(J)
             BSSUM(J) = BSSUM(J) + BSIM(J)*BSIM(J)
  175     CONTINUE
C                 OPTIONAL DETAILED PRINTOUT
C                 WRITE(10,99) I,ASIM(1),SASIM(1),BSIM(1),SBSIM(1)
  160  CONTINUE
       DO 190 J = 1,6
          AAVG(J) = ASUM(J)/RNTOT
          SDTEST = ASSUM(J) - RNTOT*AAVG(J)*AAVG(J)
          IF(SDTEST.GT.0.0) SDA(J) = DSQRT((RNTOT-1.0)*SDTEST/RNTOT)
          BAVG(J) = BSUM(J)/RNTOT
          SDTEST = BSSUM(J) - RNTOT*BAVG(J)*BAVG(J)
          IF(SDTEST.GT.0.0) SDB(J) = DSQRT((RNTOT-1.0)*SDTEST/RNTOT)
  190  CONTINUE
       WRITE(10,28)
       WRITE(10,28)
       WRITE(10,28)
       WRITE(10,28)
       WRITE(10,155)
       WRITE(10,28)
       WRITE(10,34)
       WRITE(10,28)
       WRITE(10,36) AAVG(1),SDA(1),BAVG(1),SDB(1)
       WRITE(10,38) AAVG(2),SDA(2),BAVG(2),SDB(2)
       WRITE(10,40) AAVG(3),SDA(3),BAVG(3),SDB(3)
       WRITE(10,42) AAVG(4),SDA(4),BAVG(4),SDB(4)
       WRITE(10,44) AAVG(5),SDA(5),BAVG(5),SDB(5)
       WRITE(10,46) AAVG(6),SDA(6),BAVG(6),SDB(6)

       OrthA=AAVG(5)   !XXX
       OrthB=SDA(5)    !XXX
       OrthC=BAVG(5)   !XXX
       OrthD=SDB(5)    !XXX

C
C     *     MAKE BIVARIATE NORMAL SIMULATED DATASETS AND COMPUTE AVERAGES
C     *     AND STANDARD DEVIATIONS OF REGRESSION COEFFICIENTS
C     *     (COMMENT OUT CODE FROM HERE TO `CALL BASSAV' IF THERE IS NO
C     *     ACCESS TO TULSIM MS-DOS ROUTINES)
C
C       RHO4 = RHO
C       AVG4(1) = XAVG
C       AVG4(2) = YAVG
C       SD4(1) = SIGX
C       SD4(2) = SIGY
C       CALL BASGET
C       RSUM = 0.0
C       RSSUM = 0.0
C       DO 200 J = 1,6
C          ASUM(J) = 0.0
C          ASSUM(J) = 0.0
C          BSUM(J) = 0.0
C          BSSUM(J) = 0.0
C          SDA(J) = 0.0
C          SDB(J) = 0.0
C          RSD = 0.0
C  200  CONTINUE
C       DO 210 I = 1, NSIM
C          CALL RBNORM(X4,Y4,NTOT,RHO4,AVG4,SD4)
C          RSUM = RSUM + DFLOAT(RHO4)
C          RSSUM = RSSUM + DFLOAT(RHO4*RHO4)
C          DO 220 L = 1, NTOT
C             XSIM(L) = DFLOAT(X4(L))
C             YSIM(L) = DFLOAT(Y4(L))
C  220     CONTINUE
C          CALL SIXLIN(NumMAX,NTOT,XSIM,YSIM,ASIM,SASIM,BSIM,SBSIM)
CC                 OPTIONAL DETAILED PRINTOUT
CC                 WRITE(10,99) I,ASIM(1),SASIM(1),BSIM(1),SBSIM(1)
C          DO 230 J = 1,6
C             ASUM(J) = ASUM(J) + ASIM(J)
C             ASSUM(J) = ASSUM(J) + ASIM(J)*ASIM(J)
C             BSUM(J) = BSUM(J) + BSIM(J)
C             BSSUM(J) = BSSUM(J) + BSIM(J)*BSIM(J)
C  230     CONTINUE
C  210  CONTINUE
C       DO 240 J = 1,6
C          AAVG(J) = ASUM(J)/RNSIM
C          SDTEST = ASSUM(J) - RNSIM*AAVG(J)*AAVG(J)
C          IF(SDTEST.GT.0.0) SDA(J) = DSQRT(SDTEST/(RNSIM-1.0))
C          BAVG(J) = BSUM(J)/RNSIM
C          SDTEST = BSSUM(J) - RNSIM*BAVG(J)*BAVG(J)
C          IF(SDTEST.GT.0.0) SDB(J) = DSQRT(SDTEST/(RNSIM-1.0))
C  240  CONTINUE
C       WRITE(10,28)
C       WRITE(10,28)
C       WRITE(10,28)
C       WRITE(10,28)
C       WRITE(10,205) NSIM
C       WRITE(10,28)
C       WRITE(10,34)
C       WRITE(10,28)
C       WRITE(10,36) AAVG(1),SDA(1),BAVG(1),SDB(1)
C       WRITE(10,38) AAVG(2),SDA(2),BAVG(2),SDB(2)
C       WRITE(10,40) AAVG(3),SDA(3),BAVG(3),SDB(3)
C       WRITE(10,42) AAVG(4),SDA(4),BAVG(4),SDB(4)
C       WRITE(10,44) AAVG(5),SDA(5),BAVG(5),SDB(5)
C       WRITE(10,46) AAVG(6),SDA(6),BAVG(6),SDB(6)
C       WRITE(10,28)
C       RSUM = RSUM/RNSIM
C       SDTEST = RSSUM - RNSIM*RSUM*RSUM
C       IF(SDTEST.GT.0.0) RSD = DSQRT(SDTEST/(RNSIM-1.0))
C       WRITE(10,207) RSUM,RSD
C       CALL BASSAV
C
C     *     INPUT FORMATS
C
   10  FORMAT(' ENTER THE DATA FILE : ',$)
  101  FORMAT('  NUMBER OF SIMULATIONS: ')
  102  FORMAT(BN,I6)
  103  FORMAT('  SEED (INTEGER, I6): ')
   12  FORMAT(A14)
   14  FORMAT(2F10.3)
C
C     *     OUTPUT FORMATS
C
   16  FORMAT(6X,' SLOPES:  ANALYTICAL AND SIMULATION CALCULATIONS OF')
   17  FORMAT(12X,'  LINEAR REGRESSIONS AND UNCERTAINTIES')
   28  FORMAT(' ')
   19  FORMAT(10X, ' INPUT DATA FILE:  ',A14,'  # POINTS:  ',I4)
   21  FORMAT(5X,'XAVG = ',F10.5,'  +/- ',F10.5,/,
     +        5X,'YAVG = ',F12.5,'  +/- ',F10.5,/,
     +        5X,'SXX = ',F12.5,'  SYY = ',F12.5,'  SXY = ',F12.5,/,
     +        5X,'PEARSON CORRELATION COEFFICIENT = ', F10.5)
   30  FORMAT(13X,'SIX LINEAR REGRESSIONS: ANALYTICAL RESULTS')
   32  FORMAT(7X,'A = INTERCEPT,  B = SLOPE,  SD = STANDARD DEVIATION')
   34  FORMAT(29X,'A',7X,'SD(A)',7X,'B',7X,'SD(B)')
   36  FORMAT(' OLS(Y/X)          : ',4F10.3)
   38  FORMAT(' OLS(X/Y)          : ',4F10.3)
   40  FORMAT(' OLS BISECTOR      : ',4F10.3)
   42  FORMAT(' ORTHOGONAL        : ',4F10.3)
   44  FORMAT(' REDUCED MAJ AXIS  : ',4F10.3)
   46  FORMAT(' MEAN OLS          : ',4F10.3)
  105  FORMAT(5X,'BOOTSTRAP SIMULATION RESULTS: ',I5,' SIMULATIONS WITH
     + SEED:',I5,/,
     +   '   (REGRESSION COEFFICIENT AVERAGES AND STANDARD DEVIATIONS)')
  155  FORMAT(17X,'JACKKNIFE SIMULATION RESULTS: ',/,
     +   '   (REGRESSION COEFFICIENT AVERAGES AND STANDARD DEVIATIONS)')
  205  FORMAT(7X,' NORMAL SIMULATION RESULTS: ',I5,' SIMULATIONS',/,
     +   '   (REGRESSION COEFFICIENT AVERAGES AND STANDARD DEVIATIONS)')
  207  FORMAT(' AVERAGE SIMULATED RHO AND STANDARD DEV: ',2F8.4)
   99  FORMAT(10X,I5,4D10.3)
  999  FORMAT(1X,2F10.3)
C

       CLOSE(UNIT=10)
       RETURN  ! Added by Alister Graham
c       STOP
       END
C*******************************************************************
C************************* SUBROUTINE DATSTT ***********************
C*******************************************************************
C
       SUBROUTINE DATSTT(NumMAX,N,X,Y,XAVG,VARX,YAVG,VARY,VARXY,RHO)
C
C     *  SUBROUTINE TO COMPUTE SIMPLE STATISTICAL PROPERTIES OF A  *
C     *  BIVARIATE DATA SET.  IT GIVES THE VARIANCE IN X, VARIANCE *
C     *  IN Y, AND PEARSON'S LINEAR CORRELATION COEFFICIENT.       *
C
C
       IMPLICIT REAL*8(A-H,O-Z)
       DIMENSION X(NumMAX),Y(NumMAX)
C
C
C     *         INITIALIZATIONS
C
       S1 = 0.0
       S2 = 0.0
       SXX = 0.0
       SYY = 0.0
       SXY = 0.0
       RN = DFLOAT(N)
C
C     *         COMPUTE AVERAGES AND SUMS
C
       DO 100 I=1,N
         S1 = S1 + X(I)
         S2 = S2 + Y(I)
  100  CONTINUE
       XAVG = S1/RN
       YAVG = S2/RN
       DO 200 I = 1, N
       SXX  = SXX  + (X(I) - XAVG)**2
       SYY  = SYY  + (Y(I) - YAVG)**2
       SXY  = SXY  + (X(I) - XAVG)*(Y(I) - YAVG)
 200   CONTINUE
       IF(SXY.EQ.0.0) THEN
         WRITE(*,991)
         STOP
 991   FORMAT(5X,'SXY .EQ. ZERO.  DATSTT TERMINATED.')
       ENDIF
C
C     *         COMPUTE AND RETURN RESULTS
C
       VARX = SXX/RN
       VARY = SYY/RN
       VARXY = SXY/RN
       RHO = SXY/(DSQRT(SXX*SYY))
C
       RETURN
       END
C************************************************************************
C************************* SUBROUTINE SIXLIN ****************************
C************************************************************************
C
       SUBROUTINE SIXLIN(NumMAX,N,X,Y,A,SIGA,B,SIGB)
C     *
C     *                     SIX LINEAR REGRESSIONS
C     *     WRITTEN BY T. ISOBE, G. J. BABU AND E. D. FEIGELSON
C     *               CENTER FOR SPACE RESEARCH, M.I.T.
C     *                             AND
C     *              THE PENNSYLVANIA STATE UNIVERSITY
C     *
C     *                   REV. 1.0,   SEPTEMBER 1990
C     *
C     *       THIS SUBROUTINE PROVIDES LINEAR REGRESSION COEFFICIENTS
C     *    COMPUTED BY SIX DIFFERENT METHODS DESCRIBED IN ISOBE,
C     *    FEIGELSON, AKRITAS, AND BABU 1990, ASTROPHYSICAL JOURNAL
C     *    AND BABU AND FEIGELSON 1990, SUBM. TO TECHNOMETRICS.
C     *    THE METHODS ARE OLS(Y/X), OLS(X/Y), OLS BISECTOR, ORTHOGONAL,
C     *    REDUCED MAJOR AXIS, AND MEAN-OLS REGRESSIONS.
C     *
C     *    INPUT
C     *         X(I) : INDEPENDENT VARIABLE
C     *         Y(I) : DEPENDENT VARIABLE
C     *            N : NUMBER OF DATA POINTS
C     *
C     *    OUTPUT
C     *         A(J) : INTERCEPT COEFFICIENTS
C     *         B(J) : SLOPE COEFFICIENTS
C     *      SIGA(J) : STANDARD DEVIATIONS OF INTERCEPTS
C     *      SIGB(J) : STANDARD DEVIATIONS OF SLOPES
C     *     WHERE J = 1, 6.
C     *
C     *    ERROR RETURNS
C     *         CALCULATION IS STOPPED WHEN DIVISION BY ZERO WILL
C     *         OCCUR (I.E. WHEN SXX, SXY, OR EITHER OF THE OLS
C     *         SLOPES IS ZERO).
C
C
       IMPLICIT REAL*8(A-H,O-Z)
       DIMENSION X(NumMAX),Y(NumMAX)
       DIMENSION A(6),B(6),SIGA(6),SIGB(6)
C
C
C     *         INITIALIZATIONS
C
       S1 = 0.0
       S2 = 0.0
       SXX = 0.0
       SYY = 0.0
       SXY = 0.0
       SUM1 = 0.0
       SUM2 = 0.0
       SUM3 = 0.0
       RN = N
       DO 50 J = 1, 6
         A(J)    = 0.0
         SIGA(J) = 0.0
         B(J)    = 0.0
         SIGB(J) = 0.0
   50  CONTINUE
C
C     *         COMPUTE AVERAGES AND SUMS
C
       DO 100 I=1,N
         S1 = S1 + X(I)
         S2 = S2 + Y(I)
  100  CONTINUE
       XAVG = S1/RN
       YAVG = S2/RN
       DO 200 I = 1, N
       X(I) = X(I) - XAVG
       Y(I) = Y(I) - YAVG
       SXX  = SXX  + X(I)**2
       SYY  = SYY  + Y(I)**2
       SXY  = SXY  + X(I)*Y(I)
 200   CONTINUE
       IF(SXY.EQ.0.0) THEN
         WRITE(*,992) (X(I),Y(I),I=1,N)
 992     FORMAT(2F10.3)
         WRITE(*,991)
         STOP
 991   FORMAT(5X,'SXY IS ZERO.  SIXLIN TERMINATED.')
       ENDIF
       SIGN = 1.0
       IF(SXY .LT. 0.0) SIGN = -1.0
C
C
C     *               COMPUTE THE SLOPE COEFFICIENTS
C
C
       B(1) = SXY/SXX
       B(2) = SYY/SXY
       B(3) = (B(1)*B(2) - 1.0
     +  + DSQRT((1.0 + B(1)**2)*(1.0 + B(2)**2)))/(B(1) + B(2))
       B(4) = 0.5*(B(2) - 1.0/B(1)
     +                 + SIGN*DSQRT(4.0 + (B(2) - 1.0/B(1))**2))
       B(5) = SIGN*DSQRT(B(1)*B(2))
       B(6) = 0.5*(B(1) + B(2))
C
C
C     *            COMPUTE INTERCEPT COEFFICIENTS
C
C
       DO 300 J = 1, 6
         A(J) = YAVG - B(J)*XAVG
 300   CONTINUE
C
C
C      *     PREPARE FOR COMPUTATION OF VARIANCES
C
C
       GAM1 = B(3)/((B(1) + B(2))
     +         *DSQRT((1.0 + B(1)**2)*(1.0 + B(2)**2)))
       GAM2 = B(4)/(DSQRT(4.0*B(1)**2 + (B(1)*B(2) - 1.0)**2))
       DO 400 I = 1, N
        SUM1 = SUM1 + (X(I)*(Y(I) - B(1)*X(I)))**2
        SUM2 = SUM2 + (Y(I)*(Y(I) - B(2)*X(I)))**2
        SUM3 = SUM3 + X(I)*Y(I)*(Y(I) - B(1)*X(I))*(Y(I) - B(2)*X(I))
  400  CONTINUE
       COV = SUM3/(B(1)*SXX**2)
C
C
C     *    COMPUTE VARIANCES OF THE SLOPE COEFFICIENTS
C
C
       SIGB(1) = SUM1/(SXX**2)
       SIGB(2) = SUM2/(SXY**2)
       SIGB(3) = (GAM1**2)*(((1.0 + B(2)**2)**2)*SIGB(1)
     +               + 2.0*(1.0 + B(1)**2)*(1.0 + B(2)**2)*COV
     +               + ((1.0 +B(1)**2)**2)*SIGB(2))
       SIGB(4) = (GAM2**2)*(SIGB(1)/B(1)**2 + 2.0*COV
     +               + B(1)**2*SIGB(2))
       SIGB(5) = 0.25*(B(2)*SIGB(1)/B(1)
     +                       + 2.0*COV + B(1)*SIGB(2)/B(2))
       SIGB(6) = 0.25*(SIGB(1) + 2.0*COV + SIGB(2))
C
C
C     *   COMPUTE VARIANCES OF THE INTERCEPT COEFFICIENTS
C
C
       DO 500 I = 1, N
          SIGA(1) = SIGA(1) + ((Y(I) - B(1)*X(I))
     +                        *(1.0 - RN*XAVG*X(I)/SXX))**2
          SIGA(2) = SIGA(2) + ((Y(I) - B(2)*X(I))
     +                        *(1.0 - RN*XAVG*Y(I)/SXY))**2
          SIGA(3) = SIGA(3) + ((X(I)*(Y(I)
     +            - B(1)*X(I))*(1.0 + B(2)**2)/SXX
     +            + Y(I)*(Y(I) - B(2)*X(I))*(1.0 + B(1)**2)/SXY)
     +              *GAM1*XAVG*RN - Y(I) + B(3)*X(I))**2
          SIGA(4) = SIGA(4) + ((X(I)*(Y(I) - B(1)*X(I))/SXX
     +            + Y(I)*(Y(I) - B(2)*X(I))*(B(1)**2)/SXY)*GAM2
     +               *XAVG*RN/DSQRT(B(1)**2) - Y(I) + B(4)*X(I))**2
          SIGA(5) = SIGA(5) + ((X(I)*(Y(I)
     +            - B(1)*X(I))*DSQRT(B(2)/B(1))/SXX
     +            + Y(I)*(Y(I) - B(2)*X(I))*DSQRT(B(1)/B(2))/SXY)
     +              *0.5*RN*XAVG - Y(I) + B(5)*X(I))**2
          SIGA(6) = SIGA(6) + ((X(I)*(Y(I) - B(1)*X(I))/SXX
     +            + Y(I)*(Y(I) - B(2)*X(I))/SXY)
     +              *0.5*RN*XAVG - Y(I) + B(6)*X(I))**2
  500  CONTINUE
C
C
C     *  CONVERT VARIANCES TO STANDARD DEVIATIONS
C
C
       DO 600 J = 1, 6
          SIGB(J) = DSQRT(SIGB(J))
          SIGA(J) = DSQRT(SIGA(J))/RN
  600  CONTINUE
C
C
C     *  RETURN DATA ARRAYS TO THEIR ORIGINAL FORM
C
C
       DO 900 I = 1, N
          X(I) = X(I) + XAVG
          Y(I) = Y(I) + YAVG
  900  CONTINUE
C
C
       RETURN
       END
C************************************************************************
C************************ SUBROUTINE BOOTSP *****************************
C************************************************************************
C
       SUBROUTINE BOOTSP(NumMAX,N,X,Y,XBOOT,YBOOT,ISEED)
C     * CONSTRUCTS MONTE CARLO SIMULATED DATA SET USING THE     *
C     * BOOTSTRAP ALGORITHM.                                    *
C     *                                                         *
C     * CALLS FUNCTION SUBROUTINE RANDOM                        *
C     *                                                         *
       IMPLICIT REAL*8(A-H,O-Z)
       REAL*4 RANDOM
       DIMENSION X(NumMAX),Y(NumMAX),XBOOT(NumMAX),YBOOT(NumMAX)
C
       DO 10 I = 1,N
          J = INT(RANDOM(ISEED)*FLOAT(N) + 1.0)
          XBOOT(I) = X(J)
          YBOOT(I) = Y(J)
  10   CONTINUE
       RETURN
       END
C************************************************************************
C************************ FUNCTION SUBROUTINE RANDOM ********************
C************************************************************************
C
C    * RETURNS UNIFORM RANDOM NUMBER BETWEEN 0.0 AND 1.0.           *
C    * CALLS SUBROUTINE IN55, WHICH CALLS FUNCTION IRN55,           *
C    * CODE FROM D. E. KNUTH, `THE ART OF SCIENTIFIC PROGRAMMING'   *
C    * 2ND ED., VOL. 2, PP. 170-3, 1973. SET ISEED < 0 TO START.    *
C    *                                                              *
C
      FUNCTION RANDOM(ISEED)
      DIMENSION IA(55)
      COMMON /RAND/ JRAND
C
C         INITIALIZATION
C
      IF(ISEED.LT.0) THEN
         JRAND = 55
         CALL IN55(IA,ISEED)
         ISEED = IABS(ISEED)
      ENDIF
      JRAND = JRAND + 1
      IF(JRAND .GT. 55) JRAND = IRN55(IA)
      RANDOM = FLOAT(IA(JRAND)) * 1.0E-9
      RETURN
      END
C************************************************************************
C************************* SUBROUTINE IN55 ******************************
C************************************************************************
C
      SUBROUTINE IN55(IA,IX)
      DIMENSION IA(55)
C
C    * THIS SUBROUTINE SETS IA(1), ..., IA(55) TO STARTING
C    * VALUES SUITABLE FOR LATER CALLS ON IRN55(IA).
C    * IX IS AN INTEGER "SEED" VALUE BETWEEN 0 AND 999999999.
C
      IA(55) = IX
      J = IX
      K = 1
      DO 1 I = 1, 54
        II = MOD(21*I, 55)
        IA(II) = K
        K = J - K
        IF (K .LT. 0) K = K + 1000000000
        J = IA(II)
   1  CONTINUE
C    * THE NEXT THREE LINES "WARM UP" THE GENERATOR.
      JDUM = IRN55(IA)
      JDUM = IRN55(IA)
      JDUM = IRN55(IA)
      RETURN
      END
C************************************************************************
C************************ FUNCTION IRN55 ********************************
C************************************************************************
C
      FUNCTION IRN55(IA)
      DIMENSION IA(1)
C
C     ASSUMING THAT IA(1), ..., IA(55) HAVE BEEN SET UP PROPERLY,
C     THIS SUBROUTINE RESETS THE IA ARRAY TO THE NEXT 55 NUMBERS
C     OF A PSEUDO-RANDOM SEQUENCE, AND RETURNS THE VALUE 1
C
      DO 1 I = 1, 24
        J = IA(I) - IA(I+31)
        IF (J .LT. 0) J = J + 1000000000
        IA(I) = J
   1  CONTINUE
      DO 2 I = 25, 55
        J = IA(I) - IA(I-24)
        IF (J .LT. 0) J = J + 1000000000
        IA(I) = J
   2  CONTINUE
      IRN55 = 1
      RETURN
      END
