GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
dbetai.f
Go to the documentation of this file.
1 
2 *DECK DBETAI
3  DOUBLE PRECISION FUNCTION dbetai (X, PIN, QIN)
4 C***BEGIN PROLOGUE DBETAI
5 C***PURPOSE Calculate the incomplete Beta function.
6 C***LIBRARY SLATEC (FNLIB)
7 C***CATEGORY C7F
8 C***TYPE DOUBLE PRECISION (BETAI-S, DBETAI-D)
9 C***KEYWORDS FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS
10 C***AUTHOR Fullerton, W., (LANL)
11 C***DESCRIPTION
12 C
13 C DBETAI calculates the DOUBLE PRECISION incomplete beta function.
14 C
15 C The incomplete beta function ratio is the probability that a
16 C random variable from a beta distribution having parameters PIN and
17 C QIN will be less than or equal to X.
18 C
19 C -- Input Arguments -- All arguments are DOUBLE PRECISION.
20 C X upper limit of integration. X must be in (0,1) inclusive.
21 C PIN first beta distribution parameter. PIN must be .GT. 0.0.
22 C QIN second beta distribution parameter. QIN must be .GT. 0.0.
23 C
24 C***REFERENCES Nancy E. Bosten and E. L. Battiste, Remark on Algorithm
25 C 179, Communications of the ACM 17, 3 (March 1974),
26 C pp. 156.
27 C***ROUTINES CALLED D1MACH, DLBETA, XERMSG
28 C***REVISION HISTORY (YYMMDD)
29 C 770701 DATE WRITTEN
30 C 890531 Changed all specific intrinsics to generic. (WRB)
31 C 890911 Removed unnecessary intrinsics. (WRB)
32 C 890911 REVISION DATE from Version 3.2
33 C 891214 Prologue converted to Version 4.0 format. (BAB)
34 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
35 C 920528 DESCRIPTION and REFERENCES sections revised. (WRB)
36 C***END PROLOGUE DBETAI
37  DOUBLE PRECISION x, pin, qin, alneps, alnsml, c, eps, finsum, p,
38  1 ps, q, sml, term, xb, xi, y, d1mach, dlbeta, p1
39  LOGICAL first
40  SAVE eps, alneps, sml, alnsml, first
41  DATA first /.true./
42 C***FIRST EXECUTABLE STATEMENT DBETAI
43  IF (first) THEN
44  eps = d1mach(3)
45  alneps = log(eps)
46  sml = d1mach(1)
47  alnsml = log(sml)
48  ENDIF
49  first = .false.
50 C
51  IF (x .LT. 0.d0 .OR. x .GT. 1.d0) CALL xermsg('SLATEC', 'DBETAI',
52  + 'X IS NOT IN THE RANGE (0,1)', 1, 2)
53  IF (pin .LE. 0.d0 .OR. qin .LE. 0.d0) CALL xermsg('SLATEC',
54  + 'DBETAI', 'P AND/OR Q IS LE ZERO', 2, 2)
55 C
56  y = x
57  p = pin
58  q = qin
59  IF (q.LE.p .AND. x.LT.0.8d0) go to 20
60  IF (x.LT.0.2d0) go to 20
61  y = 1.0d0 - y
62  p = qin
63  q = pin
64 C
65  20 IF ((p+q)*y/(p+1.d0).LT.eps) go to 80
66 C
67 C EVALUATE THE INFINITE SUM FIRST. TERM WILL EQUAL
68 C Y**P/BETA(PS,P) * (1.-PS)-SUB-I * Y**I / FAC(I) .
69 C
70  ps = q - aint(q)
71  IF (ps.EQ.0.d0) ps = 1.0d0
72  xb = p*log(y) - dlbeta(ps,p) - log(p)
73  dbetai = 0.0d0
74  IF (xb.LT.alnsml) go to 40
75 C
76  dbetai = exp(xb)
77  term = dbetai*p
78  IF (ps.EQ.1.0d0) go to 40
79  n = max(alneps/log(y), 4.0d0)
80  DO 30 i=1,n
81  xi = i
82  term = term * (xi-ps)*y/xi
83  dbetai = dbetai + term/(p+xi)
84  30 CONTINUE
85 C
86 C NOW EVALUATE THE FINITE SUM, MAYBE.
87 C
88  40 IF (q.LE.1.0d0) go to 70
89 C
90  xb = p*log(y) + q*log(1.0d0-y) - dlbeta(p,q) - log(q)
91  ib = max(xb/alnsml, 0.0d0)
92  term = exp(xb - ib*alnsml)
93  c = 1.0d0/(1.d0-y)
94  p1 = q*c/(p+q-1.d0)
95 C
96  finsum = 0.0d0
97  n = q
98  IF (q.EQ.dble(n)) n = n - 1
99  DO 50 i=1,n
100  IF (p1.LE.1.0d0 .AND. term/eps.LE.finsum) go to 60
101  xi = i
102  term = (q-xi+1.0d0)*c*term/(p+q-xi)
103 C
104  IF (term.GT.1.0d0) ib = ib - 1
105  IF (term.GT.1.0d0) term = term*sml
106 C
107  IF (ib.EQ.0) finsum = finsum + term
108  50 CONTINUE
109 C
110  60 dbetai = dbetai + finsum
111  70 IF (y.NE.x .OR. p.NE.pin) dbetai = 1.0d0 - dbetai
112  dbetai = max(min(dbetai, 1.0d0), 0.0d0)
113  RETURN
114 C
115  80 dbetai = 0.0d0
116  xb = p*log(max(y,sml)) - log(p) - dlbeta(p,q)
117  IF (xb.GT.alnsml .AND. y.NE.0.0d0) dbetai = exp(xb)
118  IF (y.NE.x .OR. p.NE.pin) dbetai = 1.0d0 - dbetai
119 C
120  RETURN
121  END