home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / MADTRB21.ZIP / MADTRB14.ZIP / BETAINV.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-05-17  |  9.1 KB  |  261 lines
  1. (*--------------------------------------------------------------------------*)
  2. (*                BetaInv  -- Inverse Beta Distribution                     *)
  3. (*--------------------------------------------------------------------------*)
  4.  
  5. FUNCTION BetaInv(     P, Alpha, Beta : REAL;
  6.                       MaxIter:         INTEGER;
  7.                       Dprec:           INTEGER;
  8.                   VAR Iter:            INTEGER;
  9.                   VAR Cprec:           REAL;
  10.                   VAR Ierr:            INTEGER ) : REAL;
  11.  
  12. (*--------------------------------------------------------------------------*)
  13. (*                                                                          *)
  14. (*       Function:  BetaInv                                                 *)
  15. (*                                                                          *)
  16. (*       Purpose:  Calculates inverse Beta distribution                     *)
  17. (*                                                                          *)
  18. (*       Calling Sequence:                                                  *)
  19. (*                                                                          *)
  20. (*            B := BetaInv(      P, Alpha, Beta : REAL                      *)
  21. (*                               MaxIter        : INTEGER;                  *)
  22. (*                               Dprec          : INTEGER;                  *)
  23. (*                          VAR  Iter           : INTEGER;                  *)
  24. (*                          VAR  Cprec          : REAL;                     *)
  25. (*                          VAR  Ierr           : INTEGER ) : REAL;         *)
  26. (*                                                                          *)
  27. (*                 P       --- Cumulative probability for which percentage  *)
  28. (*                             point is to be found.                        *)
  29. (*                 Alpha   --- First parameter of Beta distribution         *)
  30. (*                 Beta    --- Second parameter of Beta distribution        *)
  31. (*                 MaxIter --- Maximum iterations allowed                   *)
  32. (*                 Dprec   --- no. of digits of precision desired           *)
  33. (*                 Iter    --- no. of iterations actually performed         *)
  34. (*                 Cprec   --- no. of digits of precision actually obtained *)
  35. (*                 Ierr    --- error indicator                              *)
  36. (*                             = 0: OK                                      *)
  37. (*                             = 1: argument error                          *)
  38. (*                             = 2: evaluation error during iteration       *)
  39. (*                                                                          *)
  40. (*                 B      --- Resulting percentage point of Beta dist.      *)
  41. (*                                                                          *)
  42. (*       Calls:  CDBeta  (cumulative Beta distribution)                     *)
  43. (*                                                                          *)
  44. (*       Remarks:                                                           *)
  45. (*                                                                          *)
  46. (*            Because of the relationship between the Beta                  *)
  47. (*            distribution and the F distribution, this                     *)
  48. (*            routine can be used to find the inverse of                    *)
  49. (*            the F distribution.                                           *)
  50. (*                                                                          *)
  51. (*            The percentage point is returned as -1.0                      *)
  52. (*            if any error occurs.                                          *)
  53. (*                                                                          *)
  54. (*            Newtons' method is used to search for the percentage point.   *)
  55. (*                                                                          *)
  56. (*            The algorithm is based upon AS110 from Applied Statistics.    *)
  57. (*                                                                          *)
  58. (*--------------------------------------------------------------------------*)
  59.  
  60. VAR
  61.    Eps    : REAL;
  62.    Xim1   : REAL;
  63.    Xi     : REAL;
  64.    Xip1   : REAL;
  65.    Fim1   : REAL;
  66.    Fi     : REAL;
  67.    W      : REAL;
  68.    Cmplbt : REAL;
  69.    Adj    : REAL;
  70.    Sq     : REAL;
  71.    R      : REAL;
  72.    S      : REAL;
  73.    T      : REAL;
  74.    G      : REAL;
  75.    A      : REAL;
  76.    B      : REAL;
  77.    PP     : REAL;
  78.    H      : REAL;
  79.    A1     : REAL;
  80.    B1     : REAL;
  81.    Eprec  : REAL;
  82.  
  83.    Done   : BOOLEAN;
  84.  
  85.    Jter   : INTEGER;
  86.  
  87. LABEL 10, 30, 9000;
  88.  
  89. BEGIN (* BetaInv *)
  90.  
  91.    Ierr    := 1;
  92.    BetaInv := P;
  93.                                    (* Check validity of arguments *)
  94.  
  95.    IF( ( Alpha <= 0.0 ) OR ( Beta <= 0.0 ) ) THEN GOTO 9000;
  96.    IF( ( P > 1.0 ) OR ( P < 0.0 ) )          THEN GOTO 9000;
  97.  
  98.                                    (* Check for P = 0 or 1        *)
  99.  
  100.    IF( ( P = 0.0 ) OR ( P = 1.0 ) )          THEN
  101.       BEGIN
  102.          Iter   := 0;
  103.          Cprec  := MaxPrec;
  104.          GOTO 9000;
  105.       END;
  106.  
  107.                                   (* Set precision *)
  108.    IF Dprec > MaxPrec THEN
  109.       Dprec := MaxPrec
  110.    ELSE IF Dprec <= 0 THEN
  111.       Dprec := 1;
  112.  
  113.    Cprec  := Dprec;
  114.  
  115.    Eps    := PowTen( -2 * Dprec );
  116.  
  117.                                    (* Flip params if needed so that *)
  118.                                    (* P for evaluation is <= .5     *)
  119.    IF( P > 0.5 ) THEN
  120.       BEGIN
  121.          A      := Beta;
  122.          B      := Alpha;
  123.          PP     := 1.0 - P;
  124.       END
  125.    ELSE
  126.       BEGIN
  127.          A      := Alpha;
  128.          B      := Beta;
  129.          PP     := P;
  130.       END;
  131.                                    (* Generate initial approximation.  *)
  132.                                    (* Several different ones used,     *)
  133.                                    (* depending upon parameter values. *)
  134.    Ierr   := 0;
  135.  
  136.    Cmplbt := ALGama( A ) + ALGama( B ) - ALGama( A + B );
  137.    Fi     := Ninv( 1.0 - PP );
  138.  
  139.    IF( ( A > 1.0 ) AND ( B > 1.0 ) ) THEN
  140.       BEGIN
  141.          R      := ( Fi * Fi - 3.0 ) / 6.0;
  142.          S      := 1.0 / ( A + A - 1.0 );
  143.          T      := 1.0 / ( B + B - 1.0 );
  144.          H      := 2.0 / ( S + T );
  145.          W      := Fi * SQRT( H + R ) / H - ( T - S ) *
  146.                    ( R + 5.0 / 6.0 - 2.0 / ( 3.0 * H ) );
  147.          Xi     := A / ( A + B * EXP( W + W ) );
  148.       END
  149.    ELSE
  150.       BEGIN
  151.  
  152.          R      := B + B;
  153.          T      := 1.0 / ( 9.0 * B );
  154.          T      := R * PowerI( ( 1.0 - T + Fi * SQRT( T ) ) , 3 );
  155.  
  156.          IF( T <= 0.0 ) THEN
  157.             Xi     := 1.0 - EXP( ( LN( ( 1.0 - PP ) * B ) + Cmplbt ) / B )
  158.          ELSE
  159.             BEGIN
  160.  
  161.                T      := ( 4.0 * A + R - 2.0 ) / T;
  162.  
  163.                IF( T <= 1.0 ) THEN
  164.                   Xi := EXP( (LN( PP * A ) + Cmplbt) / PP )
  165.                ELSE
  166.                   Xi := 1.0 - 2.0 / ( T + 1.0 );
  167.  
  168.             END;
  169.  
  170.       END;
  171.                                    (* Force initial estimate to *)
  172.                                    (* reasonable range.         *)
  173.  
  174.    IF ( Xi < 0.0001 ) THEN Xi := 0.0001;
  175.    IF ( Xi > 0.9999 ) THEN Xi := 0.9999;
  176.  
  177.                                    (* Set up Newton-Raphson loop *)
  178.  
  179.    A1     := 1.0 - A;
  180.    B1     := 1.0 - B;
  181.    Fim1   := 0.0;
  182.    Sq     := 1.0;
  183.    Xim1   := 1.0;
  184.    Iter   := 0;
  185.    Done   := FALSE;
  186.  
  187.                                    (* Begin Newton-Raphson loop  *)
  188.    REPEAT
  189.  
  190.       Iter   := Iter + 1;
  191.       Done   := Done OR ( Iter > MaxIter );
  192.  
  193.       Fi     := CDBeta( Xi, A, B, Dprec+1, MaxIter, Eprec, Jter, Ierr );
  194.  
  195.       IF( Ierr <> 0 ) THEN
  196.          BEGIN
  197.             Ierr := 2;
  198.             Done := TRUE;
  199.          END
  200.       ELSE
  201.          BEGIN
  202.  
  203.             Fi     := ( Fi - PP ) * EXP( Cmplbt + A1 * LN( Xi ) +
  204.                                     B1 * LN( 1.0 - Xi ) );
  205.  
  206.             IF( ( Fi * Fim1 ) <= 0.0 ) THEN Xim1 := Sq;
  207.  
  208.             G      := 1.0;
  209.  
  210.   10:       REPEAT
  211.  
  212.                Adj    := G * Fi;
  213.                Sq     := Adj * Adj;
  214.  
  215.                IF( Sq >= Xim1 ) THEN G := G / 3.0;
  216.  
  217.             UNTIL( Sq < Xim1 );
  218.  
  219.             Xip1     := Xi - Adj;
  220.  
  221.             IF( ( Xip1 < 0.0 ) OR ( Xip1 > 1.0 ) ) THEN
  222.                BEGIN
  223.                   G      := G / 3.0;
  224.                   GOTO 10;
  225.                END;
  226.  
  227.             IF( Xim1    <= Eps  ) THEN GOTO 30;
  228.             IF( Fi * Fi <= Eps  ) THEN GOTO 30;
  229.  
  230.             IF( ( Xip1 = 0.0 ) OR ( Xip1 = 1.0 ) ) THEN
  231.                BEGIN
  232.                   G     := G / 3.0;
  233.                   GOTO 10;
  234.                END;
  235.  
  236.             IF( Xip1 <> Xi ) THEN
  237.                BEGIN
  238.                   Xi     := Xip1;
  239.                   Fim1   := Fi;
  240.                END
  241.             ELSE
  242.                Done := TRUE;
  243.  
  244.          End;
  245.  
  246.    UNTIL( Done );
  247.  
  248. 30:
  249.    BetaInv := Xi;
  250.    IF( P > 0.5 ) THEN BetaInv := 1.0 - Xi;
  251.  
  252.    IF ABS( Xi - Xim1 ) <> 0.0 THEN
  253.       Cprec  := -LogTen( ABS( Xi - Xim1 ) )
  254.    ELSE
  255.       Cprec  := MaxPrec;
  256.  
  257. 9000:
  258.    IF Ierr <> 0 THEN BetaInv := -1.0;
  259.  
  260. END   (* BetaInv *);
  261.