ON THE COMPUTATION OF THE NONCENTRAL BETA-DISTRIBUTION

Computational Statistics & Data Analysis 18 (1994) 449-455 North-Holland

449

On the computation

of the noncentral beta distribution

Cherng G. Ding

Institute of Management Science, National Chiao-Tung University, Taipei, Taiwan, Republic of China

Received September 1991 Revised April 1993

Abstract: Two recurrence formulas for computing the noncentral beta distribution function and the density are proposed. The computations can be carried to any specified accuracy using the error bounds obtained. Corresponding algorithms are provided in a step-by-step form. A close relationship between the formulas allows joint evaluation of the distribution function and the density. The noncentral beta quantile can then be efficiently computed using Newton’s method. Its algorithm is also given in a step-by-step form.

Keywords: Error bound; Gamma function; Mean-value theorem; Newton’s method; Noncentral beta distribution; Quantile; Series representation

1. Introduction

Let F(x; a, b, A) and f(x; a, b, A) denote, respectively, the noncentral beta distribution function and the density with shape parameters a, b, and noncen- tral&y parameter A, and F(x; a, b) and f(x; a, b) the central beta distribution function and the density with shape parameters a and b. It is well known that (see, e.g., Tang, 1938)

F(x; a, b, A) = 5 u,F(x; a +i, b), ₍₁₎

i=o and

f(x; a, b, A) = 2 z$f(x; a + i, b), ₍₂₎

i=O

whereO<x<l, u>O, b>O, h>O,and Ui=e -A/2(iA>‘/i!. The evaluations of F(x; a, b, A), f( x; a, b, A), and the quantile (percentage point) xP satisfying F(x,; a, b, A) =p for given p are always of interest. Lenth (1987) developed a

Correspondence to: Cherng G. Ding, Institute of Management Science, National Chiao-Tung University, 4F, 114 Chung-Hsiao W. Road, Section 1, Taipei, Taiwan, Republic of China. 0167-9473/94/$07.00 0 1994 - Elsevier Science B.V. All rights reserved

recursive algorithm for computing F(x; a, b, A), which was a modification of one given in Norton (1983). Given an accurate auxiliary routine for evaluating F(x; a, b), both of the algorithms summed up the terms in formula (1) until the corresponding error bounds were less than some desired accuracy. Singh and Relyea (1992) developed a computer program for computing F(x; a, b, h) where 2a and 2b are both integers based on the error bound given in Lenth (1987) and exact expressions for F(x; a, b). In this paper, an alternative series representation for F(x; a, b, A) in terms of the central beta densities is developed by using the method similar to that given in Ding (1992). It is easy to compute recursively for all real shape parameters without any auxiliary central beta routine. A simple recursive formula for evaluating f<x; a, b, A) is also given. The accuracy can be easily controlled by using the error bounds provided. The quantile X, is obtained by Newton’s method due to the advantage that the proposed computing formulas for F(x; a, b, A) and f(x; a, b, h) are closely related, and hence they can be evaluated jointly, rather than separately.

2. Numerical methods

Let T(a) (a > 0) denote the Gamma function. It is easy, using integration by parts, to show that (see also Abramowitz and Stegun, 1965, p. 944)

qa +

b)

qx;

a,

b) =

r(a +

l)r(b)xa(l

-x)”

+qx;

a + 1, b)

=f(x; a + 1, b)(l -X)/(U +b) +F(x, a + 1, b). ₍₃₎ The above relation implies that F(x; a + i, b) can be expressed in terms of a series given by

F(x;

a+i,

b)= f

(l-4

f( k=j

(a +

b + k)

x,u+k+l,b)= et,, i=o, l,...,

k=i

(4)

where the terms t, can be obtained recursively by r(u + b)

to = r(u + l)T(b) x”(1 -x)b,

ti=ti_,x(u+b+i-l)/(a+i), i21. ₍₅₎

Substituting (4) into (11, we have a new series representation for F(x; a, b, A), which involves the central beta densities:

F(x;a,b,A)= eui

i=O

(6)

C.G. Ding / Computation of the noncentral beta distribution 451

where the terms are evaluated by u0 = u0 = e-*1*,

Ui = Ui_I + ui, Ui = Ui_$/(2i), i2 1, ₍₇₎

tiY i 2 0, as in (5).

Considering equation (4) for i = 0, we have

F(x; a, b) = 2 t, = yt, +qx; a +n, b), ₍₈₎

k=O k=O

where F(x; a + tz, b) = Cycntk is the error of truncation after II terms. It can be shown, by the mean-value theorem, that F(x; a + IZ, b) is less than or equal to t,_, x(a+b+n-l)/[(a+n)-(a+b+n)x] if u+n>(u+b+n)x. Let

EF, be the error of truncation after II terms for the series for F(x; a, b, A) in

(6). Since ui I 1 for all i 2 0, it follows that, if a + n > (a + b + n>x,

EF,= &,+tist,,x(u+b+n-l),[(u+n)-(u+b+n)x]. ₍₉₎

i=n

Note that the error bound given above is a decreasing function of yt when a + y1> (a + b + n>x, and is a common one used to control the accuracy of evaluations of F(x; a, b, A) and F(x; a, b). The terms in (61, computed recursively through (71, can then be summed up until the error bound is not greater than a predetermined small number E. In fact, evaluation of F(x; a, b,

A) through formula (6) requires no central beta routine. Computing the error

bound is very easy because it relies only upon the factor t,_, of the last term of the truncated series Er:J ui ti.

The noncentral beta density f(x; a, b, A) expressed in (2) is an infinite weighted sum of central beta densities. Its terms can be evaluated recursively, and are related to those of (6) as follows:

f(x; a, b, A) = 2 uif(x; a + i, b) = 2 uisi, ₍₁₀₎

i=O i=O

where

r(u + b)

So = T(u)r(b)x a-‘(1 -X)b-l =ut,/x/(l -x),

si=si_,x(u+b+i-l)/(u+i-l)=ti_,(u+b+i-1)/(1-x), i.21, (11) ui and ti, i 2 0, are those in (7).

It is clear that, when a + II > (a + b + n>x, the sequence {sJ(i 2 n) is decreas- ing. Letting Ef, be the truncation error after IZ terms for the series in (lo), we have, when a + n > (a + b + 11)x,

Ef,,= &+si< &,s,=s, =x,(1-uu,_i). ₍₁₂₎

Likewise, the above error bound is a decreasing function of IZ, and is used to control the accuracy of evaluation of f(x; a, b, A). It is also easy to compute since it relies only upon factors that are computed in the course of evaluating (10).

Let G(x) = F(x; a, b, A) -p. The quantile xp is to be obtained by solving the equation G(x) = 0. Since computing formulas for F(x; a, b, A) and G’(x) =f(x; a, b, A) are closely related (the terms in (10) involve the factors ui same as those for evaluating (6) and the factors si whose computations can be based on t, in (7)), their evaluations can be combined. It is this property that makes Newton’s method for computing the root x, much more efficient because many redun- dant computations are avoided. The process is to repeat the iteration (see, e.g., Kennedy and Gentle, 1980, pp. 72-73; and Vandergraft, 1983, pp. 262-264)

xj+l

=xj-

[F(x;;

a,

b,

“)

-P]

f(xj; a, b, A) ’

j=o,

1 ‘***’

(13)

until the relative absolute difference between two successive iterations is not greater than a specified accuracy 6, i.e., 1 x~+~ -xi I I SX~+~. Since G(0) = -p,

G(1) = 1 -p, and G is strictly increasing in [0, 11, the solution xp of G(x) = 0 is unique. It is obvious that xp = 0 for p = 0 and x, = 1 for p = 1. No computation is needed for these cases. For 0 <p < 1, perform iterations with the starting value x0 = 0.5. Within each iteration, the condition that 0 <xj+r < 1 must be satisfied. If not, replace xj+r by ixj if it is nonpositive, and by +(xj + 1) if it is greater than or equal to 1. This modified rule is easy to implement, and the iterations should converge fast for most cases. Note that evaluations of F( x; a, b, A) and f(x; a, b, A) should be precise enough so that the accuracy of Newton’s solution can be warranted. Also, the number of iterations needs to be controlled.

3. Algorithms

Given an accurate routine for computing (the logarithm of) the Gamma function (see Pike and Hill, 1966; or Macleod, 1989) the formulas discussed in Section 2 lead to three simple algorithms for computing the noncentral beta distribution function, the density, and the quantile to any desired accuracy: Algorithm A (noncentrul beta distribution function). This algorithm computes the

noncentral beta distribution function F(x; a, b, A) for given abscissa x, shape parameters a, b, and noncentrality parameter A.

Al (specify the accuracy). Set EPS + E. (E denotes the desired accuracy, e.g., 10-6.)

A2 (initialize). Set n + 1, t + (T(u + b)/T(u + l)T(b))xa(l -xjb, u + e-“12, v +- u, CDF + vt.

C.G. Ding / Computation of the noncentral beta distribution 453

A4 (update the term and then accumulate). Set u 6 uh/(2n), u + u + U, t + tx(a + b + II - l)/(a + n), CDF t CDF + vt, n + n + 1, and return to step A3.

A5 (find the error bound and check for convergence). Set bound +- tx(a + b + n - l)/((a + n) - (a + b + n)x). If bound I EPS, terminate the algorithm (CDF is the answer.)

A6 (update the term and then accumulate). Set u + uh/(2n), u + u + u, t + &(a + b + n - l)/(u + n), CDF + CDF + ut, n +- n + 1, and return to step A5.

Algorithm B (noncentral beta density). This algorithm computes the noncentral beta density f(x; a, b, A) for given abscissa x, shape parameters a, b, and noncentrality parameter A.

Bl (specify the accuracy). Set EPS + E. (E denotes the desired accuracy, e.g., 10-6.)

B2 (initialize). Set II +- 1, s +- (T(u + b)/r(u)IYb))x”-‘(1 -xjb-‘, u + e-‘/=, v + u, PDF + us.

B3 (compare a + it, (a + b + n>x). If a + it > (a + b + n)x, go to step B5. B4 (update the term and then accumulate). Set u + uA/(2n), u + u + u, s + sx(u + b + n - l)/(u + II - l), PDF + PDF + us, n + n + 1, and return to

step B3.

B5 (find the error bound and check for convergence). Set bound + sx(u + b + n - l)(l - ~)/(a + II - 1). If bound I EPS, terminate the algorithm. (PDF is the answer.)

B6 (update the term and then accumulate). Set u + uA/(2n), u + u + u, s + sx(u + b + II - l)/(u + II - 1), PDF + PDF + us, II + II + 1, and return to step B5.

Algorithm C (noncentral beta quuntile). This algorithm computes the noncentral beta quantile xP such that F(x,; a, b, A) =p.

Cl (specify the accuracy and the maximum number of Newton’s iterations). Set EPS + E, DELTA + 6, ITRMAx+ N,,. (E denotes the desired accuracy (e.g., 10m6) for computing F(x; a, b, A) and f(x; a, b, A), and 6 the desired accuracy (e.g., 10-4) for computing x,. N,, is an integer (e.g., 10) used to limit the number of Newton’s iterations.)

C2 (set the constants). Set cOeff+- T(u + b)/T(u + l)T(b), u0 + e-“I’. C3 (loop on j, Newton’s iteration). First initialize x by setting x + 0.5. Then perform steps C4 through Cl0 for j = 1, 2,. . . , ITRiVZAX. (Steps C4 through Cl0 constitute one iteration.)

C4 (initialize within each iteration). Set it + 1, t + coeff*x”(l -x)~, s + at/x/(1 -x), u +- ug, v + u, CDF t vt, PDF + us.

C5 (compare a + ~1, (a + b + n)x>. If a + II > (a + b + n)x, go to step C7. C6 (update the terms and then accumulate for both of F(x; a, b, A) and f(x; a, b, A)). Set u + uA/(2n), v + v + u, s + t(u + b + II - l)/(l -x), t + tx(u +

b + n - l>/(a + n), CDF + CDF + ut, PDF + PDF + us, n + n + 1, and return

to step C5.

C7 (find the corresponding error bounds). Set bound1 + tx(a + b + n -

l>/((a + n) - (a + b + n>x), bound2 + t(a + b + n - l)(l - v)/(l -x).

C8 (check for convergence). If bound1 5 EPS and bound2 I EPS, go to step ClO.

C9 (update the terms and then accumulate for F(x; a, b, A) and/or f(x; a, b, A)). Set u 6 uh/(2n), u +u+u. If boundl<EPS, set s+sx(u+b+n- l)/(u+n-11, PDF+PDF+us, n+n+l, bound2+sx(u+b+n-Ml-

~>/(a + n - l), and return to step C8; otherwise if bound2 I EPS, set t + tx(u + b + n - l)/(u + n), CDF +CDF+ut, n+n+l, boundl+tx(u+b+n-

l)/((u + n) - (a + b + n>x), and return to step C8; otherwise set s + t(a + b +

n - l)/(l -x>, t + &(a + b + n - l>/(u + n>, CDF + CDF + vt, PDF +- PDF + us, n + n + 1, and return to step C7.

Cl0 (find new x and check for convergence of Newton’s iterations). Set

diff+ (CDF -p)/PDF. If x - diffs 0, set x + ix; otherwise if x - diff2 1, set x + $(_x + 1); otherwise set x +x - diff. If 1 diff 1 /x I DELTA, terminate the

algorithm. (x is the answer.)

Cl1 (output error message). Terminate the algorithm with the message “No convergence after iV,, iterations”.

FORTRAN codes based on the above algorithms are available upon request.

Acknowledgements

This work was partially supported by the National Science Council, Republic of China. The author is grateful to Professor Chyan Yang for helpful discussions regarding the algorithms and to a referee for helpful comments which improved the presentation of the algorithms.

References

Abramowitz, M. and I.A. Stegun, Handbook of mathematical functions (Dover, New York, 1965).

Ding, C.G., Algorithm AS 275: Computing the non-central x2 distribution function, Appl. Statist., 41 (1992) 478-482.

Kennedy, W.J. and J.E. Gentle, Statistical computing (Marcel Dekker, New York, 1980). Lenth, R.V., Algorithm AS 226: Computing noncentral beta probabilities, Appl. Statist., 36 (1987)

241-244.

Macleod, A.J., Algorithm AS 245: A robust and reliable algorithm for the logarithm of the gamma function, Appl. Statist., 38 (1989) 397-402.

Norton, V., A simple algorithm for computing the non-central F distribution, App. Statist., 32 (1983) 84-85.

Pike, M.C. and I.D. Hill, Algorithm 291: Logarithm of the gamma function, Commun. Ass. Comput. Mach., 9 (1966) 684.

C.G. Ding / Computation of the noncentral beta distribution 4.55

Singh, K.P. and G.E. Relyea, Computation of noncentral f probabilities: a computer program,

Comp. Stat. Data Anal., 13 (1992) 95-102.

Tang, P.C., The power function of the analysis of variance tests with tables and illustrations of their use, Statist. Res. Mem., 2 (1938) 126-150.

Vandergraft, J.S., Introduction to numerical computations, 2nd. ed. (Automated Sciences Group, Silver Spring, MD, 1983).