! Last change: DOS 29 Jul 2000 11:50 am ! *** copyright 2000 *** ! *** filename bvnrml.f95 *** John F. Monahan ** ! ********************** PROGRAM BVNRML ! TESTS OF P2NORM, PDIVGI FOR COMPUTING BIVARIATE NORMAL PROBS ! P2NORM USES DRESNER'S METHOD ! PDIVGI USES DIVGI'S METHOD IMPLICIT NONE REAL X,Y,R,P2,P3,DP REAL P2NORM,PDIVGI INTEGER I,J,L ! 21 FORMAT(3F6.2,3F12.8) 25 FORMAT(' X Y R DRESNER DIVGI DIFF') ! OPEN( UNIT=6, FILE='bvnrml.out') ! WRITE(6,25) ! LOTS OF DIFFERENT VALUES DO I = 1,6 DO J = 1,6 DO L = 1,7 X = ( 2.*I - 7. )/2. Y = (2.*J - 7.)/2. R = 0.3*(L-4) P2 = P2NORM(X,Y,R) P3 = PDIVGI(X,Y,R) DP = P2-P3 WRITE(6,21) X,Y,R,P2,P3,DP END DO ! LOOP ON L END DO ! LOOP ON J END DO ! LOOP ON I STOP END PROGRAM BVNRML REAL FUNCTION P2NORM(X,Y,RHO) ! COMPUTE THE BIARIATE NORMAL PROBABILITY ! P2NORM(X,Y,RHO) = PR( Z1 < X, Z2 < Y ) ! WHEN Z1, Z2 ARE BIVARIATE NORMAL WITH ZERO MEANS, UNIT VARIANCE ! AND CORRELATION RHO ! ! METHOD IS DRESNER'S THIS ROUTINE CALLS Q2NORM, PHI (NORMAL DF) AND ! PHIHKR (CORE ROUTINE) TO COVER ALL THE VARIOUS CASES ! Z. DRESNER, COMPUTATION OF THE BIVARIATE NORMAL INTEGRAL, ! MATHEMATICS OF COMPUTATION V 32, PP277-279 ! IMPLICIT NONE REAL, INTENT(IN) :: X,Y,RHO ! REAL RHONWH,RHONWK,DELTA REAL CDFN ! TAKE CARE OF SPECIAL CASE OF INDEPENDENCE IF ( RHO .EQ. 0. ) THEN P2NORM = CDFN(X)*CDFN(Y) RETURN END IF ! ( RHO .EQ. 0. ) ! GET TO CASE WHERE X*Y*RHO LE 0 IF( X*Y*RHO .LE. 0. ) THEN ! Q2NORM HANDLES THIS CASE P2NORM = Q2NORM(X,Y,RHO) RETURN END IF ! ( X*Y*RHO .LE. 0. ) ! USE FORMULA TO GET TO X*Y*RHO LE 0 RHONWH = SIGN(1.,X)*(RHO*X-Y)/SQRT(X*X - 2.*RHO*X*Y + Y*Y ) RHONWK = SIGN(1.,Y)*(RHO*Y-X)/SQRT(Y*Y - 2.*RHO*Y*X + X*X ) DELTA = ( 1. - SIGN(1.,X)*SIGN(1.,Y) )/4. ! THIS NOW CONVERTS TO TWO PIECES -- DRESNER'S (11) P2NORM = Q2NORM(X,0.,RHONWH) + Q2NORM(Y,0.,RHONWK) - DELTA RETURN ! *** Q2NORM AND PHIHKR ARE INTERNAL SUBPROGRAMS ! TO PDIVGI ONLY FOR CONVENIENCE -- THEY COULD ! STAND ALONE, BUT ARE NOT USED ELSEWHERE *** CONTAINS ! REAL FUNCTION Q2NORM(X,Y,RHO) ! COMPUTES THE BIVARIATE NORMAL PROBABILITIES ! HANDLES THE CASE WHERE X*Y*RHO LE 0. ! CALLED BY P2NORM, CALLS PHI (NORMAL DF) AND PHIHKR (CORE ROUTINE) ! METHOD IS DRESNER'S -- THIS ROUTINE CONVERTS FORMULAS ! Z. DRESNER, COMPUTATION OF THE BIVARIATE NORMAL INTEGRAL, ! MATHEMATICS OF COMPUTATION V 32, PP277-279 ! IMPLICIT NONE REAL, INTENT(IN) :: X,Y,RHO REAL CDFN,PHIHKR ! SPLIT UP THE DIFFERENT CASES IF( RHO .EQ. 0. ) THEN ! HANDLE CASE OF INDEPENDENCE Q2NORM = CDFN(X)*CDFN(Y) RETURN END IF ! ( RHO .EQ. 0. ) IF( RHO .LT. 0. ) THEN ! RHO NEGATIVE HERE IF( X .LE. 0. ) THEN ! COMPUTE DIRECTLY -- EACH ONE NEGATIVE Q2NORM = PHIHKR(X,Y,RHO) RETURN ELSE ! BOTH X AND Y POS HERE Q2NORM = CDFN(X) + CDFN(Y) - 1. + PHIHKR(-X,-Y,RHO) RETURN END IF ! ( X .LE. 0. ) ! ELSE ! NOW RHO IS POSITIVE IF( X .GT. 0. ) THEN ! (FORMULA 8) X POS, Y NEG, RHO POS Q2NORM = CDFN(Y) - PHIHKR(-X,Y,-RHO) RETURN ELSE ! X NEG, Y POS, RHO POS Q2NORM = CDFN(X) - PHIHKR(X,-Y,-RHO) RETURN END IF ! ( X .GT. 0. ) END IF ! ( RHO .LE. 0. ) END FUNCTION Q2NORM END FUNCTION P2NORM ! REAL FUNCTION PHIHKR(H,K,R) ! DRESNER'S METHOD TO COMPUTE BIVARIATE NORMAL PROBABILITIES ! THIS ROUTINE COMPUTES TAIL FOR H,K,R EACH NONPOSITIVE ! MAIN CALL IS P2NORM, WHICH CALLS Q2NORM, BOTH CALL THIS ONE ! P2NORM, Q2NORM BOTH USE FORMULAS TO GET TO H,K,R EACH LE 0 FOR THIS ! Z. DRESNER, COMPUTATION OF THE BIVARIATE NORMAL INTEGRAL, ! MATHEMATICS OF COMPUTATION V 32, PP277-279 ! IMPLICIT NONE REAL, INTENT(IN) :: H,K,R REAL, PARAMETER :: PI = 3.141593 REAL, PARAMETER :: SRHAF = .7071068 ! 1/SQRT(2) REAL SROMR2,H1,K1,XI,XJ,F,S REAL, DIMENSION(5) :: ABSC,WGHTS INTEGER I,J ! ! ABSCISSAS AND WEIGHTS FOR INTEGRATING EXP(-X*X)F(X)DX FROM 0 TO INF ! N.M.STEEN, G.D.BYRNE, AND E.M.GELBARD (1969) GAUSSIAN QUADRATURES ! MATHEMATICS OF COMPUTATION V 23, PP661-671 ! ACCURACY CAN BE IMPROVED JUST BY ADDING MORE POINTS ! DATA ABSC/ .1002422, .482814, 1.06095, 1.779729, 2.66976/ DATA WGHTS/.2484062, .3923311, .2114182, .3324666E-1, .8248533E-3/ ! SROMR2 = SQRT( DIM(1.,R*R) ) H1 = SRHAF*H/SROMR2 K1 = SRHAF*K/SROMR2 ! S = 0. ! WE ARE USING 5 POINTS IN EACH DIMENSION OF INTEGRATION DO I = 1,5 XI = ABSC(I) DO J = 1,5 XJ = ABSC(J) F = H1*(2.*XI-H1) + K1*(2.*XJ-K1) + 2.*R*(XI-H1)*(XJ-K1) ! IF F IS TOO SMALL AVOID UNDERFLOW AND COMPUTATION IF( F .GT. - 140. ) THEN F = EXP( F ) S = S + WGHTS(I)*WGHTS(J)*F END IF ! ( F .GT. -140. ) END DO ! LOOP ON J END DO ! LOOP ON I PHIHKR = SROMR2*S/PI RETURN END FUNCTION PHIHKR REAL FUNCTION PDIVGI(H,K,RHO) ! COMPUTES THE BIVARIATE NORMAL PROBABILITY USING DIVGI'S METHOD ! PDIVGI(H,K,RHO) = PR( Z1 < H, Z2 < K ) ! WHERE Z1, Z2 ARE BIVARIATE NORMAL WITH ZERO MEANS, CORRELATION RHO ! AND UNIT VARIANCES ! ! METHOD IS DIVGI'S, THIS ROUTINE CALLS RUBENW, TO TAKE CARE OF ! EXTENDING THE RANGE THROUGH DUPLICATION FORMULAS, RUBENW CALLS ! DIVGIA, WHICH USES DIVGI'S APPROXIMATION FOR THE W FUNCTION. ! ! D. R. DIVGI, "CALCULATION OF THE UNIVARIATE AND BIVARIATE NORMAL ! PROBABILITY FUNCTIONS," ANNALS OF STATISTICS (1979) PP903-910. ! REAL, INTENT(IN) :: H,K,RHO REAL RP,SROMR2,PO2MT,TMPHI RP = (H*H + K*K -2.*H*K*RHO) SROMR2 = SQRT(1. - RHO*RHO) RP = SQRT(RP) / SROMR2 ! REMEMBER DIVGI COMPUTES L FOR UPPER TAIL PROBABILITIES ! SO REVERSE SIGNS OF H AND K (RHO IS OK) PO2MT = ATAN2( -H*SROMR2, RHO*H-K) TMPHI = ATAN2( -K*SROMR2, RHO*K-H) PDIVGI = RUBENW(RP,PO2MT) + RUBENW(RP,TMPHI) IF( (H .GT. 0.) .AND. (K .GT. 0.) ) PDIVGI = PDIVGI + 1. ! RETURN ! *** RUBENW AND DIVGIA ARE INTERNAL SUBPROGRAMS ! TO PDIVGI ONLY FOR CONVENIENCE -- THEY COULD ! STAND ALONE, BUT ARE NOT USED ELSEWHERE *** CONTAINS REAL FUNCTION RUBENW(R,P) ! USES DUPLICATION FORMULAS AS PART OF DIVGI'S EMTHOD FOR ! COMPUTING BIAVARIATE NORMAL PROBABILITIES ! REAL, INTENT(IN) :: R,P REAL, PARAMETER :: PI = 3.141593 REAL, PARAMETER :: PIO2 = 1.570796 ! PI/2 REAL AP, RSP AP = ABS(P) IF( AP .GE. PIO2 ) THEN AP = PI - AP RSP = R*SIN(AP) RUBENW = CDFN(-RSP) - DIVGIA(R,AP) ELSE RUBENW = DIVGIA(R,AP) END IF ! ( AP .GE. PIO2 ) IF( P .LT. 0.) RUBENW = - RUBENW RETURN END FUNCTION RUBENW REAL FUNCTION DIVGIA(R,PSI) ! COMPUTES DIVGI'S APPROXIMATION FOR RUBEN'S W FUNCTION AS ! PART OF DIVGI'S METHOD FOR COMPUTING BIVARIATE NORMAL PROBABILITIES ! ! D. R. DIVGI "CALCULATION OF UNIVARIATE AND BIVARIATE NORMAL ! PROABABILITY FUCNTIONS," ANNALS OF STATISTICS, (1979) PP903-910 ! ! REALLY A N APPROXIMATION OF RUBEN'S W FUNCTION ! IMPLICIT NONE REAL, INTENT(IN) :: R, PSI REAL, PARAMETER :: TWOPI = 6.2831853 REAL, DIMENSION(11) :: D,CSINTP REAL CP,SP,CPKM1,BR,RK INTEGER K DATA D/ 1.253298, -.9997317, .625019, -.3281916, .1470332, & & -.5494856E-1, .1629828E-1, -.3591258E-2, .540662E-3, & & -.4890254E-4, .1984741E-5 / ! ! TAKE CARE OF R TOO LARGE DIVGIA = 0. IF( R .GT. 17. ) RETURN ! ! CREATE TABLE OF INTEGRALS OF COS**(K+1) FOR LATER ! CP = COS(PSI) SP = SIN(PSI) CPKM1 = CP CSINTP(1) = SP CSINTP(2) = (SP*CP + PSI)/2. ! USE REDUCTION FORMULA FOR INTEGRAL OF COS**K DO K = 3,11 CPKM1 = CPKM1*CP CSINTP(K) = (SP*CPKM1 + (K-1)*CSINTP(K-2) )/K END DO ! LOOP ON K RK = 1. ! BR IS FORMULA IN BRACKETS BR = PSI DO K = 1,11 RK = RK * R BR = BR - D(K)*RK*CSINTP(K) END DO ! LOOP ON K DIVGIA = EXP(-R*R/2.)*BR/TWOPI RETURN END FUNCTION DIVGIA END FUNCTION PDIVGI ! REAL FUNCTION CDFN(X) ! CDF FOR NORMAL DIST 22 MAY 74 JFM ! ** RECODED AS SEPARATE SUBPROGRAMS ** 14 JULY 1988 ** ! ** RECODED FOR READABILITY ** 25 JUNE 1996 ** ! ** RECODED FOR FORTRAN 95 ** 30 JAN 2000 ** ! ! CDFN=CDF FOR NORMAL DIST CDCN=1-CDFN ! CDFN=0.5*(1.0+ERF()/SQRT(2))) ! ! RANGE FOR X METHOD SIZE NO. ACCURACY ! (0.0,0.14) CONTINUED FRACTION 3 6.7.9 7 ! (0.14,3.54) RATIONAL FUNCTION 3,4 5663 6,MOSTLY 7 ! (3.54,11.3) RATIONAL FUNCTION 4,5 5704 5,MOSTLY 6 ! (11.3,13.15) RATIONAL FUNCTION 1,2 5721 5 ! NO. REFERS TO EQN NO. OR METHOD FROM 'COMPUTER APPROXIMATIONS' ! ED BY J F HART WILEY 1968 ! REFERENCE FOR ACCURACY 'TABLES OF ERROR FN AND ITS DERIVATIVE' ! NATL BUREAU OF STANDARDS AMS 41 1958 ! ACCURACY IS IN SIG FIGS FOR UNIVAC 1108 ! IMPLICIT NONE REAL, INTENT(IN) :: X REAL, DIMENSION(4) :: P2 REAL, DIMENSION(5) :: Q2,P3 REAL, DIMENSION(6) :: Q3 REAL, DIMENSION(2) :: P4 REAL, DIMENSION(3) :: Q4 REAL XA,XS,V DATA P2/ .100046412E+2, .842655286E+1, .346025933E+1, & & .562353612/ DATA Q2/ .100046412E+2, .197155807E+2, .157022881E+2, & & .609074879E+1, 1.0 / DATA P3/ .183983860E+2, .224003952E+2, .130613857E+2, & & .402838962E+1, .564201377/ DATA Q3/ .183983860E+2, .431607573E+2, .433645871E+2, & & .236403038E+2, .714083720E+1, 1.0/ DATA P4/ .148459208, .564187743/ DATA Q4/ .511437285, .262770594, 1.0/ ! DEFINE SOME STATEMENT FUNCTIONS ! RESCALE SINCE ORIGINAL CODING WAS FOR ERROR FUNCTION XA = ABS(X/1.414214) ! ABS( X / SQRT(2) ) XS = X*X/2. IF( XA .LE. 0.1 ) THEN ! FIRST PIECE FOR SMALLEST X -- USE CONTINUED FRACTION V = 1.12837917*XA / ( EXP(XS) * (1.0-2.0*XS/(3.0+0.8*XS)) ) CDFN = ( 1. + SIGN( V, X) )/2. RETURN END IF ! IF( XA .LE. 0.1 ) IF( XA .LE. 2.5 ) THEN ! SECOND PIECE -- USE STRAIGHT RATIONAL FUNCTION (3,4) CDFN = (P2(1)+XA*(P2(2)+XA*(P2(3)+XA*P2(4)))) / & & (2.*EXP(XS)*(Q2(1)+XA*(Q2(2)+XA*(Q2(3)+XA*(Q2(4)+XA*Q2(5)))))) IF( X .GT. 0. ) CDFN = 1. - CDFN RETURN END IF ! IF( XA .LE. 2.5 ) IF( X .GT. 5. ) THEN ! TOO FAR IN RIGHT TAIL -- RETURN ONE CDFN = 1. RETURN END IF ! IF( X .GT. 5. ) IF( XA .LE. 8. ) THEN ! THIRD PIECE -- RATIONAL FUNCTION (4,5) FOR NEAR TAIL CDFN = (P3(1)+XA*(P3(2)+XA*(P3(3)+XA*(P3(4)+XA*P3(5))))) / (2.* & & (Q3(1)+XA*(Q3(2)+XA*(Q3(3)+XA*(Q3(4)+XA*(Q3(5)+XA*Q3(6))))))) CDFN = CDFN / EXP(XS) IF( X .GT. 0. ) CDFN = 1. - CDFN RETURN END IF ! IF( XA .LE. 8. ) IF( X .GE. -9.3 ) THEN ! FOURTH (TAIL) PIECE -- ANOTHER RATIONAL FUNCTION (1,2) FOR TAIL CDFN = (P4(1)+XA*P4(2))*EXP(-XS)/(2.*(Q4(1)+XA*(Q4(2)+XA*Q4(3)))) RETURN END IF ! IF( XA .LE. 9.3 ) ! TOO FAR IN LEFT TAIL -- RETURN ZERO CDFN = 0. RETURN END FUNCTION CDFN ! *** end of filename bvnrml.f95 *********************