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