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 MIT 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 = rk ) AH, AK, the lower limits of integration. ! ! Input, real ( kind = rk ) R, the correlation between X and Y. ! ! Output, real ( kind = rk ) BIVNOR, the bivariate normal CDF. ! ! Local: ! ! integer IDIG, the number of significant digits ! to the right of the decimal point desired in the answer. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a2 real ( kind = rk ) ah real ( kind = rk ) ak real ( kind = rk ) ap real ( kind = rk ) b real ( kind = rk ) bivnor real ( kind = rk ) cn real ( kind = rk ) con real ( kind = rk ) conex real ( kind = rk ) ex real ( kind = rk ) g2 real ( kind = rk ) gauss real ( kind = rk ) gh real ( kind = rk ) gk real ( kind = rk ) gw real ( kind = rk ) h2 real ( kind = rk ) h4 integer i integer, parameter :: idig = 15 integer is real ( kind = rk ) r real ( kind = rk ) rr real ( kind = rk ) s1 real ( kind = rk ) s2 real ( kind = rk ) sgn real ( kind = rk ) sn real ( kind = rk ) sp real ( kind = rk ) sqr real ( kind = rk ) t real ( kind = rk ), parameter :: twopi = 6.283185307179587D+00 real ( kind = rk ) w2 real ( kind = rk ) wh real ( kind = rk ) 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: ! ! <