(*--------------------------------------------------------------------------*)
(*                BetaInv  -- Inverse Beta Distribution                     *)
(*--------------------------------------------------------------------------*)

FUNCTION BetaInv(     P, Alpha, Beta : REAL;
                      MaxIter:         INTEGER;
                      Dprec:           INTEGER;
                  VAR Iter:            INTEGER;
                  VAR Cprec:           REAL;
                  VAR Ierr:            INTEGER ) : REAL;

(*--------------------------------------------------------------------------*)
(*                                                                          *)
(*       Function:  BetaInv                                                 *)
(*                                                                          *)
(*       Purpose:  Calculates inverse Beta distribution                     *)
(*                                                                          *)
(*       Calling Sequence:                                                  *)
(*                                                                          *)
(*            B := BetaInv(      P, Alpha, Beta : REAL                      *)
(*                               MaxIter        : INTEGER;                  *)
(*                               Dprec          : INTEGER;                  *)
(*                          VAR  Iter           : INTEGER;                  *)
(*                          VAR  Cprec          : REAL;                     *)
(*                          VAR  Ierr           : INTEGER ) : REAL;         *)
(*                                                                          *)
(*                 P       --- Cumulative probability for which percentage  *)
(*                             point is to be found.                        *)
(*                 Alpha   --- First parameter of Beta distribution         *)
(*                 Beta    --- Second parameter of Beta distribution        *)
(*                 MaxIter --- Maximum iterations allowed                   *)
(*                 Dprec   --- no. of digits of precision desired           *)
(*                 Iter    --- no. of iterations actually performed         *)
(*                 Cprec   --- no. of digits of precision actually obtained *)
(*                 Ierr    --- error indicator                              *)
(*                             = 0: OK                                      *)
(*                             = 1: argument error                          *)
(*                             = 2: evaluation error during iteration       *)
(*                                                                          *)
(*                 B      --- Resulting percentage point of Beta dist.      *)
(*                                                                          *)
(*       Calls:  CDBeta  (cumulative Beta distribution)                     *)
(*                                                                          *)
(*       Remarks:                                                           *)
(*                                                                          *)
(*            Because of the relationship between the Beta                  *)
(*            distribution and the F distribution, this                     *)
(*            routine can be used to find the inverse of                    *)
(*            the F distribution.                                           *)
(*                                                                          *)
(*            The percentage point is returned as -1.0                      *)
(*            if any error occurs.                                          *)
(*                                                                          *)
(*            Newtons' method is used to search for the percentage point.   *)
(*                                                                          *)
(*            The algorithm is based upon AS110 from Applied Statistics.    *)
(*                                                                          *)
(*--------------------------------------------------------------------------*)

VAR
   Eps    : REAL;
   Xim1   : REAL;
   Xi     : REAL;
   Xip1   : REAL;
   Fim1   : REAL;
   Fi     : REAL;
   W      : REAL;
   Cmplbt : REAL;
   Adj    : REAL;
   Sq     : REAL;
   R      : REAL;
   S      : REAL;
   T      : REAL;
   G      : REAL;
   A      : REAL;
   B      : REAL;
   PP     : REAL;
   H      : REAL;
   A1     : REAL;
   B1     : REAL;
   Eprec  : REAL;

   Done   : BOOLEAN;

   Jter   : INTEGER;

LABEL 10, 30, 9000;

BEGIN (* BetaInv *)

   Ierr    := 1;
   BetaInv := P;
                                   (* Check validity of arguments *)

   IF( ( Alpha <= 0.0 ) OR ( Beta <= 0.0 ) ) THEN GOTO 9000;
   IF( ( P > 1.0 ) OR ( P < 0.0 ) )          THEN GOTO 9000;

                                   (* Check for P = 0 or 1        *)

   IF( ( P = 0.0 ) OR ( P = 1.0 ) )          THEN
      BEGIN
         Iter   := 0;
         Cprec  := MaxPrec;
         GOTO 9000;
      END;

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

   Cprec  := Dprec;

   Eps    := PowTen( -2 * Dprec );

                                   (* Flip params if needed so that *)
                                   (* P for evaluation is <= .5     *)
   IF( P > 0.5 ) THEN
      BEGIN
         A      := Beta;
         B      := Alpha;
         PP     := 1.0 - P;
      END
   ELSE
      BEGIN
         A      := Alpha;
         B      := Beta;
         PP     := P;
      END;
                                   (* Generate initial approximation.  *)
                                   (* Several different ones used,     *)
                                   (* depending upon parameter values. *)
   Ierr   := 0;

   Cmplbt := ALGama( A ) + ALGama( B ) - ALGama( A + B );
   Fi     := Ninv( 1.0 - PP );

   IF( ( A > 1.0 ) AND ( B > 1.0 ) ) THEN
      BEGIN
         R      := ( Fi * Fi - 3.0 ) / 6.0;
         S      := 1.0 / ( A + A - 1.0 );
         T      := 1.0 / ( B + B - 1.0 );
         H      := 2.0 / ( S + T );
         W      := Fi * SQRT( H + R ) / H - ( T - S ) *
                   ( R + 5.0 / 6.0 - 2.0 / ( 3.0 * H ) );
         Xi     := A / ( A + B * EXP( W + W ) );
      END
   ELSE
      BEGIN

         R      := B + B;
         T      := 1.0 / ( 9.0 * B );
         T      := R * PowerI( ( 1.0 - T + Fi * SQRT( T ) ) , 3 );

         IF( T <= 0.0 ) THEN
            Xi     := 1.0 - EXP( ( LN( ( 1.0 - PP ) * B ) + Cmplbt ) / B )
         ELSE
            BEGIN

               T      := ( 4.0 * A + R - 2.0 ) / T;

               IF( T <= 1.0 ) THEN
                  Xi := EXP( (LN( PP * A ) + Cmplbt) / PP )
               ELSE
                  Xi := 1.0 - 2.0 / ( T + 1.0 );

            END;

      END;
                                   (* Force initial estimate to *)
                                   (* reasonable range.         *)

   IF ( Xi < 0.0001 ) THEN Xi := 0.0001;
   IF ( Xi > 0.9999 ) THEN Xi := 0.9999;

                                   (* Set up Newton-Raphson loop *)

   A1     := 1.0 - A;
   B1     := 1.0 - B;
   Fim1   := 0.0;
   Sq     := 1.0;
   Xim1   := 1.0;
   Iter   := 0;
   Done   := FALSE;

                                   (* Begin Newton-Raphson loop  *)
   REPEAT

      Iter   := Iter + 1;
      Done   := Done OR ( Iter > MaxIter );

      Fi     := CDBeta( Xi, A, B, Dprec+1, MaxIter, Eprec, Jter, Ierr );

      IF( Ierr <> 0 ) THEN
         BEGIN
            Ierr := 2;
            Done := TRUE;
         END
      ELSE
         BEGIN

            Fi     := ( Fi - PP ) * EXP( Cmplbt + A1 * LN( Xi ) +
                                    B1 * LN( 1.0 - Xi ) );

            IF( ( Fi * Fim1 ) <= 0.0 ) THEN Xim1 := Sq;

            G      := 1.0;

  10:       REPEAT

               Adj    := G * Fi;
               Sq     := Adj * Adj;

               IF( Sq >= Xim1 ) THEN G := G / 3.0;

            UNTIL( Sq < Xim1 );

            Xip1     := Xi - Adj;

            IF( ( Xip1 < 0.0 ) OR ( Xip1 > 1.0 ) ) THEN
               BEGIN
                  G      := G / 3.0;
                  GOTO 10;
               END;

            IF( Xim1    <= Eps  ) THEN GOTO 30;
            IF( Fi * Fi <= Eps  ) THEN GOTO 30;

            IF( ( Xip1 = 0.0 ) OR ( Xip1 = 1.0 ) ) THEN
               BEGIN
                  G     := G / 3.0;
                  GOTO 10;
               END;

            IF( Xip1 <> Xi ) THEN
               BEGIN
                  Xi     := Xip1;
                  Fim1   := Fi;
               END
            ELSE
               Done := TRUE;

         End;

   UNTIL( Done );

30:
   BetaInv := Xi;
   IF( P > 0.5 ) THEN BetaInv := 1.0 - Xi;

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

9000:
   IF Ierr <> 0 THEN BetaInv := -1.0;

END   (* BetaInv *);