blob d10f573c2667925539ea35fd65c2b82f6f2f6b63
1 #include "xmath.h"
3 #define twopi 6.283195307179587
4 #define con (twopi / 2.0) * 10.0e-10
6 extern void normbase(double *, double *);
8 double bivnor(double ah, double ak, double r)
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)
15 Tranlated from FORTRAN to ratfor using struct; from ratfor to C by hand.
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;
21 temp = -ah;
22 normbase(&temp, &gh);
23 gh = gh / 2.0;
24 temp = -ak;
25 normbase(&temp, &gk);
26 gk = gk / 2.0;
28 b = 0;
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;
45 else if (ak==0) {
46 b = atan(r/sqr)/twopi+.25;
47 goto label50;
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;
80 b = b-(gw+g2)/2+gw*g2;
83 h2 = wh*wh;
84 a2 = wk*wk;
85 h4 = h2*.5;
86 ex = 0;
87 if (h4<150.0) {
88 ex = exp(-h4);
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;
112 if (is>=0)
113 break;
114 } while(ak!=0);
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;
124 label50:
125 if (b<0)
126 b = 0;
127 if (b>1)
128 b = 1;
130 return(b);