2.

3.
4.

5.

References

Abramowitz, Milton and Stegun,~ Irene A., Handbook of Mathematical

Functions, NBS Applied Mathematics Series 55, 1970, p.

Bosten, Nancy E. and Battiste, E.L., “Remark on Algorithm 179",
Communications of the ACM, Vol. 17, No. 3, 1974, pp. 156-7.

IMSL, Inc., Houston, TX. [MDBETA]

Koo, Joo 0., “Algorithm to Evaluate Incomplete Beta Function Ratio",
1984 Proceedings of the Statistical Computing Section, American

Statistical Association, p. 195-7.

Majumber, K.L. and Bhattacharjee, G.P., "The Incomplete Beta Integral",
Algorithm AS63, Applied Statistics, Vol. 22, No. 3, 1973, pp. 409-11.

NBS Core Math Library (CMLIB). [BETAI]
Numerical Algorithms Group (NAG), Downers Grove, IL. [G01BDE]

Reeve, Charles P., "Accurate Computation of the Log of the Gamma
Function", SED Note 86-1, October 1986.


Cc

GO TO 10
20 Vv = YE*U
YXEPS = YX*EPS

C--- INCREMENT PQ

c

DO 30 J = 0, JMAX
IF (V.LE.YXEPS) GO TO 40
CDFX = CDFE+V
PQ = PQ+1.0
V = (PQ+QP-1.0)*XY¥*V/PQ
30 CONTINUE
IFLAG = 2
ENDIF
40 IF (.NOT.LL) CDFX = 1.0-CDFX
ENDIF
RETURN
END


Table 1
Comparison of the number of terms (n+m) required for
convergence of the beta c.d.f. using two different
cutoff criteria (w). The integers n and m are as in
(4) and (5).

cutoff w=0 w=20
value x cod. f. (x n m  ontm n m  ont+m
— 0.10 . T 0 0 0 0 0
0.20 0.0000 0 0 0 0 0 0
0.30 0.0000 0 0 0 0 0 0
0.40 0.0001 tt} 1 1 0 1 1
0.50 0.0009 0 4 4 0 4 4
0.5882 eet
(w=20) 0.60 0.0060 0 9 9 9 1 10
0.70 0.0282 0 16 «#16 9 0 9
0.80 0.1073 0 32 32 7 0 7
0.90 0.3486 0 78 78 6 0 6
0.9091 +e en
(w=0) 0.95 0.5988 4 0 4 4 0 4
p=10 q=l_ e=10"°
cutoff w=0 w=20
value x c.d.f. (x) n mont n m  ontm
0.10 0..0000000 0 0 0" 0 0 0
0.20 0.00000010 0 2 2 0 2 2
0.30 0.00000590 0 6 6 0 6 6
0.40 0.00010485 0 ll 11 0 11 i121
0.50 0.00097656 0 #17 #17 Qo 17 #17
0.5882 ee a a a
(w=20) 0.60 0.00604662 Q 27 27 9 11 20
0.70 0.02824753 Q 42 42 9 6 15
0.80 0.10737418 0 73 73 9 2 1
0.90 0. 34867844 0 165 165 9 0 9
0.9091 eee ee
(w=0) 0.95 0.59873694 7 0 7 7 0 7


aaanaaaaagaagaaAgaAgaAaAAaAaAAAAaaNngAaaaaAnanagaAa

eanagaaaaaaaaaaaaaaaaa

SUBROUTINE CDFBET (X,P,Q,EPS,IFLAG,CDFX)

CDFBET WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING
DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG,
MARYLAND 20899

FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE BETA
DISTRIBUTION (ALSO KNOWN AS THE INCOMPLETE BETA RATIO) TO A
SPECIFIED ACCURACY (TRUNCATION ERROR IN THE INFINITE SERIES).
THE ALGORITHM, DESCRIBED IN REFERENCE 2, IS A MODIFICATION OF
THE ALGORITHM OF REFERENCE 1. THREE FEATURES HAVE BEEN ADDED:

1) A PRECISE METHOD OF MEETING THE TRUNCATION ACCURACY,

2) A CONSTANT W USED IN DETERMINING FOR WHICH XK VALUES THE
RELATION I(X,P,Q) = 1 - I(1-X,Q,P) IS TO BE USED, AND

3) A CONSTANT UFLO >= THE UNDERFLOW LIMIT ON THE COMPUTER.

SUBPROGRAMS CALLED: GAMLOG (LOG OF GAMMA FUNCTION)

CURRENT VERSION COMPLETED OCTOBER 24, 1986

REFERENCES :

1) MAJUMDER, K.L. AND BHATTACHARJEE, G.P., ‘THE INCOMPLETE BETA
INTEGRAL’, ALGORITHM AS 63, APPLIED STATISTICS, VOL. 22, NO. 3,
1973, PP. 409-411.

2) REEVE, CHARLES P., ‘AN ALGORITHM FOR COMPUTING THE BETA C.D.F.

TO A SPECIFIED ACCURACY’, STATISTICAL ENGINEERING DIVISION
NOTE 86-3, OCTOBER 1986.

DEFINITION OF PASSED PARAMETERS: ss—“<ts—‘s*~“‘“‘éTMSTMS*S*TMS
* X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL)
* P = FIRST PARAMETER OF THE BETA FUNCTION (>0) (REAL)
* Q = SECOND PARAMETER OF THE BETA FUNCTION (>0) (REAL)

* EPS

THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (>0) (REAL)

IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION:
Q ->» NO ERRORS DETECTED
1 ->» EITHER P OR @ OR EPS IS <= UFLO
2 -> NUMBER OF TERMS EVALUATED IN THE INFINITE SERIES
EXCEEDS JMAX

CDFX = THE C.D.F. EVALUATED AT X (REAL)

* INDICATES PARAMETERS REQUIRING INPUT VALUES


National Bureau of Standards N-86-3
Statistical Engineering Division October 1986 _
Gaithersburg, MD 20899 a

An Algorithm for Computing the Beta C.D.F. to a Specified Accuracy Vie

Charles P. Reeve

The cumulative distribution function (c.d.f.) of the beta distribution (also
called the incomplete beta ratio) is

x
I(xsp.q) =f tP7*(1-t)Tat/B(p,q)  (O<x<13p>05q>0) (1)
0

1
where B(p,q) = tP-1(1-t)% lat. [The left hand side of (1) is not standard
0

statistical notation but is used here for convenience.] This function is very

important in statistical computing because the c.d.f.! t, noncentral
FE ative binomial distributions are easily

obtainable from it. For example,

{1 + sign(x)1(x?/Lx?4v],1/2,v/2)}/2
I(v,x/Lv,xtv,],v,/2,v,/2),

ony P{t(v)<x}
F-> P{F(v, .v,)<x}

Bm~>  P{Bi(n,p)<x} {i -Pan-xsx#1) (x=0,1,2,.+. n-1

txen } . and

~B—> P{NBI(x,p)<y}

I(p,x,ytl) — (x=1,2,3,+-+5y=0,1,2,...).

Recent papers by Bosten and Battiste [2], Koo[4], and Majumber and
Bhattacharjee [5] have served as a basis. for computer algorithms currently
available at NBS [3,6,7]. Half-integer arguments are required in [4], andthe
implementation of [2] in CMLIB [6] was found to give erroneous results for
certain large values of p and g.

The algorithm described in this note is a modification of the algorithm in [5]
which, in the author's opinion, has two shortcomings:

1) The criterion for truncating the infinite series is based on the
first omitted term rather the sum of all the omitted terms, and

2) When x is near p/(p+q), for certain p and q, a very Jarge number
of_terms must be summed before achieving the required accuracy.

The modified algorithm, which has been strengthened in these two areas, is
based on the two recurrence relations



LOGICAL LL
DOUBLE PRECISION DP,DQ,GAMLOG
DATA wer. ¥, UFLO / §000,20.0,1E-100 /
CDFX = 0.0 “
Cc
C--- CHECK FOR VALIDITY OF ARGUMENTS P, Q, AND EPS
Cc
IF (P.LE.UFLO.OR.Q.LE.UFLO.OR.EPS.LE.UFLO) THEN
IFLAG = 1
RETURN
ENDIF
IFLAG = 0
Cc
C-~-- CHECK FOR SPECIAL CASES OF X
Cc
IF (X.LE.0.0) RETURN
IF (X.GE.1.0) THEN
CDFX = 1.0
ELSE
Cc
C--- SWITCH ARGUMENTS IF NECESSARY
Cc
LL = P+W.GE. (P+Q+2.0*W)*X
IF (LL) THEN

YX_= _1.0-XY

Py
OUR
ooh wa

C--- EVALUATE THE BETA P.D.F. AND CHECK FOR UNDERFLOW
Cc
DP = DBLE(PQ-1.0)*DLOG(DBLE(XY) )-GAMLOG( PQ)
na = DBLE(QP-1. 0) *DLOG( DBLE (YX) )~ —-GAMLOG(@P)

Pi
e (PDFL.LT.ALOGCUFLO) ) THEN
ELSE
U = EXP(PDFL)*XY/PQ
R = KY/YX
10 IF (QP.LE.1.0) GO TO 20
c
C--- INCREMENT PQ AND DECREMENT QP
c

IF (U.LE.EPS*(1.0-(PQ+QP)*XY/(PQ+1.0))) GO TO 40
CDFX = CDFX+U

PQ = PQ+1.0

QP = QP-1.0

U = QP*R*U/PQ

yeh : a-w* “

& Cb)

band %


1

P(y_y) 9°
I(x,p.q) = ie + I(x,ptl,q-1) — (p>03q>1) (2)

oa ”

Peqexy .
1(x,p.q) = 2SRe%LS + T(xptlsq) — (p205920) (3)

given by equations 26.5.15 and 26.5:16 of Abramowitz and Stegun[1]. In both
cases the lower tail area of one beta distribution is expressed in terms of
the lower tail area of another beta distribution. Applying (2) n times yields

n pti-1 1- q-i
I(x,p.q) = ah ai B atic qi + I(x,ptn,q-n) (4)
provided that O<n<q. Applying (3) mtimes to the rightmost term in (4),
n xPtinl yy ard m xPtmtj-2eq_ yarn

I(x,p.q) = oh (pt i-T)B(pHi-1,q-te1y + sh (penty-T)B(penej-lvqeny (5)

+ I(x,ptn+m,q-n)

where m0. For a given x<1, p, q, and lute accuracy limit e>0 there will
be non-negative integers n<q and m such that I(x,ptntm,q-n)<e (the proof is
left to the reader). The summation of ntm terms will then be an acceptable
approximation to I(x,p,q). Note that the second summation in (5) is needed
only if the first summation is carried to completion (0<q-n<1) without

achieving the desired absolute accuracy.

An upper bound for I(x,p,q) is computed by applying (3) repeatedly to obtain

P q 2 : ;
x wae + +gtl) ... (ptgti-1
1(x,p.q) = pB(p.q [1 + ah ptl)(p+2) ... (pti x]. (6)

If q<l then from (6)

P(1-x)4 ei
I(x,psq) < ae (l+]x
PB (Psa dy

(7)
< Pag!
pB(p,q)(I=x) ”
when x<l. If q>1 then
p -y)4 2 as
1(x,p.q) < a gh [1+ 2 Gab]
(8)

< xP(1-x)4
pB(p.q)[1- (ptq)x/(p+1)]



when x<(p+1)/(p+q). Combining the bounds in (7) and (8) yields the bound
for the remainder term in (5),

“ p+n+m qen -
x (1-x)
T(x,p+ném.q-n) § Tormem)B(pentm,q-n) (1-x/P) (9)

where r=min[1,(ptntm+1)/(ptqtm)] provided that x<r. When the right hand side
of (9) becomes <e as n or m increases, the series which is being summed is
truncated.

Although the above algorithm works for O<x<1, the number of terms required may
become very larae as 7 In this case a more efficient procedure js to
compute I(1-x,q,p) and then use the identity

I(x,p,q) = 1 - I(1-x,q,p)-

This is done in [5] when xp/(p+q). However, when x is slightly less than
this cutoff value, a large number of terms may be required for convergence.
The author conducted a study using cutoff values of the form (ptw)/(p+qt2w)
and found that the value w=20 did the best job of minimizing the number of
terms required for convergence for a broad range of p and q values. The
increased efficiency using this new cutoff value is illustrated in table 1 for
the case p=10 and q=l.

The above algorithm with w=20 has been implemented in the FORTRAN subroutine
CDFBET. If x<0 then 0 is returned, and if x>1 then 1 is returned. If x is in
the extreme tails of the distribution (near 0 or 1), exponential underflow
might occur when computing the first term (i=1) in (5). This is prevented by
testing the exponent to see if it is less than a pre-set underflow limit. If
it is then I(x,p,q) is set to 0 or 1 as appropriate. An external routine,
such as Reeve[8] must be provided for computing the log of the gamma function

jin_double_ precision.
A listing of CDFBET is an appendix to this note. It is invoked by

CALL CDFBET(X,P,Q,EPS,IFLAG, CDFX)

where the variable names are defined in the program documentation. The
returned value of CDFX is valid only if IFLAG=0 on return. In passing «
(variable name EPS) to CDFBET the user should realize that accuracy is limited
by the number of digits carried in a single precision variable, and that
roundoff error may affect the last one or two of these digits.