SUBROUTINE MECDF(NDIM, D, RHO, PROB, IER) INTEGER NDIM, IER REAL D(*), PROB, RHO(*) C----------------------------------------------------------------- C 6/29/90 C This subroutine is designed to calculate the MVN CDF C using the Mendell-Elston procedure as described in C Kamakura (1989). The current version is set up to go C as high as 19 dimensions (=> 20 MNP alternatives) C NOTE: Equation (15) in Kamakura has an error. C C Specifically, assume that Z is a set of random variables C with a standard normal distribution with correlations C stored in RHO (in packed form). Then this subroutine C calculates Prob[Z(1)>D(1);...; Z(NDIM) > D(NDIM)]. C----------------------------------------------------------------- REAL ALNORM, PHI EXTERNAL ALNORM, PHI INTEGER MAXALT, NMAX PARAMETER (MAXALT=20, NMAX=MAXALT-1) INTEGER I, IM1, IR, J, JM1, K, KM1 REAL PROBI, TMP REAL R(NMAX,NMAX,0:NMAX-1), SIG(NMAX,0:NMAX-1), 1 U(NMAX), UUMZ(NMAX-1), Z(NMAX,0:NMAX-1) REAL ONE, ZERO PARAMETER (ONE=1.E0, ZERO=0.E0) C----------------------------------------------------------------- C Test dimension IER = 0 IF (NDIM.GT.NMAX) THEN IER = -1 RETURN ENDIF C Set up arrays IR = 0 DO 10 I = 1, NDIM Z(I,0) = D(I) DO 10 J = 1, I-1 IR = IR + 1 R(J,I,0) = RHO(IR) 10 CONTINUE PROB = ALNORM(Z(1,0), .TRUE.) IF (PROB.LE.ZERO) THEN IER = 1 RETURN ENDIF U(1) = PHI(Z(1,0), ZERO)/PROB UUMZ(1) = U(1)*(U(1)-Z(1,0)) C Main loop DO 40 I = 2, NDIM IM1 = I-1 DO 30 J = 1, IM1 JM1 = J-1 DO 20 K = 1, JM1 KM1 = K-1 TMP = R(J,I,KM1)-R(K,J,KM1)*R(K,I,KM1)*UUMZ(K) R(J,I,K) = TMP/SIG(J,K)/SIG(I,K) 20 CONTINUE SIG(I,J) = SQRT(ONE - UUMZ(J)*R(J,I,JM1)**2) Z(I,J) = (Z(I,JM1)-U(J)*R(J,I,JM1))/SIG(I,J) 30 CONTINUE PROBI = ALNORM(Z(I,IM1), .TRUE.) IF (PROBI.LE.ZERO) THEN IER = I RETURN ENDIF PROB = PROB * PROBI IF (I.LT.NDIM) THEN U(I) = PHI(Z(I,IM1), ZERO)/PROBI UUMZ(I) = U(I)*(U(I)-Z(I,IM1)) ENDIF 40 CONTINUE END C--------------------------------------------------- REAL FUNCTION PHI(X, Y) REAL X, Y REAL ARG REAL HALF, SQ2P, XLOW, ZERO PARAMETER (HALF = 0.5E0, SQ2P = 0.91893853320467274E0, 1 XLOW = -87.E0, ZERO = 0.E0) PHI = ZERO ARG = -HALF * X * X - SQ2P - Y IF (ARG .GT. XLOW) PHI = EXP(ARG) END C--------------------------------------------------- REAL FUNCTION ALNORM(X,UPPER) REAL X LOGICAL UPPER C C ALGORITHM AS 66 BY I.D. HILL C C Full reference: Hill, I.D. (1973), "Algorithm AS 66: The Normal C Integral," Applied Statistics, 22, 424-427. C LOGICAL UP REAL Y, Z REAL CON, HALF, LTONE, ONE, UTZERO, ZERO PARAMETER (CON=1.28E0, HALF=0.5E0, LTONE=5.E0, ONE=1.E0, 1 UTZERO=12.5E0, ZERO=0.E0) UP=UPPER Z=X IF(Z.GE.ZERO) GO TO 10 UP=.NOT.UP Z=-Z 10 IF(Z .LE. LTONE .OR. UP .AND. Z .LE. UTZERO) GO TO 20 ALNORM = ZERO GO TO 40 20 Y=HALF*Z*Z IF(Z.GT.CON) GO TO 30 ALNORM = HALF - Z * (0.398942280444E0 - 0.399903438504E0*Y/ 1 (Y + 5.75885480458E0 - 29.8213557808E0/ 2 (Y + 2.62433121679E0 + 48.6959930692E0/ 3 (Y + 5.92885724438E0)))) GO TO 40 30 ALNORM = 0.398942280385E0 * EXP(-Y)/ 1 (Z - 3.8052E-8 + 1.00000615302E0/ 2 (Z + 3.98064794E-4 + 1.98615381364E0/ 3 (Z - 0.151679116635E0 + 5.29330324926E0/ 4 (Z + 4.8385912808E0 - 15.1508972451E0/ 5 (Z + 0.742380924027E0 + 30.789933034E0/ 6 (Z + 3.99019417011E0)))))) 40 IF(.NOT.UP) ALNORM = ONE - ALNORM END