lib/bivnor.c

   1 #include "xmath.h"
   2
   3 #define twopi 6.283195307179587
   4 #define con (twopi / 2.0) * 10.0e-10
   5
   6 extern void normbase(double *, double *);
   7
   8 double bivnor(double ah, double ak, double r)
   9 {
  10   /*
  11     based on alg 4628 comm. acm oct 73
  12     gives the probability that a bivariate normal exceeds (ah,ak).
  13     gh and gk are .5 times the right tail areas of ah, ak under a n(0,1)
  14
  15     Tranlated from FORTRAN to ratfor using struct; from ratfor to C by hand.
  16   */
  17   double a2, ap, b, cn, conex, ex, g2, gh, gk, gw, h2, h4, rr, s1, s2,
  18     sgn, sn, sp, sqr, t, temp, w2, wh, wk;
  19   int is;
  20
  21   temp = -ah;
  22   normbase(&temp, &gh);
  23   gh = gh / 2.0;
  24   temp = -ak;
  25   normbase(&temp, &gk);
  26   gk = gk / 2.0;
  27
  28   b = 0;
  29
  30   if (r==0)
  31     b = 4*gh*gk;
  32   else {
  33     rr = 1-r*r;
  34     if (rr<0)
  35       return(-1.0);
  36     if (rr!=0) {
  37       sqr = sqrt(rr);
  38       if (ah!=0) {
  39         b = gh;
  40         if (ah*ak<0)
  41           b = b-.5;
  42         else if (ah*ak==0)
  43           goto label10;
  44       }
  45       else if (ak==0) {
  46         b = atan(r/sqr)/twopi+.25;
  47         goto label50;
  48       }
  49       b = b+gk;
  50       if (ah==0)
  51         goto label20;
  52     label10:
  53       wh = -ah;
  54       wk = (ak/ah-r)/sqr;
  55       gw = 2*gh;
  56       is = -1;
  57       goto label30;
  58     label20:
  59       do {
  60         wh = -ak;
  61         wk = (ah/ak-r)/sqr;
  62         gw = 2*gk;
  63         is = 1;
  64       label30:
  65         sgn = -1;
  66         t = 0;
  67         if (wk!=0) {
  68           if (fabs(wk)>=1) {
  69             if (fabs(wk)==1) {
  70               t = wk*gw*(1-gw)/2;
  71               goto label40;
  72             } else {
  73               sgn = -sgn;
  74               wh = wh*wk;
  75               normbase(&wh, &g2);
  76               wk = 1/wk;
  77               if (wk<0) {
  78                 b = b+.5;
  79               }
  80               b = b-(gw+g2)/2+gw*g2;
  81             }
  82           }
  83           h2 = wh*wh;
  84           a2 = wk*wk;
  85           h4 = h2*.5;
  86           ex = 0;
  87           if (h4<150.0) {
  88             ex = exp(-h4);
  89           }
  90           w2 = h4*ex;
  91           ap = 1;
  92           s2 = ap-ex;
  93           sp = ap;
  94           s1 = 0;
  95           sn = s1;
  96           conex = fabs(con/wk);
  97           do {
  98             cn = ap*s2/(sn+sp);
  99             s1 = s1+cn;
 100             if (fabs(cn)<=conex)
 101               break;
 102             sn = sp;
 103             sp = sp+1;
 104             s2 = s2-w2;
 105             w2 = w2*h4/sp;
 106             ap = -ap*a2;
 107           } while (1);
 108           t = (atan(wk)-wk*s1)/twopi;
 109         label40:
 110           b = b+sgn*t;
 111         }
 112         if (is>=0)
 113           break;
 114       } while(ak!=0);
 115     }
 116     else if (r>=0)
 117       if (ah>=ak)
 118         b = 2*gh;
 119       else
 120         b = 2*gk;
 121     else if (ah+ak<0)
 122       b = 2*(gh+gk)-1;
 123   }
 124  label50:
 125   if (b<0)
 126     b = 0;
 127   if (b>1)
 128     b = 1;
 129
 130   return(b);
 131 }