dbetai.f

Go to the documentation of this file.

*DECK DBETAI
      DOUBLE PRECISION FUNCTION dbetai (X, PIN, QIN)
C***BEGIN PROLOGUE  DBETAI
C***PURPOSE  Calculate the incomplete Beta function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7F
C***TYPE      DOUBLE PRECISION (BETAI-S, DBETAI-D)
C***KEYWORDS  FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C   DBETAI calculates the DOUBLE PRECISION 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 DOUBLE PRECISION.
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  D1MACH, DLBETA, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  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   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
C***END PROLOGUE  DBETAI
      DOUBLE PRECISION x, pin, qin, alneps, alnsml, c, eps, finsum, p,
     1  ps, q, sml, term, xb, xi, y, d1mach, dlbeta, p1
      LOGICAL first
      SAVE eps, alneps, sml, alnsml, first
      DATA first /.true./
C***FIRST EXECUTABLE STATEMENT  DBETAI
      IF (first) THEN
         eps = d1mach(3)
         alneps = log(eps)
         sml = d1mach(1)
         alnsml = log(sml)
      ENDIF
      first = .false.
C
      IF (x .LT. 0.d0 .OR. x .GT. 1.d0) CALL xermsg('SLATEC', 'DBETAI',
     +   'X IS NOT IN THE RANGE (0,1)', 1, 2)
      IF (pin .LE. 0.d0 .OR. qin .LE. 0.d0) CALL xermsg('SLATEC',
     +   'DBETAI', 'P AND/OR Q IS LE ZERO', 2, 2)
C
      y = x
      p = pin
      q = qin
      IF (q.LE.p .AND. x.LT.0.8d0) go to 20
      IF (x.LT.0.2d0) go to 20
      y = 1.0d0 - y
      p = qin
      q = pin
C
 20   IF ((p+q)*y/(p+1.d0).LT.eps) go to 80
C
C EVALUATE THE INFINITE SUM FIRST.  TERM WILL EQUAL
C Y**P/BETA(PS,P) * (1.-PS)-SUB-I * Y**I / FAC(I) .
C
      ps = q - aint(q)
      IF (ps.EQ.0.d0) ps = 1.0d0
      xb = p*log(y) - dlbeta(ps,p) - log(p)
      dbetai = 0.0d0
      IF (xb.LT.alnsml) go to 40
C
      dbetai = exp(xb)
      term = dbetai*p
      IF (ps.EQ.1.0d0) go to 40
      n = max(alneps/log(y), 4.0d0)
      DO 30 i=1,n
        xi = i
        term = term * (xi-ps)*y/xi
        dbetai = dbetai + term/(p+xi)
 30   CONTINUE
C
C NOW EVALUATE THE FINITE SUM, MAYBE.
C
 40   IF (q.LE.1.0d0) go to 70
C
      xb = p*log(y) + q*log(1.0d0-y) - dlbeta(p,q) - log(q)
      ib = max(xb/alnsml, 0.0d0)
      term = exp(xb - ib*alnsml)
      c = 1.0d0/(1.d0-y)
      p1 = q*c/(p+q-1.d0)
C
      finsum = 0.0d0
      n = q
      IF (q.EQ.dble(n)) n = n - 1
      DO 50 i=1,n
        IF (p1.LE.1.0d0 .AND. term/eps.LE.finsum) go to 60
        xi = i
        term = (q-xi+1.0d0)*c*term/(p+q-xi)
C
        IF (term.GT.1.0d0) ib = ib - 1
        IF (term.GT.1.0d0) term = term*sml
C
        IF (ib.EQ.0) finsum = finsum + term
 50   CONTINUE
C
 60   dbetai = dbetai + finsum
 70   IF (y.NE.x .OR. p.NE.pin) dbetai = 1.0d0 - dbetai
      dbetai = max(min(dbetai, 1.0d0), 0.0d0)
      RETURN
C
 80   dbetai = 0.0d0
      xb = p*log(max(y,sml)) - log(p) - dlbeta(p,q)
      IF (xb.GT.alnsml .AND. y.NE.0.0d0) dbetai = exp(xb)
      IF (y.NE.x .OR. p.NE.pin) dbetai = 1.0d0 - dbetai
C
      RETURN
      END