C C f2py -m mvnprd -h mvnprd.pyf mvnprd.f only: mvnprd C edit mvnprd.pyf with input and output and then C f2py mvnprd.pyf mvnprd.f -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 C C f2py --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 -m mvnprd -c mvnprd.f C Altarnative: compile mvnprd and link to it through mvnprd_interface.f C C gfortran -fPIC -c mvnprd.f C f2py -m mvnprdmod -c mvnprd.o mvnprd_interface.f --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 C C df -c mvnprd.f C f2py -m mvnprdmod -c mvnprd.obj mvnprd_interface.f --fcompiler=compaqv --compiler=mingw32 -lmsvcr71 ! This is a MEX-file for MATLAB. ! and contains a mex-interface to Charles W. Dunnett's programs ! ,MVNPRD and MVSTUD subroutines for computing multivariate normal ! or student T probabilities with product correlation structure. The ! file should compile without errors on (Fortran77) standard Fortran ! compilers. * * The mex-interface was written by * Per Andreas Brodtkorb * Norwegian Defence Research Establishment * P.O. Box 115m * N-3191 Horten * Norway * Email: Per.Brodtkorb@ffi.no * * 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. 27104) * * MVNPRDMEX Computes multivariate normal or student T probability * with product correlation structure. * * CALL [value,bound,inform] = mvnprdmex(RHO,A,B,D,NDF,abseps,IERC,HNC) * * RHO REAL, array of coefficients defining the correlation * coefficient by: * correlation(I,J) = RHO(I)*RHO(J) for J/=I * where * 1 < RHO(I) < 1 * A REAL, array of lower integration limits. * B REAL, array of upper integration limits. * NOTE: any values greater the 37, are considered as * infinite values. * D Real array of means * NDF Degrees of freedom, NDF<=0 gives normal probabilities * ABSEPS REAL absolute error tolerance. * IERC INTEGER 1 if strict error control based on fourth * derivative * 0 if intuitive error control based on halving the * intervals * HINC REAL start interval width of simpson rule * * OUTPUT: * VALUE REAL estimated value for the integral * BOUND REAL bound on the error of the approximation * INFORM INTEGER, termination status parameter: * 0, if normal completion with ERROR < EPS; * 1, if N > 100 or N < 1. * 2, IF any abs(rho)>=1 * 4, if ANY(B(I)<=A(i)) * 5, if number of terms computed exceeds maximum number of * evaluation points * 6, if fault accurs in normal subroutines * 7, if subintervals are too narrow or too many * 8, if bounds exceeds abseps * * * MVNPRDMEX calculates multivariate normal or student T probability * with product correlation structure for rectangular regions. * The accuracy is up to around single precision, i.e., about 1e-7. * * This file was successfully compiled for matlab 5.3 * using Compaq Visual Fortran 6.1, and Windows 2000. * The example here uses Fortran77 source. * First, you will need to modify your mexopts.bat file. * To find it, issue the command prefdir(1) from the Matlab command line, * the directory it answers with will contain your mexopts.bat file. * Open it for editing. The first section will look like: * *rem ******************************************************************** *rem General parameters *rem ******************************************************************** *set MATLAB=%MATLAB% *set DF_ROOT=C:\Program Files\Microsoft Visual Studio *set VCDir=%DF_ROOT%\VC98 *set MSDevDir=%DF_ROOT%\Common\msdev98 *set DFDir=%DF_ROOT%\DF98 *set PATH=%MSDevDir%\bin;%DFDir%\BIN;%VCDir%\BIN;%PATH% *set INCLUDE=%DFDir%\INCLUDE;%DFDir%\IMSL\INCLUDE;%INCLUDE% *set LIB=%DFDir%\LIB;%VCDir%\LIB * * then you are ready to compile this file at the matlab prompt using the * following command: * mex -O mvnprdmex.f C The rest of this file contains: C 1. A Readme file provided by Charles Dunnett, the author of AS 251. C 2. The published algorithm AS 251 together with the two other AS algorithms C which it calls (AS 66 and AS 241). C 3. A driver program (MVTIN) for either multivariate normal or t. C 4. MVSTUD for calculating multivariate t probabilities. C *************************************************************************** C Date: Mon, 10 Apr 1995 16:49:10 +0059 (EDT) C From: "Charles W. Dunnett" C Subject: Readme for AS251 (extended version incl. multivariate t) C MVTIN is a driver program for computing multivariate normal or t C probability integrals over arbitrary rectangular regions. The C correlation structure is assumed to be of product form, rho_ij = C b_i x b_j, where -1 < b_i < +1. C It requires the following:- C 1. MVNPRD (published as algorithm AS 251 in Applied C Statistics (1989), 38: 564-579; see also the correction C note in Applied Statistics (1993), 42: 709), C 2. ALNORM and PPND7 (published as algorithms AS 66 and C AS 241, respectively, in Applied Statistics, and C 3. MVSTUD ( which Studentizes MVNPRD). SUBROUTINE MVNPRD(A, B, BPD, EPS, N, INF, IERC, HINC, PROB, BOUND, * IFAULT) implicit none 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 = 100) DOUBLE PRECISION 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) DOUBLE PRECISION 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, TSUM, EXCESS, ERROR, * PROB1, SAFE, ONEP5,X2880 INTEGER N, IERC, IFAULT, I, NTM, NMAX, LVL, NR, NDIM DOUBLE PRECISION 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 NR = NINT(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)) TSUM = ESTL + ESTR DLG = ABS(CONTRB - TSUM) 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) TSUM = 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) = TSUM 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 TSUM = TSUM + PSUM(LVL) LVL = LVL - 1 IF (LVL .GT. 0) GO TO 80 CONTRB = TSUM 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) implicit none 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 = 100) DOUBLE PRECISION 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(*) DOUBLE PRECISION 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 DOUBLE PRECISION 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) implicit none C C ALGORITHM AS 251.3 APPL.STATIST. (1989), VOL.38, NO.3 C C C COMPUTE DERIVATIVES OF NORMAL CDF'S. C DOUBLE PRECISION FF(4) DOUBLE PRECISION U, U2, BP, HALF, ONE, THREE, SQ2PI, T1, T2, T3 $ ,ZERO, UMAX, SMALL 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) implicit none 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 DOUBLE PRECISION 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) implicit none 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 DOUBLE PRECISION 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 DOUBLE PRECISION FUNCTION ALNORM(X, UPPER) implicit none 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 DOUBLE PRECISION 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 DOUBLE PRECISION FUNCTION PPND7 (P, IFAULT) implicit none 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 DOUBLE PRECISION 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) implicit none 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 DOUBLE PRECISION :: HNC,PROB,BND INTEGER :: NN, MAXDF,I,IERC,NDF,N,IFLT PARAMETER (NN=100, MAXDF = 150) integer :: INF(*) DOUBLE PRECISION :: A(*),B(*),BPD(*),D(*),F(3), & AA(NN),BB(NN) DOUBLE PRECISION :: ERB2, ERRB, AX,BX,XX DOUBLE PRECISION,SAVE :: ZERO,HALF,TWO,THREE,FOUR INTEGER :: NF !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 / !external float 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 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 = DBLE(NF) / DBLE(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) implicit none C C STUDENTIZES A MULTIVARIATE INTEGRAL USING SIMPSON'S RULE. C double precision :: A,B,BPD,D, * FV,F1T,F2T,F3T, * LDIR,PSUM,ESTT,ERRR,GV,G1T,G2T, * G3T,GSUM double precision :: PROB, BOUNDA, BOUNDG,sTART, DAX, ERB2,ERRB, $ EPS1, F0,G0, HNC, ERROR,DA,Z3,WT, CONTRG,DX,Z2,Z4,ESTL,ESTR, $ ESTGL,ESTGR,CONTRB,TSUM,SUMG,DLG,ERRL,EXCESS,BND double precision :: ZERO,HALF,ONE,ONEP5,TWO,FOUR,SIX,DXMIN INTEGER :: IFLAG, IER, NDF,N, INF, IERC,LVL,IFLT 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)) TSUM = ESTL + ESTR SUMG = ESTGL + ESTGR DLG = ABS(CONTRB - TSUM) 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 TSUM = TSUM + PSUM(LVL) SUMG = SUMG + GSUM(LVL) LVL = LVL - 1 IF (LVL .GT. 0) GO TO 60 CONTRB = TSUM CONTRG = SUMG LVL = 1 DLG = ERROR GO TO 80 C C RESTORE SAVED INFORMATION TO PROCESS RIGHT HALF. C 70 PSUM(LVL) = TSUM 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 DOUBLE PRECISION FUNCTION SDIST(Y,N) implicit none C C COMPUTE Y**(N/2 - 1) EXP(-Y) / GAMMA(N/2) C C (Revised: 1994-01-19) C DOUBLE PRECISION :: Y,XN,TEST, ZERO, HALF, ONE, X23,SQRTPI INTEGER :: N,JJ,JK,JKP,J 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 = DBLE(N) * HALF TEST = LOG(Y) - Y / DBLE(JKP) IF ( TEST .LT. X23 ) THEN SDIST = ZERO RETURN ENDIF SDIST = LOG ( SDIST ) DO 10 J = 1, JKP XN = XN - ONE SDIST = SDIST + TEST - LOG(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) implicit none double precision :: ZERO, ONE, TWO, SMALL, Z, arg, term , f0, g0 $ ,df,ERB2,HNC,BND,PROB INTEGER NN,NDF,N, I, IER,IERC,IFLT PARAMETER (NN=100) DOUBLE precision :: A,B,H,HL,BPD,D, SDIST integer :: INF DIMENSION A(NN),B(NN),H(*),HL(*),BPD(*),INF(*),D(*) DATA ZERO, ONE, TWO, SMALL / 0.0, 1.0, 2.0, 1.0E-08 / external SDIST F0 = ZERO G0 = ZERO IF (Z .LE. -ONE .OR. Z .GE. ONE) RETURN DF = DBLE(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. 27104) C * * * * * * * * * * * * * * * * * * * * * * * * * * * *