betai.f

Go to the documentation of this file.

*DECK BETAI
      REAL FUNCTION betai (X, PIN, QIN)
C***BEGIN PROLOGUE  BETAI
C***PURPOSE  Calculate the incomplete Beta function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7F
C***TYPE      SINGLE PRECISION (BETAI-S, DBETAI-D)
C***KEYWORDS  FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C   BETAI calculates the REAL incomplete beta function.
C
C   The incomplete beta function ratio is the probability that a
C   random variable from a beta distribution having parameters PIN and
C   QIN will be less than or equal to X.
C
C     -- Input Arguments -- All arguments are REAL.
C   X      upper limit of integration.  X must be in (0,1) inclusive.
C   PIN    first beta distribution parameter.  PIN must be .GT. 0.0.
C   QIN    second beta distribution parameter.  QIN must be .GT. 0.0.
C
C***REFERENCES  Nancy E. Bosten and E. L. Battiste, Remark on Algorithm
C                 179, Communications of the ACM 17, 3 (March 1974),
C                 pp. 156.
C***ROUTINES CALLED  ALBETA, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
C***END PROLOGUE  BETAI
      LOGICAL FIRST
      SAVE eps, alneps, sml, alnsml, first
      DATA first /.true./
C***FIRST EXECUTABLE STATEMENT  BETAI
      IF (first) THEN
         eps = r1mach(3)
         alneps = log(eps)
         sml = r1mach(1)
         alnsml = log(sml)
      ENDIF
      first = .false.
C
      IF (x .LT. 0. .OR. x .GT. 1.0) CALL xermsg('SLATEC', 'BETAI',
     +   'X IS NOT IN THE RANGE (0,1)', 1, 2)
      IF (pin .LE. 0. .OR. qin .LE. 0.) CALL xermsg('SLATEC', 'BETAI',
     +   'P AND/OR Q IS LE ZERO', 2, 2)
C
      y = x
      p = pin
      q = qin
      IF (q.LE.p .AND. x.LT.0.8) go to 20
      IF (x.LT.0.2) go to 20
      y = 1.0 - y
      p = qin
      q = pin
C
 20   IF ((p+q)*y/(p+1.).LT.eps) go to 80
C
C EVALUATE THE INFINITE SUM FIRST.
C TERM WILL EQUAL Y**P/BETA(PS,P) * (1.-PS)I * Y**I / FAC(I)
C
      ps = q - aint(q)
      IF (ps.EQ.0.) ps = 1.0
      xb = p*log(y) -  albeta(ps, p) - log(p)
      betai = 0.0
      IF (xb.LT.alnsml) go to 40
C
      betai = exp(xb)
      term = betai*p
      IF (ps.EQ.1.0) go to 40
C
      n = max(alneps/log(y), 4.0e0)
      DO 30 i=1,n
        term = term*(i-ps)*y/i
        betai = betai + term/(p+i)
 30   CONTINUE
C
C NOW EVALUATE THE FINITE SUM, MAYBE.
C
 40   IF (q.LE.1.0) go to 70
C
      xb = p*log(y) + q*log(1.0-y) - albeta(p,q) - log(q)
      ib = max(xb/alnsml, 0.0e0)
      term = exp(xb - ib*alnsml)
      c = 1.0/(1.0-y)
      p1 = q*c/(p+q-1.)
C
      finsum = 0.0
      n = q
      IF (q.EQ.REAL(n)) n = n - 1
      DO 50 i=1,n
        IF (p1.LE.1.0 .AND. term/eps.LE.finsum) go to 60
        term = (q-i+1)*c*term/(p+q-i)
C
        IF (term.GT.1.0) ib = ib - 1
        IF (term.GT.1.0) term = term*sml
C
        IF (ib.EQ.0) finsum = finsum + term
 50   CONTINUE
C
 60   betai = betai + finsum
 70   IF (y.NE.x .OR. p.NE.pin) betai = 1.0 - betai
      betai = max(min(betai, 1.0), 0.0)
      RETURN
C
 80   betai = 0.0
      xb = p*log(max(y,sml)) - log(p) - albeta(p,q)
      IF (xb.GT.alnsml .AND. y.NE.0.) betai = exp(xb)
      IF (y.NE.x .OR. p.NE.pin) betai = 1.0 - betai
      RETURN
C
      END