betai.f
2.6 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
function betai (x, pin, qin)
c august 1980 version. 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
external albeta, r1mach
data eps, alneps, sml, alnsml / 4*0.0 /
c
if (eps.ne.0.) go to 10
eps = r1mach(3)
alneps = alog(eps)
sml = r1mach(1)
alnsml = alog(sml)
c
10 if (x.lt.0. .or. x.gt.1.0) call seteru (
1 35hbetai x is not in the range (0,1), 35, 1, 2)
if (pin.le.0. .or. qin.le.0.) call seteru (
1 29hbetai 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.8) go to 20
if (x.lt.0.2) go to 20
y = 1.0 - y
p = qin
q = pin
c
20 if ((p+q)*y/(p+1.).lt.eps) go to 80
c
c evaluate the infinite sum first.
c term will equal y**p/beta(ps,p) * (1.-ps)i * y**i / fac(i)
c
ps = q - aint(q)
if (ps.eq.0.) ps = 1.0
xb = p*alog(y) - albeta(ps, p) - alog(p)
betai = 0.0
if (xb.lt.alnsml) go to 40
c
betai = exp (xb)
term = betai*p
if (ps.eq.1.0) go to 40
c
n = amax1 (alneps/alog(y), 4.0)
do 30 i=1,n
term = term*(float(i)-ps)*y/float(i)
betai = betai + term/(p+float(i))
30 continue
c
c now evaluate the finite sum, maybe.
c
40 if (q.le.1.0) go to 70
c
xb = p*alog(y) + q*alog(1.0-y) - albeta(p,q) - alog(q)
ib = amax1 (xb/alnsml, 0.0)
term = exp (xb - float(ib)*alnsml)
c = 1.0/(1.0-y)
p1 = q*c/(p+q-1.)
c
finsum = 0.0
n = q
if (q.eq.float(n)) n = n - 1
do 50 i=1,n
if (p1.le.1.0 .and. term/eps.le.finsum) go to 60
term = (q-float(i-1))*c*term/(p+q-float(i))
c
if (term.gt.1.0) ib = ib - 1
if (term.gt.1.0) term = term*sml
c
if (ib.eq.0) finsum = finsum + term
50 continue
c
60 betai = betai + finsum
70 if (y.ne.x .or. p.ne.pin) betai = 1.0 - betai
betai = amax1 (amin1 (betai, 1.0), 0.0)
return
c
80 betai = 0.0
xb = p*alog(amax1(y,sml)) - alog(p) - albeta(p,q)
if (xb.gt.alnsml .and. y.ne.0.) betai = exp (xb)
if (y.ne.x .or. p.ne.pin) betai = 1.0 - betai
return
c
end