C FORTRAN program maxsd.for Version 2.0 4/11/01 C Tests toxicity data on dose levels d_1,...,d_k compared with C a zero dose d_0 to determine the estimated maximum safe dose. C C by Charles W. Dunnett C Mathematics & Statistics Department C McMaster University C Hamilton ON, L8S 4K1, Canada C dunnett@mcmaster.ca C C Please direct all comments and questions about this program to C the author at the above address, or to C C Ajit C. Tamhane C Department of Industrial Engineering C & Management Sciences C Northwestern University C Evanston, IL 60208 C ajit@iems.northwestern.edu C C This program may be freely used and distributed for non-commercial C purposes. C C This FORTRAN program, MAXSD, requires the program MVSTUD and C other subroutines from Applied Statistics algorithm 251 (Dunnett C 1989) to compute adjusted p-values; these programs are included C with MAXSD. These latter programs, including descriptive information, C can also be obtained separately from the Statlib at Carnegie Mellon C University by sending an e-mail to statlib@lib.stat.cmu.edu and writing: C C send 251 from apstat. C C MAXSD enters data either manually from the keyboard or from a disk C file in the following format: C C line 1: k, lambda, alpha C line 2: n_0, n_1, ..., n_k (sample sizes) C line 3: ybar_0, ybar_1, ..., ybar_k (sample means) C line 4: variance estimate, d.f. C C Here, lambda is a parameter set by the user, representing the C proportion of the expected response at the zero dose which is C considered to be unsafe. Hypotheses H_0i, ..., H_0k are tested, C where H_0i: mu_i <= lambda mu_0 means the i-th dose is unsafe. C The H_0i are tested in order from i=1 to i=k, stopping the first C time an H_0i fails to be rejected, in which case the preceding C dose level is declared to be the maximum safe dose. The H_0i are C tested with familywise error rate (FWE) <= alpha using contrasts C in the sample means as test statistics. C C Users can choose between Helmert contrasts and pairwise contrasts C to use in the testing, with a choice of 2 types of pairwise contrasts C denoted SD1 and SD2. The program performs all computations and displays C the results which identify the estimated maximum safe dose. C C A typical value of lambda is 0.85, which means that a 15% decrease C in response from the zero dose is considered to be unsafe. C C More detailed information is given in the paper by Tamhane et al C (2001), which can be obtained from the authors. C C REFERENCES C C Dunnett, C. W. (1989), "Algorithm AS 251: Multivariate Normal C Probability Integrals with Product Correlation Structure," C Applied Statistics, 38, 564 - 579. C C Tamhane, A. C., Dunnett, C. W., Green, J. W. and Wetherington, J. D. C (2001), "Multiple Test Procedures for Identifying the Maximum Safe C Dose," Journal of the American Statistical Association, 96. To C appear. C C PROGRAM MAXSD C C OUTPUT: C T-STATISTICS C ADJUSTED P-VALUES FOR EACH DOSE C ESTIMATED MAXSD FOR SPECIFIED ALPHA C C OUTPUT IS STORED IN A FILE TO BE NAMED BY THE USER C AT THE PROMPT C C PROGRAMMED BY: C.W. DUNNETT, MCMASTER UNIVERSITY C DIMENSION X(10), YBAR(10), XN(10), REVYB(10), REVXN(10) * , CTR(10), XSE(10) CHARACTER*64 FNAME, DNAME COMMON C(10,10) DATA ZERO, ONE / 0.0, 1.0 / J0 = 0 WRITE (*, ' (A\)') ' Enter file name for storing output: ' READ (*, '(A)') FNAME OPEN (4, FILE = FNAME ) 5 WRITE (*, 410) READ (*, *) MCTRST C IF (MCTRST .LT. 1 .OR. MCTRST .GT. 4) GO TO 5 IF (MCTRST .LT. 1 .OR. MCTRST .GT. 3) GO TO 5 WRITE (4, 350) WRITE (*, 370) READ (*, *) IPRNT WRITE (*, 400) READ (*, *) IDATA IF (IDATA .NE. 1 .AND. IDATA .NE. 2) STOP IF (IDATA .EQ. 2) GO TO 30 10 WRITE (*, 200) READ (*, *) K, XLAM, ALPHA IF (K .LE. 0) STOP IF (XLAM .LE. ZERO .OR. XLAM .GT. ONE) THEN WRITE (*, 440) GO TO 10 ENDIF WRITE (4, 210) K, ALPHA WRITE (*, 220) READ (*, *) XN0, (XN(I), I = 1, K) WRITE (4, 230) XN0, (XN(I), I = 1, K) WRITE (*, 270) READ (*, *) YBAR0, (YBAR(I), I = 1, K) IF (YBAR0 .LE. 0.0) STOP WRITE (4, 260) YBAR0, (YBAR(I), I = 1, K) WRITE (*, 240) READ (*, *) SIGMA2, NDF SIGMA = SQRT ( SIGMA2 ) DF = ZERO IF (NDF .GT. 0) DF = FLOAT (NDF) WRITE (4, 240) WRITE (4, 280) SIGMA2, NDF WRITE (4, 490) MCTRST IF (MCTRST .EQ. 1) WRITE (4, 491) IF (MCTRST .EQ. 2) WRITE (4, 492) IF (MCTRST .EQ. 3) WRITE (4, 493) IF (MCTRST .EQ. 4) WRITE (4, 494) GO TO 39 C C Enter data from disk file C 30 WRITE (*, ' (A\)') ' Enter file name for input data: ' READ (*, '(A)') DNAME OPEN (3, FILE = DNAME ) JOBN = 0 35 READ (3, *) K, XLAM, ALPHA IF (K .LE. 0) STOP JOBN = JOBN + 1 IF (JOBN .GT. 1) WRITE (4, 310) READ (3, *) XN0, (XN(I), I = 1, K) WRITE (*, 230) XN0, (XN(I), I = 1, K) WRITE (4, 230) XN0, (XN(I), I = 1, K) READ (3, *) YBAR0, (YBAR(I), I = 1, K) WRITE (*, 260) YBAR0, (YBAR(I), I = 1, K) WRITE (4, 260) YBAR0, (YBAR(I), I = 1, K) READ (3, *) SIGMA2, NDF SIGMA = SQRT ( SIGMA2 ) DF = ZERO IF (NDF .GT. 0) DF = FLOAT (NDF) WRITE (*, 280) SIGMA2, NDF WRITE (4, 280) SIGMA2, NDF WRITE (*, 210) K, ALPHA WRITE (4, 210) K, ALPHA 38 WRITE (4, 490) MCTRST IF (MCTRST .EQ. 1) WRITE (4, 491) IF (MCTRST .EQ. 2) WRITE (4, 492) IF (MCTRST .EQ. 3) WRITE (4, 493) IF (MCTRST .EQ. 4) WRITE (4, 494) 39 WRITE (4, 450) XLAM WRITE (*, 450) XLAM C C REVERSE LABELS C 40 CONTINUE DO 45 I = 1, K IX = K + 1 - I REVYB(I) = YBAR(IX) REVXN(I) = XN(IX) 45 CONTINUE C C DETERMINE P-VALUES FOR EACH DOSE, C AND IXX = ESTIMATED MAXSD FOR SPECIFIED ALPHA LEVEL, C WHERE MAXSD = MIN. DOSE WITH RESPONSE < LAMBDA * MU_0 C MINUS 1, USING SD3S PROCEDURE C PVAL = ZERO IXMUD = 0 C C IXMUD = minimum dose level that is unsafe C (i.e., p-value for tmax exceeds alpha) C DO 130 I = 1, K IX = K + 1 - I GO TO (51, 52, 52, 52), MCTRST 51 CALL CONTRAST ( IX, YBAR0, REVYB, X, XLAM, XN0, REVXN * , SIGMA, RHOB, CTR, XSE) GO TO 55 52 CALL PAIR2 ( IX, YBAR0, REVYB, X, XLAM, XN0, REVXN, * SIGMA, RHOB, CTR, XSE) 55 CONTINUE IF (IPRNT .EQ. 0) GO TO 90 C C PRINT MATRIX OF CONTRAST COEFFICIENTS C WRITE (4, 360) I WRITE (4, 330) J0, (J1, J1 = K - IX + 1, K) WRITE (4, 235) 90 LOW = 1 IF (MCTRST .EQ. 3) LOW = IX IF (IPRNT .EQ. 0) GO TO 115 DO 110 IT = LOW, IX WRITE (4, 235) IT, C(IT,1),(C(IT,IX+3-J1), J1=2,IX+1) 110 CONTINUE WRITE (4, 250) XLAM * YBAR0, (YBAR(J1), J1 = I, K) 115 TMAX = -100. WRITE (*, 340) I, (X(J1), J1 = 1, IX) WRITE (4, 340) I IF (RHOB .LT. ZERO) RHOB = ZERO LX = IX IF (MCTRST .EQ. 3) LX = 1 IF (MCTRST .EQ. 4) LX = K WRITE (4, 380) DO 118 J1 = LOW, IX SE = XSE(J1) CALL SUBTINT(LX,RHOB,X(J1),PROB,BND,IFLT, NDF) PROB = ONE - PROB WRITE (4, 390 ) J1, CTR(J1), SE, X(J1), PROB 118 CONTINUE DO 120 J1 = LOW, IX IF (X(J1) .GT. TMAX) TMAX = X(J1) 120 CONTINUE IF (RHOB .LT. ZERO) RHOB = ZERO LX = IX IF (MCTRST .EQ. 3) LX = 1 IF (MCTRST .EQ. 4) LX = K CALL SUBTINT(LX,RHOB,TMAX,PROB,BND,IFLT, NDF) PROB = ONE - PROB IF (PROB .LT. PVAL) PROB = PVAL PVAL = PROB IF (PVAL .GT. ALPHA .AND. IXMUD .EQ. 0) IXMUD = I WRITE (*, 320) TMAX, RHOB, PVAL WRITE (4, 320) TMAX, RHOB, PVAL 130 CONTINUE IF (IXMUD .EQ. 0) IXMUD = K + 1 IXMAXSD = IXMUD - 1 WRITE (*, 300) IXMAXSD, ALPHA WRITE (4, 300) IXMAXSD, ALPHA IF (IDATA .EQ. 2) GO TO 35 WRITE (*, 420) READ (*, *) IX1 IF (IX1 .LE. 0 .OR. IX1 .GT. 3) STOP GO TO (10, 140, 150), IX1 140 WRITE (*, 410) READ (*, *) MCTRST IF (MCTRST .LT. 1 .OR. MCTRST .GT. 4) GO TO 140 GO TO 38 150 WRITE (*, 440) READ (*, *) XLAM IF (XLAM .LE. ZERO .OR. XLAM .GT. ONE) GO TO 150 GO TO 39 200 FORMAT (/, 1X, 'Enter k, lambda and alpha') 210 FORMAT (/, 1X, 'k =', I4, 5X, 'alpha =', F6.3) 220 FORMAT (/, 1X, 'Enter n_0, n_1, ..., n_k') 230 FORMAT (/, 1X, 'Sample sizes:', 10F9.0) 235 FORMAT (6X, I4, 1X, 10F8.0) 330 FORMAT (5X, ' j', 10I8) 240 FORMAT (/, 1X, 'Enter variance estimate, and d.f.' * /, 1X, ' (enter 0 for inf. d.f.)' ) 250 FORMAT (11X, 10F8.4) 260 FORMAT (1X, 'Means: ', 10F9.4) 270 FORMAT (/, 1X, 'Enter ybar_0, ybar_1, ... , ybar_k') 300 FORMAT (/,1X,'Est. MAXSD =',I2,' (alpha =', F4.2, ')', /) 280 FORMAT (5X, 'sigma^2 =', F17.5, 4X, 'd.f. =', I6) 290 FORMAT (/, 1X, 'Enter 1 for new problem, 0 to stop') 310 FORMAT (//, 1X, 'NEW PROBLEM:-') 340 FORMAT (/, 1X, 't-values for dose', I2,': ', 10(F8.3,1X)) 320 FORMAT (3X, 't_max =', F8.3, ' rhob =', F7.4 * , ' Adj. p-value =', F5.3, /) 390 FORMAT ( 6X, I4, F11.4, F11.5, F11.3, F12.3) 400 FORMAT (/, 1X, 'Enter 1 to enter data from keyboard' * /, 1X, ' or 2 to enter data from disk file' * /, 1X, ' or 0 to cancel') 410 FORMAT (/, 1X, 'Enter 1 for SD1 Helmert contrasts' * /, 1X, ' or 2 for SD1 pairwise contrasts' * /, 1X, ' or 3 for SD2 pairwise contrasts') C * /, 1X, ' or 4 for SS1 pairwise contrasts') 420 FORMAT (/, 1X, 'Enter 1 to start a new problem' * /, 1X, ' or 2 to change the method' * /, 1X, ' or 3 to change the value of lambda' * /, 1X, ' or 0 to stop') 440 FORMAT (/, 1X, 'Enter lambda: 0 < lambda <= 1') 450 FORMAT (/, 1X, 'lambda = ', F4.2) 350 FORMAT (1X, 'Formulation: MAXSD, or Max. Safe Dose:-' ) 360 FORMAT (/, 1X, 'Contrast coefficients for testing dose', I2, * ':-' /) 370 FORMAT (/, 1X, 'Enter 1 to print contrasts, 0 to omit') 380 FORMAT (/, 9X, 'j contrast', 7X, 's.e.', 9X, 't', 9X, 'prob',/) 490 FORMAT (/, 1X, 'Method Type = ', I3) 491 FORMAT ( 1X, '(SD1HC - SD1 Helmert contrasts)') 492 FORMAT ( 1X, '(SD1PC - SD1 prws. contrasts)') 493 FORMAT ( 1X, '(SD2PC - SD2 prws. contrasts)') 494 FORMAT ( 1X, '(SS1PC - SS1 prws. contrasts)') END SUBROUTINE CONTRAST (K, W0, RW, X, XLAM, XN0, RXN, SIGMA, RHOB, * CTR, XSE) C C THE CONTRASTS FOR K = 5 ARE AS FOLLOWS C FOR DECREASING DOSE RESPONSE:- C contrast d_0 d_1 d_2 d_3 d_4 d_5 C 1 -5 1 1 1 1 1 C 2 -4 1 1 1 1 0 C 3 -3 1 1 1 0 0 C 4 -2 1 1 0 0 0 C 5 -1 1 0 0 0 0 C DIMENSION X(*), RW(*), RXN(*), SS(10,10), CTR(10), RHO(10,10) * , XSE(*) COMMON C(10,10) DATA ZERO, ONE / 0.0, 1.0 / DO 20 I = 1, K IX = K + 1 - I FLI = -FLOAT(IX) X(I) = FLI * XLAM * W0 SSS = (FLI * XLAM) ** 2 / XN0 C(I, 1) = FLI DO 10 J = 1, K IF ( J .GE. I ) THEN X(I) = X(I) + RW(J) C(I, J + 1) = ONE SSS = SSS + ONE / RXN(J) ELSE C(I, J + 1) = ZERO ENDIF 10 CONTINUE C C CTR = CONTRAST VALUE C XSE = CONTRAST STD. ERR. C X = CONTRAST T-VALUE C CTR(I) = X(I) XSE(I) = SQRT(SSS) * SIGMA X(I) = X(I) / XSE(I) 20 CONTINUE C C COMPUTE RHOB = AVG. RHO C RHOB = ZERO XNR = ZERO IF (K .EQ. 1) RETURN DO 40 I = 1, K DO 40 J = 1, I SS(I, J) = XLAM * XLAM * C(I,1) * C(J,1) / XN0 DO 40 M = 2, K+1 SS(I, J) = SS(I, J) + C(I,M) * C(J,M) / RXN(M-1) 40 CONTINUE DO 50 I = 2, K DO 50 J = 1, I-1 RHO(I, J) = SS(I, J) / SQRT(SS(I, I) * SS(J, J)) RHOB = RHOB + RHO(I, J) XNR = XNR + ONE 50 CONTINUE RHOB = RHOB / XNR WRITE (4, 100) DO 60 I = 2, K WRITE (4, 90) (RHO(I, J), J = 1, I-1) 60 CONTINUE RETURN 90 FORMAT ( 6X, 10F8.5) 100 FORMAT (/, 1X, 'Correlation coefficients:- ', /) END SUBROUTINE SUBTINT (N, RHO, H, PROB, BOUND, IFAULT, NDF) C C PROBABILITY INTEGRAL FOR MULTIVARIATE NORMAL OR MULTIV. T C WITH COMMON CORRELATION C DIMENSION A(50), B(50), D(50), BPD(50), INF(50) DATA ZERO, EPS / 0.0, 0.00001/ HINC = ZERO IERC = 0 XN = FLOAT(N) ITYPE = 1 SQTRHO = ZERO IF (RHO .GT. ZERO) SQTRHO = SQRT(RHO) DO 10 I = 1, N INF(I) = ITYPE BPD(I) = SQTRHO D(I) = ZERO A(I) = H B(I) = -H 10 CONTINUE C COMPUTE THE MULTIVARIATE T PROBABILITY CALL MVSTUD(NDF,A,B,BPD,EPS,N,INF,D,IERC,HINC,PROB, * BOUND,IFAULT) RETURN END SUBROUTINE PAIR2 (K, W0, RW, X, XLAM, XN0, RXN, SIGMA, RHOB, * CTR, XSE) C C THE CONTRASTS FOR K = 5 ARE AS FOLLOWS C FOR MAXSD:- C contrast d_0 d_1 d_2 d_3 d_4 d_5 C 1 -1 0 0 0 0 1 C 2 -1 0 0 0 1 0 C 3 -1 0 0 1 0 0 C 4 -1 0 1 0 0 0 C 5 -1 1 0 0 0 0 C DIMENSION X(*), RW(*), RXN(*), SS(10,10), CTR(10), XSE(*) COMMON C(10,10) DATA ZERO, ONE / 0.0, 1.0 / DO 20 I = 1, K IX = 1 FLI = -FLOAT(IX) X(I) = FLI * XLAM * W0 SSS = (FLI * XLAM) ** 2 / XN0 C(I, 1) = FLI DO 10 J = 1, K IF ( J .EQ. I ) THEN X(I) = X(I) + RW(J) C(I, J + 1) = ONE SSS = SSS + ONE / RXN(J) ELSE C(I, J + 1) = ZERO ENDIF 10 CONTINUE CTR(I) = X(I) XSE(I) = SQRT(SSS) * SIGMA X(I) = X(I) / XSE(I) 20 CONTINUE C C COMPUTE RHOB = AVG. RHO C RHOB = ZERO XNR = ZERO IF (K .EQ. 1) RETURN DO 40 I = 1, K DO 40 J = 1, I SS(I, J) = XLAM * XLAM * C(I,1) * C(J,1) / XN0 DO 40 M = 2, K+1 SS(I, J) = SS(I, J) + C(I,M) * C(J,M) / RXN(M-1) 40 CONTINUE DO 50 I = 2, K DO 50 J = 1, I-1 RHOB = RHOB + SS(I, J) / SQRT(SS(I, I) * SS(J, J)) XNR = XNR + ONE 50 CONTINUE RHOB = RHOB / XNR RETURN END SUBROUTINE MVNPRD(A, B, BPD, EPS, N, INF, IERC, HINC, PROB, BOUND, * IFAULT) C C ALGORITHM AS 251.1 APPL.STATIST. (1989), VOL.38, NO.3 C C FOR A MULTIVARIATE NORMAL VECTOR WITH CORRELATION STRUCTURE C DEFINED BY RHO(I,J) = BPD(I) * BPD(J), COMPUTES THE PROBABILITY C THAT THE VECTOR FALLS IN A RECTANGLE IN N-SPACE WITH ERROR C LESS THAN EPS. C INTEGER NN PARAMETER (NN = 50) REAL A(*), B(*), BPD(*), ESTT(22), FV(5), FD(5), F1T(22), * F2T(22), F3T(22), G1T(22), G3T(22), PSUM(22), H(NN), HL(NN), * BB(NN) INTEGER INF(*), INFT(NN), LDIR(22) REAL ZERO, HALF, ONE, TWO, FOUR, SIX, PT1, PT24, ONEP5, * X2880, SMALL, DXMIN, SQRT2, PROB, ERRL, BI, START, * Z, HINC, ADDN, EPS2, EPS1, EPS, ZU, Z2, Z3, Z4, Z5, ZZ, * ERFAC, EL, EL1, BOUND, PART0, PART2, PART3, FUNC0, FUNC2, * FUNCN, WT, CONTRB, DLG, DX, DA, ESTL, ESTR, SUM, EXCESS, ERROR, * PROB1, SAFE INTEGER N, IERC, IFAULT, I, NTM, NMAX, LVL, NR, NDIM REAL ALNORM, PPND7 EXTERNAL ALNORM, PPND7 DATA ZERO, HALF, ONE, TWO, FOUR, SIX /0.0, 0.5, 1.0, 2.0, * 4.0, 6.0/ DATA PT1, PT24, ONEP5, X2880 /0.1, 0.24, 1.5, 2880.0/ DATA SMALL, DXMIN, SQRT2 /1.0E-10, 0.0000001, 1.41421356237310/ C C CHECK FOR INPUT VALUES OUT OF RANGE. C PROB = ZERO BOUND = ZERO IFAULT = 1 IF (N .LT. 1 .OR. N .GT. NN) RETURN DO 10 I = 1, N BI = ABS(BPD(I)) IFAULT = 2 IF (BI .GE. ONE) RETURN IFAULT = 3 IF (INF(I) .LT. 0 .OR. INF(I) .GT. 2) RETURN IFAULT = 4 IF (INF(I) .EQ. 2 .AND. A(I) .LE. B(I)) RETURN 10 CONTINUE IFAULT = 0 PROB = ONE C C CHECK WHETHER ANY BPD(I) = 0. C NDIM = 0 DO 20 I = 1, N IF (BPD(I) .NE. ZERO) THEN NDIM = NDIM + 1 H(NDIM) = A(I) HL(NDIM) = B(I) BB(NDIM) = BPD(I) INFT(NDIM) = INF(I) ELSE C C IF ANY BPD(I) = 0, THE CONTRIBUTION TO PROB FOR THAT C VARIABLE IS COMPUTED FROM A UNIVARIATE NORMAL. C IF (INF(I) .LT. 1) THEN PROB = PROB * (ONE - ALNORM(B(I), .FALSE.)) ELSE IF (INF(I) .EQ. 1) THEN PROB = PROB * ALNORM(A(I), .FALSE.) ELSE PROB = PROB * (ALNORM(A(I), .FALSE.) - * ALNORM(B(I), .FALSE.)) END IF IF (PROB .LE. SMALL) PROB = ZERO END IF 20 CONTINUE IF (NDIM .EQ. 0 .OR. PROB .EQ. ZERO) RETURN C C IF NOT ALL BPD(I) = 0, PROB IS COMPUTED BY SIMPSON'S RULE. C BUT FIRST, INITIALIZE THE VARIABLES. C Z = ZERO IF (HINC .LE. ZERO) HINC = PT24 ADDN = -ONE DO 30 I = 1, NDIM IF (INFT(I) .EQ. 2 .OR. * (INFT(I) .NE. INFT(1) .AND. BB(I) * BB(1) .GT. ZERO) .OR. * (INFT(I) .EQ. INFT(1) .AND. BB(I) * BB(1) .LT. ZERO)) * ADDN = ZERO 30 CONTINUE C C THE VALUE OF ADDN IS TO BE ADDED TO THE PRODUCT EXPRESSIONS IN C THE INTEGRAND TO INSURE THAT THE LIMITING VALUE IS ZERO. C PROB1 = ZERO NTM = 0 NMAX = 400 IF (IERC .EQ. 0) NMAX = NMAX * 2 CALL PFUNC (Z, H, HL, BB, NDIM, INFT, ADDN, SAFE, FUNC0, NTM, * IERC, PART0) EPS2 = EPS * PT1 * HALF C C SET UPPER BOUND ON Z AND APPORTION EPS. C ZU = -PPND7(EPS2, IFAULT) / SQRT2 IF (IFAULT .NE. 0) THEN IFAULT = 6 RETURN END IF NR = IFIX(ZU / HINC) + 1 ERFAC = ONE IF (IERC .NE. 0) ERFAC = X2880 / HINC ** 5 EL = (EPS - EPS2) / FLOAT(NR) * ERFAC EL1 = EL C C START COMPUTATIONS FOR THE INTERVAL (Z, Z + HINC). C 40 ERROR = ZERO LVL = 0 FV(1) = PART0 FD(1) = SAFE START = Z DA = HINC Z3 = START + HALF * DA CALL PFUNC(Z3, H, HL, BB, NDIM, INFT, ADDN, FD(3), FUNCN, NTM, * IERC, FV(3)) Z5 = START + DA CALL PFUNC(Z5, H, HL, BB, NDIM, INFT, ADDN, FD(5), FUNC2, NTM, * IERC, FV(5)) PART2 = FV(5) SAFE = FD(5) WT = DA / SIX CONTRB = WT * (FV(1) + FOUR * FV(3) + FV(5)) DLG = ZERO IF (IERC .NE. 0) THEN CALL WMAX(FD(1), FD(3), FD(5), DLG) IF (DLG .LE. EL) GO TO 90 DX = DA GO TO 60 END IF LVL = 1 LDIR(LVL) = 2 PSUM(LVL) = ZERO C C BISECT INTERVAL. IF IERC = 1, COMPUTE ESTIMATE ON LEFT C HALF; IF IERC = 0, ON BOTH HALVES. C 50 DX = HALF * DA WT = DX / SIX Z2 = START + HALF * DX CALL PFUNC(Z2, H, HL, BB, NDIM, INFT, ADDN, FD(2), FUNCN, NTM, * IERC,FV(2)) ESTL = WT * (FV(1) + FOUR * FV(2) + FV(3)) IF (IERC .EQ. 0) THEN Z4 = START + ONEP5 * DX CALL PFUNC(Z4, H, HL, BB, NDIM, INFT, ADDN, FD(4), FUNCN, * NTM, IERC, FV(4)) ESTR = WT * (FV(3) + FOUR * FV(4) + FV(5)) SUM = ESTL + ESTR DLG = ABS(CONTRB - SUM) EPS1 = EL / TWO ** (LVL - 1) ERRL = DLG ELSE FV(3) = FV(2) FD(3) = FD(2) CALL WMAX(FD(1), FD(3), FD(5), DLG) ERRL = DLG / TWO ** (5 * LVL) SUM = ESTL EPS1 = EL * (TWO ** LVL) ** 4 END IF C C STOP SUBDIVIDING INTERVAL WHEN ACCURACY IS SUFFICIENT, C OR IF INTERVAL TOO NARROW OR SUBDIVIDED TOO OFTEN. C IF (DLG .LE. EPS1 .OR. DLG .LT. SMALL) GO TO 70 IF (IFAULT .EQ. 0 .AND. NTM .GE. NMAX) IFAULT = 5 IF (ABS(DX) .LE. DXMIN .OR. LVL .GT. 21) IFAULT = 7 IF (IFAULT .NE. 0) GO TO 70 C C RAISE LEVEL. STORE INFORMATION FOR RIGHT HALF AND APPLY C SIMPSON'S RULE TO LEFT HALF. C 60 LVL = LVL + 1 LDIR(LVL) = 1 F1T(LVL) = FV(3) F3T(LVL) = FV(5) DA = DX FV(5) = FV(3) IF (IERC .EQ. 0) THEN F2T(LVL) = FV(4) ESTT(LVL) = ESTR CONTRB = ESTL FV(3) = FV(2) ELSE G1T(LVL) = FD(3) G3T(LVL) = FD(5) FD(5) = FD(3) END IF GO TO 50 C C ACCEPT APPROXIMATE VALUE FOR INTERVAL. C RESTORE SAVED INFORMATION TO PROCESS C RIGHT HALF INTERVAL. C 70 ERROR = ERROR + ERRL 80 IF (LDIR(LVL) .EQ. 1) THEN PSUM(LVL) = SUM LDIR(LVL) = 2 IF (IERC .EQ. 0) DX = DX * TWO START = START + DX DA = HINC / TWO ** (LVL - 1) FV(1) = F1T(LVL) IF (IERC .EQ. 0) THEN FV(3) = F2T(LVL) CONTRB = ESTT(LVL) ELSE FV(3) = F3T(LVL) FD(1) = G1T(LVL) FD(5) = G3T(LVL) END IF FV(5) = F3T(LVL) GO TO 50 END IF SUM = SUM + PSUM(LVL) LVL = LVL - 1 IF (LVL .GT. 0) GO TO 80 CONTRB = SUM LVL = 1 DLG = ERROR 90 PROB1 = PROB1 + CONTRB BOUND = BOUND + DLG EXCESS = EL - DLG EL = EL1 IF (EXCESS .GT. ZERO) EL = EL1 + EXCESS IF ((FUNC0 .GT. ZERO .AND. FUNC2 .LE. FUNC0) .OR. * (FUNC0 .LT. ZERO .AND. FUNC2 .GE. FUNC0)) THEN ZZ = -SQRT2 * Z5 PART3 = ABS(FUNC2) * ALNORM(ZZ, .FALSE.) + BOUND / ERFAC IF (PART3 .LE. EPS .OR. NTM .GE. NMAX .OR. Z5 .GE. ZU) GOTO 100 END IF Z = Z5 PART0 = PART2 FUNC0 = FUNC2 IF (Z .LT. ZU .AND. NTM .LT. NMAX) GO TO 40 100 PROB = (PROB1 - ADDN * HALF) * PROB BOUND = PART3 IF (NTM .GE. NMAX .AND. IFAULT .EQ. 0) IFAULT = 5 IF (BOUND .GT. EPS .AND. IFAULT .EQ. 0) IFAULT = 8 RETURN END SUBROUTINE PFUNC(Z, A, B, BPD, N, INF, ADDN, DERIV, FUNCN, NTM, * IERC, RESULT) C C ALGORITHM AS 251.2 APPL.STATIST. (1989), VOL.38, NO.3 C C C COMPUTE FUNCTION IN INTEGRAND AND ITS 4TH DERIVATIVE. C INTEGER NN PARAMETER (NN = 50) REAL A(*), B(*), BPD(*), FOU(NN), FOU1(4, NN), TMP(4), GOU(NN), * GOU1(4, NN), FF(4), GF(4), TERM(4), GERM(4) INTEGER INF(*) REAL ZERO, ONE, TWO, THREE, FOUR, SIX, EIGHT, TWELVE, SIXTN, * SMALL, Z, U, U1, U2, BI, HI, HLI, BP, ADDN, DERIV, FUNCN, * RESULT, RSLT1, RSLT2, DEN, SQRT2, SQRTPI, PHI, PHI1, PHI2, * PHI3, PHI4, FRM, GRM INTEGER N, NTM, IERC, INFI, I, J, K, M, L, IK REAL ALNORM EXTERNAL ALNORM DATA ZERO, ONE, TWO, THREE, FOUR, SIX, EIGHT, TWELVE, SIXTN, * SMALL /0.0, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 12.0, 16.0, 0.1E-12/ DATA SQRT2, SQRTPI /1.41421356237310, 1.77245385090552/ DERIV = ZERO NTM = NTM + 1 RSLT1 = ONE RSLT2 = ONE BI = ONE HI = A(1) + ONE HLI = B(1) + ONE INFI = -1 DO 60 I = 1, N IF (BPD(I) .EQ. BI .AND. A(I) .EQ. HI .AND. B(I) .EQ. HLI .AND. * INF(I) .EQ. INFI) THEN FOU(I) = FOU(I - 1) GOU(I) = GOU(I - 1) DO 10 IK = 1, 4 FOU1(IK, I) = FOU1(IK, I - 1) GOU1(IK, I) = GOU1(IK, I - 1) 10 CONTINUE ELSE BI = BPD(I) HI = A(I) HLI = B(I) INFI = INF(I) IF (BI .EQ. ZERO) THEN IF (INFI .LT. 1) THEN FOU(I) = ONE - ALNORM(HLI, .FALSE.) ELSE IF (INFI .EQ. 1) THEN FOU(I) = ALNORM(HI, .FALSE.) ELSE FOU(I) = ALNORM(HI, .FALSE.) - ALNORM(HLI, .FALSE.) END IF GOU(I) = FOU(I) DO 20 IK = 1, 4 FOU1(IK, I) = ZERO GOU1(IK, I) = ZERO 20 CONTINUE ELSE DEN = SQRT(ONE - BI * BI) BP = BI * SQRT2 / DEN IF (INFI .LT. 1) THEN U = -HLI / DEN + Z * BP FOU(I) = ALNORM(U, .FALSE.) CALL ASSIGN (U, BP, FOU1(1, I)) BP = -BP U = -HLI / DEN + Z * BP GOU(I) = ALNORM(U, .FALSE.) CALL ASSIGN (U, BP, GOU1(1, I)) ELSE IF (INFI .EQ. 1) THEN U = HI / DEN + Z * BP GOU(I) = ALNORM(U, .FALSE.) CALL ASSIGN (U, BP, GOU1(1, I)) BP = -BP U = HI / DEN + Z * BP FOU(I) = ALNORM(U, .FALSE.) CALL ASSIGN (U, BP, FOU1(1, I)) ELSE U2 = -HLI / DEN + Z * BP CALL ASSIGN (U2, BP, FOU1(1, I)) BP = -BP U1 = HI / DEN + Z * BP CALL ASSIGN (U1, BP, TMP(1)) FOU(I) = ALNORM(U1, .FALSE.) + ALNORM(U2, .FALSE.) - ONE DO 30 IK = 1, 4 FOU1(IK, I) = FOU1(IK, I) + TMP(IK) 30 CONTINUE IF (-HLI .EQ. HI) THEN GOU(I) = FOU(I) DO 40 IK = 1, 4 GOU1(IK, I) = FOU1(IK, I) 40 CONTINUE ELSE U2 = -HLI / DEN + Z * BP CALL ASSIGN (U2, BP, GOU1(1, I)) BP = -BP U1 = HI / DEN + Z * BP GOU(I) = ALNORM(U1, .FALSE.) + ALNORM(U2, .FALSE.)-ONE CALL ASSIGN (U1, BP, TMP(1)) DO 50 IK = 1, 4 GOU1(IK, I) = GOU1(IK, I) + TMP(IK) 50 CONTINUE END IF END IF END IF END IF RSLT1 = RSLT1 * FOU(I) RSLT2 = RSLT2 * GOU(I) IF (RSLT1 .LE. SMALL) RSLT1 = ZERO IF (RSLT2 .LE. SMALL) RSLT2 = ZERO 60 CONTINUE FUNCN = RSLT1 + RSLT2 + ADDN RESULT = FUNCN * EXP(-Z * Z) / SQRTPI C C IF 4TH DERIVATIVE IS NOT WANTED, STOP HERE. C OTHERWISE, PROCEED TO COMPUTE 4TH DERIVATIVE. C IF (IERC .EQ. 0) RETURN DO 70 IK = 1, 4 FF(IK) = ZERO GF(IK) = ZERO 70 CONTINUE DO 100 I = 1, N FRM = ONE GRM = ONE DO 80 J = 1, N IF (J .EQ. 1) GO TO 80 FRM = FRM * FOU(J) GRM = GRM * GOU(J) IF (FRM .LE. SMALL) FRM = ZERO IF (GRM .LE. SMALL) GRM = ZERO 80 CONTINUE DO 90 IK = 1, 4 FF(IK) = FF(IK) + FRM * FOU1(IK, I) GF(IK) = GF(IK) + GRM * GOU1(IK, I) 90 CONTINUE 100 CONTINUE IF (N .LE. 2) GO TO 230 DO 130 I = 1, N DO 120 J = I + 1, N TERM(2) = FOU1(1, I) * FOU1(1, J) GERM(2) = GOU1(1, I) * GOU1(1, J) TERM(3) = FOU1(2, I) * FOU1(1, J) GERM(3) = GOU1(2, I) * GOU1(1, J) TERM(4) = FOU1(3, I) * FOU1(1, J) GERM(4) = GOU1(3, I) * GOU1(1, J) TERM(1) = FOU1(2, I) * FOU1(2, J) GERM(1) = GOU1(2, I) * GOU1(2, J) DO 110 K = 1, N IF (K .EQ. I .OR. K .EQ. J) GO TO 110 CALL TOOSML (1, TERM, FOU(K)) CALL TOOSML (1, GERM, GOU(K)) 110 CONTINUE FF(2) = FF(2) + TWO * TERM(2) FF(3) = FF(3) + TWO * TERM(3) * THREE FF(4) = FF(4) + TWO * (TERM(4) * FOUR + TERM(1) * THREE) GF(2) = GF(2) + TWO * GERM(2) GF(3) = GF(3) + TWO * GERM(3) * THREE GF(4) = GF(4) + TWO * (GERM(4) * FOUR + GERM(1) * THREE) 120 CONTINUE 130 CONTINUE DO 170 I = 1, N DO 160 J = I + 1, N DO 150 K = J + 1, N TERM(3) = FOU1(1, I) * FOU1(1, J) * FOU1(1, K) TERM(4) = FOU1(2, I) * FOU1(1, J) * FOU1(1, K) GERM(3) = GOU1(1, I) * GOU1(1, J) * GOU1(1, K) GERM(4) = GOU1(2, I) * GOU1(1, J) * GOU1(1, K) IF (N .GT. 3) THEN DO 140 M = 1, N IF (M .EQ. I .OR. M .EQ. J .OR. M .EQ. K) GO TO 140 CALL TOOSML (3, TERM, FOU(M)) CALL TOOSML (3, GERM, GOU(M)) 140 CONTINUE END IF FF(3) = FF(3) + SIX * TERM(3) FF(4) = FF(4) + SIX * TERM(4) * SIX GF(3) = GF(3) + SIX * GERM(3) GF(4) = GF(4) + SIX * GERM(4) * SIX 150 CONTINUE 160 CONTINUE 170 CONTINUE IF (N .LE. 3) GO TO 230 DO 220 I = 1, N DO 210 J = I + 1, N DO 200 K = J + 1, N DO 190 M = K + 1, N TERM(4) = FOU1(1, I) * FOU1(1, J) * FOU1(1, K) * FOU1(1, M) GERM(4) = GOU1(1, I) * GOU1(1, J) * GOU1(1, K) * GOU1(1, M) IF (N .GT. 4) THEN DO 180 L = 1, N IF (L .EQ. I .OR. L .EQ. J .OR. L .EQ. K .OR. L .EQ. M)GOTO 180 CALL TOOSML (4, TERM, FOU(L)) CALL TOOSML (4, GERM, GOU(L)) 180 CONTINUE END IF FF(4) = FF(4) + FOUR * SIX * TERM(4) GF(4) = GF(4) + FOUR * SIX * GERM(4) 190 CONTINUE 200 CONTINUE 210 CONTINUE 220 CONTINUE C 230 CONTINUE PHI = EXP(-Z * Z) / SQRTPI PHI1 = -TWO * Z * PHI PHI2 = (FOUR * Z ** 2 - TWO) * PHI PHI3 = (-EIGHT * Z ** 3 + TWELVE * Z) * PHI PHI4 = (SIXTN * Z ** 2 * (Z ** 2 - THREE) + TWELVE) * PHI DERIV = PHI * (FF(4) + GF(4)) + FOUR * PHI1 * (FF(3) + GF(3)) * + SIX * PHI2 * (FF(2) + GF(2)) + FOUR * PHI3 * (FF(1) + GF(1)) * + PHI4 * FUNCN RETURN END SUBROUTINE ASSIGN (U, BP, FF) C C ALGORITHM AS 251.3 APPL.STATIST. (1989), VOL.38, NO.3 C C C COMPUTE DERIVATIVES OF NORMAL CDF'S. C REAL FF(4) REAL U, U2, BP, HALF, ONE, THREE, SQ2PI, T1, T2, T3 INTEGER I DATA HALF, ONE, THREE, SQ2PI /0.5, 1.0, 3.0, 2.50662827463100/ DATA ZERO, UMAX, SMALL /0.0, 8.0, 0.1E-07/ IF (ABS(U) .GT. UMAX) THEN DO 10 I = 1, 4 FF(I) = ZERO 10 CONTINUE ELSE U2 = U * U T1 = BP * EXP(-HALF * U2) / SQ2PI T2 = BP * T1 T3 = BP * T2 FF(1) = T1 FF(2) = -U * T2 FF(3) = (U2 - ONE) * T3 FF(4) = (THREE - U2) * U * BP * T3 DO 20 I = 1, 4 IF(ABS(FF(I)) .LT. SMALL) FF(I) = ZERO 20 CONTINUE END IF RETURN END SUBROUTINE WMAX(W1, W2, W3, DLG) C C ALGORITHM AS 251.4 APPL.STATIST. (1989), VOL.38, NO.3 C C C LARGEST ABSOLUTE VALUE OF QUADRATIC FUNCTION FITTED C TO THREE POINTS. C REAL W1, W2, W3, DLG, QUAD, QLIM, QMIN, ONE, TWO, B2C DATA ONE, TWO, QMIN /1.0, 2.0, 0.00001/ DLG = MAX( ABS(W1), ABS(W3) ) QUAD = W1 - W2 * TWO + W3 QLIM = MAX( ABS(W1 - W3) / TWO , QMIN) IF (ABS(QUAD) .LE. QLIM) RETURN B2C = (W1 - W3) / QUAD / TWO IF (ABS(B2C) .GE. ONE) RETURN DLG = MAX( DLG, ABS(W2 - B2C * QUAD * B2C / TWO) ) RETURN END SUBROUTINE TOOSML (N, FF, F) C C ALGORITHM AS 251.5 APPL.STATIST. (1989), VOL.38, NO.3 C C C MULTIPLY FF(I) BY F FOR I = N TO 4. SET TO ZERO IF TOO SMALL. C REAL FF(4), F, ZERO, SMALL INTEGER N, I DATA ZERO, SMALL /0.0, 0.1E-12/ DO 10 I = N, 4 FF(I) = FF(I) * F IF (ABS(FF(I)) .LE. SMALL) FF(I) = ZERO 10 CONTINUE RETURN END REAL FUNCTION ALNORM(X, UPPER) C C ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, P.424 C C EVALUATES THE TAIL AREA OF THE STANDARDIZED NORMAL CURVE C FROM X TO INFINITY IF UPPER IS .TRUE. OR C FROM MINUS INFINITY TO X IF UPPER IS .FALSE. C REAL LTONE, UTZERO, ZERO, HALF, ONE, CON, A1, A2, A3, $ A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7, B8, B9, $ B10, B11, B12, X, Y, Z, ZEXP LOGICAL UPPER, UP C C LTONE AND UTZERO MUST BE SET TO SUIT THE PARTICULAR COMPUTER C (SEE INTRODUCTORY TEXT) C DATA LTONE, UTZERO /7.0, 18.66/ DATA ZERO, HALF, ONE, CON /0.0, 0.5, 1.0, 1.28/ DATA A1, A2, A3, $ A4, A5, A6, $ A7 $ /0.398942280444, 0.399903438504, 5.75885480458, $ 29.8213557808, 2.62433121679, 48.6959930692, $ 5.92885724438/ DATA B1, B2, B3, $ B4, B5, B6, $ B7, B8, B9, $ B10, B11, B12 $ /0.398942280385, 3.8052E-8, 1.00000615302, $ 3.98064794E-4, 1.98615381364, 0.151679116635, $ 5.29330324926, 4.8385912808, 15.1508972451, $ 0.742380924027, 30.789933034, 3.99019417011/ C ZEXP(Z) = EXP(Z) C UP = UPPER Z = X IF (Z .GE. ZERO) GOTO 10 UP = .NOT. UP Z = -Z 10 IF (Z .LE. LTONE .OR. UP .AND. Z .LE. UTZERO) GOTO 20 ALNORM = ZERO GOTO 40 20 Y = HALF * Z * Z IF (Z .GT. CON) GOTO 30 C ALNORM = HALF - Z * (A1 - A2 * Y / (Y + A3 - A4 / (Y + A5 + $ A6 / (Y + A7)))) GOTO 40 C 30 ALNORM = B1 * ZEXP(-Y) / (Z - B2 + B3 / (Z + B4 + B5 / (Z - $ B6 + B7 / (Z + B8 - B9 / (Z + B10 + B11 / (Z + B12)))))) C 40 IF (.NOT. UP) ALNORM = ONE - ALNORM RETURN END REAL FUNCTION PPND7 (P, IFAULT) C C ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3 C C PRODUCES THE NORMAL DEVIATE Z CORRESPONDING TO A GIVEN LOWER C TAIL AREA OF P; Z IS ACCURATE TO ABOUT 1 PART IN 10**7. C C THE HASH SUMS BELOW ARE THE SUMS OF THE MANTISSAS OF THE C COEFFICIENTS. THEY ARE INCLUDED FOR USE IN CHECKING C TRANSCRIPTION. C INTEGER IFAULT REAL ZERO, ONE, HALF, SPLIT1, SPLIT2, CONST1, CONST2, * A0, A1, A2, A3, B1, B2, B3, C0, C1, C2, C3, D1, D2, * E0, E1, E2, E3, F1, F2, P, Q, R PARAMETER (ZERO = 0.0E0, ONE = 1.0E0, HALF = 0.5E0, * SPLIT1 = 0.425E0, SPLIT2 = 5.0E0, * CONST1 = 0.180625E0, CONST2 = 1.6E0) C C COEFFICIENTS FOR P CLOSE TO 1/2 PARAMETER (A0 = 3.38713 27179E0, * A1 = 5.04342 71938E1, * A2 = 1.59291 13202E2, * A3 = 5.91093 74720E1, * B1 = 1.78951 69469E1, * B2 = 7.87577 57664E1, * B3 = 6.71875 63600E1) C HASH SUM AB 32.31845 77772 C C COEFFICIENTS FOR P NEITHER CLOSE TO 1/2 NOR 0 OR 1 PARAMETER (C0 = 1.42343 72777E0, * C1 = 2.75681 53900E0, * C2 = 1.30672 84816E0, * C3 = 1.70238 21103E-1, * D1 = 7.37001 64250E-1, * D2 = 1.20211 32975E-1) C HASH SUM CD 15.76149 29821 C C COEFFICIENTS FOR P NEAR 0 OR 1 PARAMETER (E0 = 6.65790 51150E0, * E1 = 3.08122 63860E0, * E2 = 4.28682 94337E-1, * E3 = 1.73372 03997E-2, * F1 = 2.41978 94225E-1, * F2 = 1.22582 02635E-2) C HASH SUM EF 19.40529 10204 C IFAULT = 0 Q = P - HALF IF (ABS(Q) .LE. SPLIT1) THEN R = CONST1 - Q * Q PPND7 = Q * (((A3 * R + A2) * R + A1) * R + A0) / * (((B3 * R + B2) * R + B1) * R + ONE) RETURN ELSE IF (Q .LT. 0) THEN R = P ELSE R = ONE - P ENDIF IF (R .LE. ZERO) THEN IFAULT = 1 PPND7 = ZERO RETURN ENDIF R = SQRT(-LOG(R)) IF (R .LE. SPLIT2) THEN R = R - CONST2 PPND7 = (((C3 * R + C2) * R + C1) * R + C0) / * ((D2 * R + D1) * R + ONE) ELSE R = R - SPLIT2 PPND7 = (((E3 * R + E2) * R + E1) * R + E0) / * ((F2 * R + F1) * R + ONE) ENDIF IF (Q .LT. 0) PPND7 = -PPND7 RETURN ENDIF END SUBROUTINE MVSTUD(NDF,A,B,BPD,ERRB,N,INF,D,IERC,HNC,PROB, * BND,IFLT) C C COMPUTE MULTIVARIATE STUDENT INTEGRAL, C USING MVNPRD (DUNNETT, APPL. STAT., 1989) C IF RHO(I,J) = BPD(I)*BPD(J). C C IF RHO(I,J) HAS GENERAL STRUCTURE, USE C MULNOR (SCHERVISH, APPL. STAT., 1984) AND REPLACE C CALL MVNPRD(A,B,BPD,EPS,N,INF,IERC,HNC,PROB,BND,IFLT) C BY CALL MULNOR(A,B,SIG,EPS,N,INF,PROB,BND,IFLT). C C AUTHOR: C.W. DUNNETT, MCMASTER UNVERSITY C C BASED ON ADAPTIVE SIMPSON'S RULE ALGORITHM C DESCRIBED IN SHAMPINE & ALLEN: "NUMERICAL C COMPUTING", (1974), PAGE 240. C C PARAMETERS ARE SAME AS IN ALGORITHM AS 251 C IN APPL. STAT. (1989), VOL. 38: 564-579 C WITH THE FOLLOWING ADDITIONS: C NDF INTEGER INPUT DEGREES OF FREEDOM C D REAL ARRAY INPUT NON-CENTRALITY VECTOR C (PUT NDF = 0 FOR INFINITE D.F.) C INTEGER NN PARAMETER (NN=50) DIMENSION A(*),B(*),BPD(*),INF(*),D(*),F(3),AA(NN),BB(NN) DATA ZERO,HALF,TWO,THREE,FOUR / 0.0, 0.5, 2.0, 3.0, 4.0 / DO 10 I = 1, N AA(I) = A(I) - D(I) BB(I) = B(I) - D(I) 10 CONTINUE IF (NDF .LE. 0) THEN CALL MVNPRD(AA,BB,BPD,ERRB,N,INF,IERC,HNC,PROB,BND,IFLT) RETURN ENDIF BND = ZERO IFLT = 0 MAXDF = 150 ERB2 = ERRB C C CHECK IF D.F. EXCEED MAXDF; IF YES, THEN PROB C IS COMPUTED BY QUADRATIC INTERPOLATION ON 1./D.F. C IF (NDF .LE. MAXDF) GO TO 20 CALL MVNPRD(AA,BB,BPD,ERB2,N,INF,IERC,HNC,F(1),BND,IFLT) NF = MAXDF / 2 CALL SIMPSN(NF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,F(3),BND,IFLT) NF = NF * 2 CALL SIMPSN(NF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,F(2),BND,IFLT) XX = FLOAT(NF) / FLOAT(NDF) AX = F(3) - F(2)*TWO + F(1) BX = F(2)*FOUR - F(3) - F(1)*THREE PROB = F(1) + XX * (AX * XX + BX) * HALF RETURN 20 CALL SIMPSN (NDF,A,B,BPD,ERB2,N,INF,D,IERC,HNC,PROB,BND,IFLT) RETURN END SUBROUTINE SIMPSN (NDF,A,B,BPD,ERRB,N,INF,D,IERC,HNC,PROB, * BND,IFLT) C C STUDENTIZES A MULTIVARIATE INTEGRAL USING SIMPSON'S RULE. C DIMENSION A(*),B(*),BPD(*),INF(*),D(*), * FV(5),F1T(30),F2T(30),F3T(30), * LDIR(30),PSUM(30),ESTT(30),ERRR(30),GV(5),G1T(30),G2T(30), * G3T(30),GSUM(30) DATA ZERO,HALF,ONE,ONEP5,TWO,FOUR,SIX,DXMIN /0.0,0.5,1.0,1.5, * 2.0,4.0,6.0,0.000004/ PROB = ZERO BOUNDA = ZERO BOUNDG = ZERO IFLAG = 0 IER = 0 START = -ONE DAX = ONE ERB2 = ERRB * HALF EPS1 = ERB2 * HALF CALL FUN (ZERO,NDF,A,B,BPD,ERB2,N,INF,D,F0,G0,IERC,HNC,IER) 10 FV(1) = ZERO GV(1) = ZERO ERROR = ZERO DA = DAX LVL = 1 Z3 = START + HALF*DA CALL FUN(Z3,NDF,A,B,BPD,ERB2,N,INF,D,FV(3),GV(3),IERC,HNC,IER) FV(5) = F0 GV(5) = G0 WT = ABS(DA) / SIX CONTRB = WT * (FV(1) + FOUR * FV(3) + FV(5)) CONTRG = WT * (GV(1) + FOUR * GV(3) + GV(5)) LDIR(LVL) = 2 PSUM(LVL) = ZERO GSUM(LVL) = ZERO C C BISECT INTERVAL; COMPUTE ESTIMATES FOR EACH HALF. C 20 DX = HALF * DA WT = ABS(DX) / SIX Z2 = START + HALF * DX CALL FUN(Z2,NDF,A,B,BPD,ERB2,N,INF,D,FV(2),GV(2),IERC,HNC,IER) Z4 = START + ONEP5 * DX CALL FUN(Z4,NDF,A,B,BPD,ERB2,N,INF,D,FV(4),GV(4),IERC,HNC,IER) ESTL = WT * (FV(1) + FOUR * FV(2) + FV(3)) ESTR = WT * (FV(3) + FOUR * FV(4) + FV(5)) ESTGL = WT * (GV(1) + FOUR * GV(2) + GV(3)) ESTGR = WT * (GV(3) + FOUR * GV(4) + GV(5)) SUM = ESTL + ESTR SUMG = ESTGL + ESTGR DLG = ABS(CONTRB - SUM) ERRL = DLG C C STOP BISECTING WHEN ACCURACY SUFFICIENT, OR IF C INTERVAL TOO NARROW OR BISECTED TOO OFTEN. C 30 IF (DLG .LE. EPS1) GO TO 50 IF (ABS(DX) .LE. DXMIN .OR. LVL .GE. 30) GO TO 40 C C RAISE LEVEL. STORE INFORMATION FOR RIGHT HALF C AND APPLY SIMPSON'S RULE TO LEFT HALF. C LVL = LVL + 1 LDIR(LVL) = 1 F1T(LVL) = FV(3) F2T(LVL) = FV(4) F3T(LVL) = FV(5) G1T(LVL) = GV(3) G2T(LVL) = GV(4) G3T(LVL) = GV(5) DA = DX FV(5) = FV(3) FV(3) = FV(2) GV(5) = GV(3) GV(3) = GV(2) ESTT(LVL) = ESTR CONTRB = ESTL CONTRG = ESTGL EPS1 = EPS1 * HALF ERRR(LVL) = EPS1 GO TO 20 C C ACCEPT APPROXIMATE VALUE FOR INTERVAL. C 40 IFLAG = 11 50 ERROR = ERROR + ERRL 60 IF (LDIR(LVL) .EQ. 1) GO TO 70 SUM = SUM + PSUM(LVL) SUMG = SUMG + GSUM(LVL) LVL = LVL - 1 IF (LVL .GT. 0) GO TO 60 CONTRB = SUM CONTRG = SUMG LVL = 1 DLG = ERROR GO TO 80 C C RESTORE SAVED INFORMATION TO PROCESS RIGHT HALF. C 70 PSUM(LVL) = SUM GSUM(LVL) = SUMG LDIR(LVL) = 2 DA = DAX / TWO**(LVL-1) START = START + DX * TWO FV(1) = F1T(LVL) FV(3) = F2T(LVL) FV(5) = F3T(LVL) GV(1) = G1T(LVL) GV(3) = G2T(LVL) GV(5) = G3T(LVL) CONTRB = ESTT(LVL) EXCESS = EPS1 - DLG EPS1 = ERRR(LVL) IF (EXCESS .GT. ZERO) EPS1 = EPS1 + EXCESS GO TO 20 80 PROB = PROB + CONTRB BOUNDG = BOUNDG + CONTRG BOUNDA = BOUNDA + DLG IF (Z4 .LE. ZERO) GO TO 90 IF (IFLT .EQ. 0) IFLT = IER IF (IFLT .EQ. 0) IFLT = IFLAG BOUNDA = BOUNDA + BOUNDG IF (BND .LT. BOUNDA) BND = BOUNDA RETURN 90 EPS1 = ERB2 * HALF EXCESS = EPS1 - BND IF (EXCESS .GT. ZERO) EPS1 = EPS1 + EXCESS START = ONE DAX = -ONE GO TO 10 END FUNCTION SDIST(Y,N) C C COMPUTE Y**(N/2 - 1) EXP(-Y) / GAMMA(N/2) C C (Revised: 1994-01-19) C DATA ZERO, HALF, ONE, X23 / 0.0, 0.5, 1.0, -23.0 / DATA SQRTPI / 1.77245385090552 / SDIST = ZERO IF (Y .LE. ZERO) RETURN JJ = N/2 - 1 JK = 2 * JJ - N + 2 JKP = JJ - JK SDIST = ONE IF (JK .LT. 0) SDIST = SDIST / SQRT(Y) / SQRTPI IF (JKP .EQ. 0) GO TO 20 XN = FLOAT(N) * HALF TEST = ALOG(Y) - Y / FLOAT(JKP) IF ( TEST .LT. X23 ) THEN SDIST = ZERO RETURN ENDIF SDIST = ALOG ( SDIST ) DO 10 J = 1, JKP XN = XN - ONE SDIST = SDIST + TEST - ALOG(XN) 10 CONTINUE IF ( SDIST .LT. X23 ) THEN SDIST = ZERO ELSE SDIST = EXP( SDIST ) ENDIF RETURN 20 SDIST = SDIST * EXP(-Y) RETURN END SUBROUTINE FUN (Z,NDF,H,HL,BPD,ERB2,N,INF,D,F0,G0,IERC * ,HNC,IER) INTEGER NN PARAMETER (NN=50) DIMENSION A(NN),B(NN),H(*),HL(*),BPD(*),INF(*),D(*) DATA ZERO, ONE, TWO, SMALL / 0.0, 1.0, 2.0, 1.0E-08 / F0 = ZERO G0 = ZERO IF (Z .LE. -ONE .OR. Z .GE. ONE) RETURN DF = FLOAT(NDF) ARG = (ONE + Z) / (ONE - Z) TERM = ARG * DF * TWO / (ONE-Z)**2 * SDIST(DF/TWO*ARG*ARG,NDF) IF (TERM .LE. SMALL) RETURN DO 10 I = 1, N A(I) = ARG * H(I) - D(I) B(I) = ARG * HL(I) - D(I) 10 CONTINUE CALL MVNPRD (A,B,BPD,ERB2,N,INF,IERC,HNC,PROB,BND,IFLT) IF (IER .EQ. 0) IER = IFLT G0 = TERM * BND F0 = TERM * PROB RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * C Charles Dunnett C Dept. of Mathematics and Statistics C McMaster University C Hamilton, Ontario L8S 4K1 C Canada C E-mail: dunnett@mcmaster.ca C Tel.: 905-525-9140 (Ext. 23414) C * * * * * * * * * * * * * * * * *