function bivnor ( ah, ak, r ) !*****************************************************************************80 ! !! BIVNOR computes the bivariate normal CDF. ! ! Discussion: ! ! BIVNOR computes the probability for two normal variates X and Y ! whose correlation is R, that AH <= X and AK <= Y. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 13 April 2012 ! ! Author: ! ! Original FORTRAN77 version by Thomas Donnelly. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Thomas Donnelly, ! Algorithm 462: Bivariate Normal Distribution, ! Communications of the ACM, ! October 1973, Volume 16, Number 10, page 638. ! ! Parameters: ! ! Input, real ( kind = 8 ) AH, AK, the lower limits of integration. ! ! Input, real ( kind = 8 ) R, the correlation between X and Y. ! ! Output, real ( kind = 8 ) BIVNOR, the bivariate normal CDF. ! ! Local Parameters: ! ! Local, integer ( kind = 4 ) IDIG, the number of significant digits ! to the right of the decimal point desired in the answer. ! implicit none real ( kind = 8 ) a2 real ( kind = 8 ) ah real ( kind = 8 ) ak real ( kind = 8 ) ap real ( kind = 8 ) b real ( kind = 8 ) bivnor real ( kind = 8 ) cn real ( kind = 8 ) con real ( kind = 8 ) conex real ( kind = 8 ) ex real ( kind = 8 ) g2 real ( kind = 8 ) gauss real ( kind = 8 ) gh real ( kind = 8 ) gk real ( kind = 8 ) gw real ( kind = 8 ) h2 real ( kind = 8 ) h4 integer ( kind = 4 ) i integer ( kind = 4 ), parameter :: idig = 15 integer ( kind = 4 ) is real ( kind = 8 ) r real ( kind = 8 ) rr real ( kind = 8 ) s1 real ( kind = 8 ) s2 real ( kind = 8 ) sgn real ( kind = 8 ) sn real ( kind = 8 ) sp real ( kind = 8 ) sqr real ( kind = 8 ) t real ( kind = 8 ), parameter :: twopi = 6.283185307179587D+00 real ( kind = 8 ) w2 real ( kind = 8 ) wh real ( kind = 8 ) wk gauss ( t ) = ( 1.0D+00 + erf ( t / sqrt ( 2.0D+00 ) ) ) / 2.0D+00 ! GAUSS is a univariate lower normal tail area calculated here from the ! central error function ERF. ! b = 0.0D+00 gh = gauss ( - ah ) / 2.0D+00 gk = gauss ( - ak ) / 2.0D+00 if ( r == 0.0D+00 ) then b = 4.0D+000 * gh * gk b = max ( b, 0.0D+00 ) b = min ( b, 1.0D+00 ) bivnor = b return end if rr = ( 1.0D+00 + r ) * ( 1.0D+00 - r ) if ( rr < 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BIVNOR - Fatal error!' write ( *, '(a)' ) ' 1 < |R|.' stop end if if ( rr == 0.0D+00 ) then if ( r < 0.0D+00 ) then if ( ah + ak < 0.0D+00 ) then b = 2.0D+00 * ( gh + gk ) - 1.0D+00 end if else if ( ah - ak < 0.0D+00 ) then b = 2.0D+00 * gk else b = 2.0D+00 * gh end if end if b = max ( b, 0.0D+00 ) b = min ( b, 1.0D+00 ) bivnor = b return end if sqr = sqrt ( rr ) if ( idig == 15 ) then con = twopi * 1.0D-15 / 2.0D+00 else con = twopi / 2.0D+00 do i = 1, idig con = con / 10.0D+00 end do end if ! ! (0,0) ! if ( ah == 0.0D+00 .and. ak == 0.0D+00 ) then b = 0.25D+00 + asin ( r ) / twopi b = max ( b, 0.0D+00 ) b = min ( b, 1.0D+00 ) bivnor = b return end if ! ! (0,nonzero) ! if ( ah == 0.0D+00 .and. ak /= 0.0D+00 ) then b = gk wh = -ak wk = ( ah / ak - r ) / sqr gw = 2.0D+00 * gk is = 1 ! ! (nonzero,0) ! else if ( ah /= 0.0D+00 .and. ak == 0.0D+00 ) then b = gh wh = -ah wk = ( ak / ah - r ) / sqr gw = 2.0D+00 * gh is = -1 ! ! (nonzero,nonzero) ! else if ( ah /= 0.0 .and. ak /= 0.0 ) then b = gh + gk if ( ah * ak < 0.0D+00 ) then b = b - 0.5D+00 end if wh = - ah wk = ( ak / ah - r ) / sqr gw = 2.0D+00 * gh is = -1 end if do sgn = -1.0D+00 t = 0.0D+00 if ( wk /= 0.0D+00 ) then if ( abs ( wk ) == 1.0D+00 ) then t = wk * gw * ( 1.0D+00 - gw ) / 2.0D+00 b = b + sgn * t else if ( 1.0D+00 < abs ( wk ) ) then sgn = -sgn wh = wh * wk g2 = gauss ( wh ) wk = 1.0D+00 / wk if ( wk < 0.0D+00 ) then b = b + 0.5D+00 end if b = b - ( gw + g2 ) / 2.0D+00 + gw * g2 end if h2 = wh * wh a2 = wk * wk h4 = h2 / 2.0D+00 ex = exp ( - h4 ) w2 = h4 * ex ap = 1.0D+00 s2 = ap - ex sp = ap s1 = 0.0D+00 sn = s1 conex = abs ( con / wk ) do cn = ap * s2 / ( sn + sp ) s1 = s1 + cn if ( abs ( cn ) <= conex ) then exit end if sn = sp sp = sp + 1.0D+00 s2 = s2 - w2 w2 = w2 * h4 / sp ap = - ap * a2 end do t = ( atan ( wk ) - wk * s1 ) / twopi b = b + sgn * t end if end if if ( 0 <= is ) then exit end if if ( ak == 0.0D+00 ) then exit end if wh = -ak wk = ( ah / ak - r ) / sqr gw = 2.0D+00 * gk is = 1 end do b = max ( b, 0.0D+00 ) b = min ( b, 1.0D+00 ) bivnor = b return end subroutine bivariate_normal_cdf_values ( n_data, x, y, r, fxy ) !*****************************************************************************80 ! !! BIVARIATE_NORMAL_CDF_VALUES returns some values of the bivariate normal CDF. ! ! Discussion: ! ! FXY is the probability that two variables A and B, which are ! related by a bivariate normal distribution with correlation R, ! respectively satisfy A <= X and B <= Y. ! ! Mathematica can evaluate the bivariate normal CDF via the commands: ! ! <