The incomplete beta integral
Citation for published version (APA):
Newby, M. J. (1991). The incomplete beta integral. (TH Eindhoven. THE/BDK/ORS, Vakgroep ORS : rapporten; Vol. 9106). Technische Universiteit Eindhoven.
Document status and date: Published: 01/01/1991
Document Version:
Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)
Please check the document version of this publication:
• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.
• The final author version and the galley proof are versions of the publication after peer review.
• The final published version features the final layout of the paper including the volume, issue and page numbers.
Link to publication
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal.
If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:
www.tue.nl/taverne
Take down policy
If you believe that this document breaches copyright please contact us at: openaccess@tue.nl
providing details and we will investigate your claim.
.RW 16
~DK
9106
Technische Universiteit Eindhoven Faculteit Technische Bedrijfskunde
THE INCOMPLETE BETA INTEGRAL
Dr. M.J. Newby
THE INCOMPLETE BETA INTEGRAL
by Martin Newbyt
Faculty of Industrial Engineering and Management, Eindhoven University of Technology,
and
Frits Philips Institute for Quality Management, Eindhoven
Keywords: beta integral; beta distribution.
LANGUAGE Fortran 77
DESCRIPTION AND PURPOSE
This note presents a revised version of the algorithm for the incomplete beta integral x
I 1
J
p-l q-lx(p,q) ::::::; B(p,q) oU (1 - u) du O:s;x:s;l, p>O, pO.
NUMERICAL METHOD
The algorithm is computationally equivalent to AS63 (Majumder and Bhattacharjee, 1973) but is clearer and the code relates directly to the mathematical expression of the algorithm. The algorithm is based on results given in Abramowitz and Stegun (26.5.4 and 26.5.5, 1972) which are essentially the same as the reduction formula of Soper (1921). Set r ::::::; x/(l-x),
then I x(p,q) (26.5,5) and I x( p+n, q-n) x p +n (I-x) q-n {
~
}= -;-(
p-+-n-;"") Bn+-( p-+'-n..!-, q---n-:-) 1+
4-t
dj / 1 =1 (26.5.4)t Address for correspondence: Faculty of Industrial Engineering and Management, Eindhoven University of Technology, Den Dolech 2, PO Box 513, 5600 MB Eindhoven, The Netherlands
where Moreover d p+q 1 = p+n=t=I x p+n ( l_x)q-n (p+n)B(p+n,q-n)
Thus the integral can be written
p+q
----, d· p+n+~ 1
t-The length of the first series is determined by taking n
=
[q + (p + q)(l - x)], and ifn = 0 the first summation in (1) is empty. The choice of tail also follows Majumder and Bhattacharjee. The algorithm is expressed as a function subprogram which returns the value through the function name and an error flag as one of the arguments. The arguments required are the values of x, p, and qj a value for the logarithm of the beta function B(p,q) is also required. The logarithm of the beta function can be passed to the subprogram or calculated within the subprogram as necessary. The subprogram used for the logarithm of the gamma function is based on an asymptotic expansion in Abramowitz and Stegun (6.1.41, 1972) and a reduction formula.
The code is further divided for ease of understanding into a function BET AO which
calculates the sum (1) which is called by an interface routine BET AIN that handles error
checking, the simple special cases
and
and the choice of the optimal tail (Majumder and Bhattacharjee, 1973). BET AO performs no
error checking but could be called directly by another program unit. There appears to be no benefit from considering the special case in which q is an integer when the first sum
tenninates after q-l terms and the second series is not needed; this case is implicitly handled by the tests in BET AO.
STRUCTURE
REAL FUNCTION BETAIN(X,P,Q,B,IFAIL)
Formal parameters
x
Real input: the value of the upper limit xP Real input: the parameter p of the beta integral, O<P
Q
B
IFAIL
Real input: the parameter q of the beta integral, O<Q
Real input/output: the logarithm of the beta function. If B>O it is used as the value of B(p,q). If B'5.0 the value of B(p,q) is evaluated within the subprogram and is available through B for later use Integer output: error flag, IF AIL=O indicates success
IF AIL=! P'5.0 or Q'5.0
IF AIL=2 X<O or X>!
REAL FUNCTION BETAO(X,P,Q,B,TOL)
Formal parameters
X Real input: the value of the upper limit x
P Real input: the parameter p of the beta integral, O<P
Q Real input: the parameter q of the beta integral, O<Q
B Real input: the logarithm of the beta function.
TOL Real input: controls the accuracy, set in BET AIN to 1.0E-1O
REAL FUNCTION ALGAMA(X,IFAULT)
Formal Parameters
X Real input: argument of function
IFAULT Integer output: error flag, IFAULT=O, success
IFAULT=!, error
ACCURACY
The results obtained agree with those given by algorithm AS63, but the algorithm is clearer and protected against most problems of overflow or underflow. Moreover, through the use of symbolic constants established in a PARAMETER statement and generic functions only the declarations have to be changed to change from single to double precision
REFERENCES
Abramowitz, M and Stegun I (1972) Handbook of Mathematical Functions, Dover, New York. Cran, G. W., Martin, K. J. and Thomas, G. E. (1977) Remark AS RI9 and algorithm AS 109,
Applied Statistics, 26, 111-114.
Majumder, K. L. and Bhattacharjee, G. P. (1973) Algorithm AS 63. The Incomplete Beta Integral, Applied Statistics, 22, 409-411.
Soper, H. E. (1921) The Numerical Evaluation of the Incomplete Beta-Function. In Tracts for Computers No.7. Cambridge: Cambridge University Press.
REAL FUNCfION BETAIN(X,P,Q,B,IFAIL} IMPLICIT REAL (A-H,O-Z)
C********************************************************************* C
C BETA DISTRIBUTION FUNCfION C
C********************************************************************* INTEGER NOUGHT,IFAIL,IFAULT,UNITY
REAL W,X,P,Q,B,TOL,ONE,ZERO,
1 PPLUSQ
PARAMETER (ONE == 1.0EO,
1 ZERO == O.OEO, 2 NOUGHT == 0, 3 UNITY = 1, 4 TOL = 1.0E-1O) IFAIL == NOUGHT C********************************************************************* C
C CHECK THE PARAMETERS - X IN [0,1] AND P &; Q >0
C
C********************************************************************* IF (P .LE. ZERO .OR. Q .LE. ZERO) THEN
IF AIL == UNITY RETURN ENDIF BETAIN= X W=X PPLUSQ = P
+
QIF ( W .LT. ZERO .OR. W .GT. ONE ) THEN IFAIL == 2
GO TO 1000
C********************************************************************* C
C THE SIMPLE SPECIAL CASES WHEN ONE OF P OR Q IS 1, C OR X IS 0 OR 1
C
C********************************************************************* ELSEIF ( W .EQ. ZERO .OR. W .EQ. ONE) THEN
GO TO 1000
ELSEIF ( Q .EQ. ONE ) THEN BETAIN= W**P
ELSEIF ( P .EQ. ONE) THEN BETAIN= ONE-(ONE-W)**Q
C********************************************************************* C
C NOW USE THE EXPANSIONS FROM ABRAMOWITZ & STEGUN C CHOOSE THE BEST TAIL
C
C********************************************************************* ELSEIF ( P .GT. W*PPLUSQ ) THEN
IF ( B .LE. ZERO ) THEN
B = ALGAMA(P ,IF AULT)+ALGAMA( Q,IF AULT)-ALGAMA(PPLUSQ,IFAULT) ENDIF
BETAIN= BETAO(W,P,Q,B,TOL) ELSE
IF ( B .LE. ZERO ) THEN
B = ALGAMA(P,IFAULT)+ALGAMA(Q,IFAULT)-ALGAMA(PPLUSQ,IFAULT) ENDIF BETAIN= ONE-BETAO(ONE-W,Q,P,B,TOL) ENDIF 1000 CONTINUE RETURN END
REAL FUNCfION BETAO(X,P,Q,B,TOL)
C*********************************************************************
C
C FUNCfION BETAO(X,P,Q,B,TOL) - BETA DISTRIBUTION FUNCfION C USES THE SERIES EXPANSIONS FROM ABRAMOWITZ & STEGUN C RESULTS 26.5.4 AND 26.5.5
C
C********************************************************************* IMPLICIT REAL (A-H,O-Z)
INTEGER NOUGHT,NS,I
REAL P,Q,B,TOL,ONE,ZERO,U,CU,X,CHECK,
1 SO,TO,CA,CB,R,PPLUSQ,AO
PARAMETER (ONE
=
1.0EO,1 ZERO
=
O.OEO, 2 NOUGHT = 0) U=X CU = ONE - U C******************~;*************************************************** CC USE LOGARITHMS TO REDUCE CHANCE OF OVERFLOW C C********************************************************************** AO = P*LOG(U) + (Q-ONE)*LOG(CU)-LOG(P)-B SO = ONE TO
=
ONE PPLUSQ = P + Q NS = INT( Q + CU*PPLUSQ ) C********************************************************************** CC EXPANSION 26.5.5 FROM ABRAMOWITZ AND STEGUN C
C EVEN IF NS=O AT LEAST ONE REDUCfION GETS CORRECfLY DONE C BECAUSE THE LOOP DOES NOT EXECUTE IN STRICf FORTRAN 77 C C********************************************************************** CA
=
Q-ONE CB=
P+ONE R=
U/CU DO 10,I=1,NS TO=
TO*R*CA/CB SO=
SO + TO CA=
CA - ONE CB=
CB + ONE 10 CONTINUE C********************************************************************** CC CORRECfION TERM FOR EXPANSION 26.5.4
C
C**********************************************************************
TO
=
TO*U*CA/CB CA = PPLUSQ CB = CB+ONE CHECK = ABS(TO) C********************************************************************** CC EXPANSION 26.5.4 EXPRESSED AS A DO UNTIL LOOP C
C********************************************************************** 20 CONTINUE
IF ( (CHECK .LE. TOL) .AND. ( CHECK .LE. TOL*ABS(SO)) ) THEN GO TO 30 ENDIF SO
