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