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=<A9> *
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
PARAMETER (NMAX = 500)
IMPLICIT REAL*8 (A-H,O-Z)
CHARACTER*80 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)
C
C * READ THE FILE NAME IN WHICH A DATA SET EXISTS.
C
WRITE(*,10)
READ(*,12) FILE
C
OPEN (UNIT=50, FILE=FILE, STATUS='OLD', FORM='FORMATTED')
C
C * READ THE DATA SET AND OPEN OUTPUT FILE
C
I=0
15 I=I+1
X(I) = 0.0
Y(I) = 0.0
c READ(50,14,END=25) X(I),Y(I)
READ(50,*,END=25) X(I),Y(I)
GOTO 15
25 CONTINUE
NTOT=I-1
CLOSE (UNIT = 50)
OPEN (UNIT=10, FILE='SLOPES.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(NMAX,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(NMAX,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
C * MAKE BOOTSTRAP SIMULATED DATASETS, AND COMPUTE AVERAGES
C * AND STANDARD DEVIATIONS OF REGRESSION COEFFICIENTS
C
WRITE(*,101)
READ(*,102) NSIM
WRITE(*,103)
READ(*,102) ISEED
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(NMAX,NTOT,X,Y,XSIM,YSIM,ISEED)
CALL SIXLIN(NMAX,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(NMAX,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)
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(NMAX,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(' DATA FILE : ',$)
101 FORMAT(' NUMBER OF SIMULATIONS: ')
102 FORMAT(BN,I6)
103 FORMAT(' SEED (INTEGER, I6): ')
12 FORMAT(A80)
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: ',A9,' # 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
STOP
END
C*******************************************************************
C************************* SUBROUTINE DATSTT ***********************
C*******************************************************************
C
SUBROUTINE DATSTT(NMAX,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(NMAX),Y(NMAX)
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(NMAX,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(NMAX),Y(NMAX)
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(NMAX,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(NMAX),Y(NMAX),XBOOT(NMAX),YBOOT(NMAX)
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(1)
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