dbetai.f
2.87 KB
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
double precision function dbetai (x, pin, qin)
c august 1980 edition. w. fullerton, c3, los alamos scientific lab.
c based on bosten and battiste, remark on algorithm 179, comm. acm,
c v 17, p 153, (1974).
c
c input arguments --
c x upper limit of integration. x must be in (0,1) inclusive.
c p first beta distribution parameter. p must be gt 0.0.
c q second beta distribution parameter. q must be gt 0.0.
c betai the incomplete beta function ratio is the probability that a
c random variable from a beta distribution having parameters
c p and q will be less than or equal to x.
c
double precision x, pin, qin, alneps, alnsml, c, eps, finsum, p,
1 ps, p1, q, sml, term, xb, xi, y, dint, d1mach, dlbeta,
2 dexp, dlog
external d1mach, dlbeta
data eps, alneps, sml, alnsml / 4*0.0d0 /
c
if (eps.ne.0.0d0) go to 10
eps = d1mach(3)
alneps = dlog (eps)
sml = d1mach(1)
alnsml = dlog (sml)
c
10 if (x.lt.0.d0 .or. x.gt.1.d0) call seteru (
1 35hdbetai x is not in the range (0,1), 35, 1, 2)
if (pin.le.0.d0 .or. qin.le.0.d0) call seteru (
1 29hdbetai a and/or b is le zero, 29, 2, 2)
c
y = x
p = pin
q = qin
if (q.le.p .and. x.lt.0.8d0) go to 20
if (x.lt.0.2d0) go to 20
y = 1.0d0 - y
p = qin
q = pin
c
20 if ((p+q)*y/(p+1.d0).lt.eps) go to 80
c
c evaluate the infinite sum first. term will equal
c y**p/beta(ps,p) * (1.-ps)-sub-i * y**i / fac(i) .
c
ps = q - dint(q)
if (ps.eq.0.d0) ps = 1.0d0
xb = p*dlog(y) - dlbeta(ps,p) - dlog(p)
dbetai = 0.0d0
if (xb.lt.alnsml) go to 40
c
dbetai = dexp (xb)
term = dbetai*p
if (ps.eq.1.0d0) go to 40
n = max1 (sngl(alneps/dlog(y)), 4.0)
do 30 i=1,n
xi = i
term = term * (xi-ps)*y/xi
dbetai = dbetai + term/(p+xi)
30 continue
c
c now evaluate the finite sum, maybe.
c
40 if (q.le.1.0d0) go to 70
c
xb = p*dlog(y) + q*dlog(1.0d0-y) - dlbeta(p,q) - dlog(q)
ib = max1 (sngl(xb/alnsml), 0.0)
term = dexp (xb - dble(float(ib))*alnsml )
c = 1.0d0/(1.d0-y)
p1 = q*c/(p+q-1.d0)
c
finsum = 0.0d0
n = q
if (q.eq.dble(float(n))) n = n - 1
do 50 i=1,n
if (p1.le.1.0d0 .and. term/eps.le.finsum) go to 60
xi = i
term = (q-xi+1.0d0)*c*term/(p+q-xi)
c
if (term.gt.1.0d0) ib = ib - 1
if (term.gt.1.0d0) term = term*sml
c
if (ib.eq.0) finsum = finsum + term
50 continue
c
60 dbetai = dbetai + finsum
70 if (y.ne.x .or. p.ne.pin) dbetai = 1.0d0 - dbetai
dbetai = dmax1 (dmin1 (dbetai, 1.0d0), 0.0d0)
return
c
80 dbetai = 0.0d0
xb = p*dlog(dmax1(y,sml)) - dlog(p) - dlbeta(p,q)
if (xb.gt.alnsml .and. y.ne.0.0d0) dbetai = dexp(xb)
if (y.ne.x .or. p.ne.pin) dbetai = 1.0d0 - dbetai
c
return
end