(*-------------------------------------------------------------------------*)
(*               CDBeta  -- Cumulative Beta Distribution                   *)
(*-------------------------------------------------------------------------*)

FUNCTION CDBeta(     X,     Alpha, Beta: REAL;
                     Dprec, MaxIter    : INTEGER;
                 VAR Cprec             : REAL;
                 VAR Iter              : INTEGER;
                 VAR Ifault            : INTEGER  ) : REAL;

(*-------------------------------------------------------------------------*)
(*                                                                         *)
(*       Function:  CDBeta                                                 *)
(*                                                                         *)
(*       Purpose:   Evaluates CPDF of Incomplete Beta Function             *)
(*                                                                         *)
(*       Calling Sequence:                                                 *)
(*                                                                         *)
(*            P     := CDBeta(     X, Alpha, Beta: REAL;                   *)
(*                                 Dprec, Maxitr : INTEGER;                *)
(*                             VAR Cprec         : REAL;                   *)
(*                             VAR Iter          : INTEGER;                *)
(*                             VAR Ifault        : INTEGER  ) : REAL;      *)
(*                                                                         *)
(*                 X      --- Upper percentage point of PDF                *)
(*                 Alpha  --- First shape parameter                        *)
(*                 Beta   --- Second shape parameter                       *)
(*                 Dprec  --- Number of digits of precision required       *)
(*                 Maxitr --- Maximum number of iterations                 *)
(*                 Cprec  --- Actual resulting precision                   *)
(*                 Iter   --- Iterations actually used                     *)
(*                 Ifault --- error indicator                              *)
(*                            = 0:  no error                               *)
(*                            = 1:  argument error                         *)
(*                                                                         *)
(*                 P      --- Resultant probability                        *)
(*                                                                         *)
(*       Calls:                                                            *)
(*                                                                         *)
(*            ALGama                                                       *)
(*                                                                         *)
(*       Method:                                                           *)
(*                                                                         *)
(*            The continued fraction expansion as given by                 *)
(*            Abramowitz and Stegun (1964) is used.  This                  *)
(*            method works well unless the minimum of (Alpha, Beta)        *)
(*            exceeds about 70000.                                         *)
(*                                                                         *)
(*            An error in the input arguments results in a returned        *)
(*            probability of -1.                                           *)
(*                                                                         *)
(*-------------------------------------------------------------------------*)

VAR
   Epsz : REAL;
   A    : REAL;
   B    : REAL;
   C    : REAL;
   F    : REAL;
   Fx   : REAL;
   Apb  : REAL;
   Zm   : REAL;
   Alo  : REAL;
   Ahi  : REAL;
   Blo  : REAL;
   Bhi  : REAL;
   Bod  : REAL;
   Bev  : REAL;
   Zm1  : REAL;
   D1   : REAL;
   Aev  : REAL;
   Aod  : REAL;

   Ntries : INTEGER;

   Qswap  : BOOLEAN;
   Qdoit  : BOOLEAN;
   Qconv  : BOOLEAN;

LABEL   20, 9000;

BEGIN (* CdBeta *)

                                   (* Initialize *)
   IF Dprec > MaxPrec THEN
      Dprec := MaxPrec
   ELSE IF Dprec <= 0 THEN
      Dprec := 1;

   Cprec  := Dprec;

   Epsz   := PowTen( -Dprec );

   X      := X;
   A      := Alpha;
   B      := Beta;
   QSwap  := FALSE;
   CDBeta := -1.0;
   Qdoit  := TRUE;
                                   (* Check arguments *)
                                   (* Error if:       *)
                                   (*    X <= 0       *)
                                   (*    A <= 0       *)
                                   (*    B <= 0       *)

   Ifault := 1;

   IF( X <= 0.0 ) THEN GOTO 9000;

   IF( ( A <= 0.0 ) OR ( B <= 0.0 ) ) THEN GOTO 9000;

   CDBeta := 1.0;
   Ifault := 0;
                                   (* If X >= 1, return 1.0 as prob *)

   IF( X >= 1.0 ) THEN GOTO 9000;

                                   (* If X > A / ( A + B ) then swap *)
                                   (* A, B for more efficient eval.  *)

   IF( X > ( A / ( A + B ) ) ) THEN
      BEGIN
         X      := 1.0 - X;
         A      := Beta;
         B      := Alpha;
         QSwap  := TRUE;
      END;

                                   (* Check for extreme values *)

   IF( ( X = A ) OR ( X = B ) ) THEN GOTO 20;
   IF( A = ( ( B * X ) / ( 1.0 - X ) ) ) THEN GOTO 20;
   IF( ABS( A - ( X * ( A + B ) ) ) <= Epsz ) THEN GOTO 20;

         C     := ALGama( A + B ) + A * LN( X ) +
                  B * LN( 1.0 - X ) - ALGama( A ) - ALGama( B ) -
                  LN( A - X * ( A + B ) );

         IF( ( C < -36.0 ) AND QSwap ) THEN GOTO 9000;

            CDBeta := 0.0;
            IF(  C < -180.0 ) THEN GOTO 9000;

                                   (*  Set up continued fraction expansion *)
                                   (*  evaluation.                         *)

20:
      Apb    := A + B;
      Zm     := 0.0;
      Alo    := 0.0;
      Bod    := 1.0;
      Bev    := 1.0;
      Bhi    := 1.0;
      Blo    := 1.0;

      Ahi    := EXP( ALGama( Apb ) + A * LN( X ) +
                     B * LN( 1.0 - X ) - ALGama( A + 1.0 ) -
                     ALGama( B ) );

      F      := Ahi;

      Iter   := 0;

                                   (* Continued fraction loop begins here. *)
                                   (* Evaluation continues until maximum   *)
                                   (* iterations are exceeded, or          *)
                                   (* convergence achieved.                *)

      Qconv  := FALSE;

      REPEAT

         Fx     := F;

         Zm1    := Zm;
         Zm     := Zm + 1.0;
         D1     := A + Zm + Zm1;
         Aev    := -( A + Zm1 ) * ( Apb + Zm1 ) * X / D1 / ( D1 - 1.0 );
         Aod    := Zm * ( B - Zm ) * X / D1 / ( D1 + 1.0 );
         Alo    := Bev * Ahi + Aev * Alo;
         Blo    := Bev * Bhi + Aev * Blo;
         Ahi    := Bod * Alo + Aod * Ahi;
         Bhi    := Bod * Blo + Aod * Bhi;

         IF ABS( Bhi ) < Rsmall THEN Bhi := 0.0;

         IF( Bhi <> 0.0 ) THEN
            BEGIN
               F     := Ahi / Bhi;
               Qconv := ( ABS( ( F - Fx ) / F ) < Epsz );
            END;

         Iter  := Iter + 1;

      UNTIL ( ( Iter > MaxIter ) OR Qconv ) ;

                                   (* Arrive here when convergence    *)
                                   (* achieved, or maximum iterations *)
                                   (* exceeded.                       *)

   IF ( Qswap ) THEN
      CDBeta := 1.0 - F
   ELSE
      CDBeta := F;

                                   (* Calculate precision of result *)

   IF ABS( F - Fx ) <> 0.0 THEN
      Cprec := -LogTen( ABS( F - Fx ) )
   ELSE
      Cprec := MaxPrec;

9000:  (* Error exit *)

END   (* CDBeta *);