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 * * * * * * * * * * * * * * * * *




