home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / MADTRB21.ZIP / MADTRB14.ZIP / CINV.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-05-17  |  6.8 KB  |  174 lines
  1. (*--------------------------------------------------------------------------*)
  2. (*                 Cinv --- Inverse chi-square routine                      *)
  3. (*--------------------------------------------------------------------------*)
  4.  
  5. FUNCTION Cinv( P, V: REAL; VAR Ifault: INTEGER ) : REAL;
  6.  
  7. (*--------------------------------------------------------------------------*)
  8. (*                                                                          *)
  9. (*      Function:  Cinv                                                     *)
  10. (*                                                                          *)
  11. (*      Purpose:  Compute inverse chi-square (percentage point)             *)
  12. (*                                                                          *)
  13. (*      Calling Sequence:                                                   *)
  14. (*                                                                          *)
  15. (*         Cp := Cinv( P, V: REAL; VAR Ifault: INTEGER ) : REAL;            *)
  16. (*                                                                          *)
  17. (*            P      --- tail probability to get percentage point for       *)
  18. (*            V      --- degrees of freedom                                 *)
  19. (*            Ifault --- error indicator                                    *)
  20. (*                       = 0:  no error                                     *)
  21. (*                       = 1:  probability too small to accurately compute  *)
  22. (*                             percentage point.                            *)
  23. (*                       = 2:  degrees of freedom given as <= 0.            *)
  24. (*                       = 3:  chi-square probability evaluation failed.    *)
  25. (*                                                                          *)
  26. (*            Cp    --- resulting chi-square percentage point               *)
  27. (*                                                                          *)
  28. (*      Calls:    GammaIn                                                   *)
  29. (*                Ninv                                                      *)
  30. (*                ALGama                                                    *)
  31. (*                                                                          *)
  32. (*      Remarks:                                                            *)
  33. (*                                                                          *)
  34. (*         This is basically AS 91 from Applied Statistics.                 *)
  35. (*         An initial approximation is improved using a Taylor series       *)
  36. (*         expansion.                                                       *)
  37. (*                                                                          *)
  38. (*--------------------------------------------------------------------------*)
  39.  
  40. CONST
  41.    E       = 1.0E-8;
  42.    Dprec   = 8;
  43.    MaxIter = 100;
  44.  
  45. VAR
  46.    XX   : REAL;
  47.    C    : REAL;
  48.    Ch   : REAL;
  49.    Q    : REAL;
  50.    P1   : REAL;
  51.    P2   : REAL;
  52.    T    : REAL;
  53.    X    : REAL;
  54.    B    : REAL;
  55.    A    : REAL;
  56.    G    : REAL;
  57.    S1   : REAL;
  58.    S2   : REAL;
  59.    S3   : REAL;
  60.    S4   : REAL;
  61.    S5   : REAL;
  62.    S6   : REAL;
  63.    Cprec: REAL;
  64.  
  65.    Iter : INTEGER;
  66.  
  67. LABEL 9000;
  68.  
  69. BEGIN (* Cinv *)
  70.                                    (* Test arguments for validity *)
  71.    Cinv   := -1.0;
  72.    Ifault := 1;
  73.  
  74.    IF ( P < E ) OR ( P > ( 1.0 - E ) ) THEN GOTO 9000;
  75.  
  76.    Ifault := 2;
  77.    IF( V <= 0.0 ) THEN GOTO 9000;
  78.  
  79.                                    (* Initialize *)
  80.    P      := 1.0 - P;
  81.    XX     := V / 2.0;
  82.    G      := ALGama( XX );
  83.    Ifault := 0;
  84.    C      := XX - 1.0;
  85.                                    (* Starting approx. for small chi-square *)
  86.  
  87.    IF( V < ( -1.24 * LN( P ) ) ) THEN
  88.       BEGIN
  89.  
  90.          Ch := Power( P * XX * EXP( G + XX * LnTwo ) , ( 1.0 / XX ) );
  91.  
  92.          IF Ch < E THEN
  93.             BEGIN
  94.                Cinv := Ch;
  95.                GOTO 9000;
  96.             END;
  97.       END
  98.    ELSE                         (* Starting approx. for v <= .32 *)
  99.       IF ( V <= 0.32 ) THEN
  100.          BEGIN
  101.  
  102.             Ch := 0.4;
  103.             A  := LN( 1.0 - P );
  104.  
  105.             REPEAT
  106.                Q  := Ch;
  107.                P1 := 1.0 + Ch * ( 4.67 + Ch );
  108.                P2 := Ch * ( 6.73 + Ch * ( 6.66 + Ch ) );
  109.                T  := -0.5 + ( 4.67 + 2.0 * Ch ) / P1 -
  110.                      ( 6.73 + Ch * ( 13.32 + 3.0 * Ch ) ) / P2;
  111.                Ch := Ch - ( 1.0 - EXP( A + G + 0.5 * Ch + C * LnTwo ) *
  112.                             P2 / P1 ) / T;
  113.             UNTIL( ABS( Q / Ch - 1.0 ) <= 0.01 );
  114.  
  115.          END
  116.       ELSE
  117.          BEGIN
  118.                                    (* Starting approx. using Wilson and *)
  119.                                    (* Hilferty estimate                 *)
  120.             X  := Ninv( P );
  121.             P1 := 2.0 / ( 9.0 * V );
  122.             Ch := V * PowerI( ( X * SQRT( P1 ) + 1.0 - P1 ) , 3 );
  123.  
  124.                                    (* Starting approx. for P ---> 1     *)
  125.  
  126.             IF ( Ch > ( 2.2 * V + 6.0 ) ) THEN
  127.                Ch := -2.0 * ( LN( 1.0 - P ) - C * LN( 0.5 * Ch ) + G );
  128.  
  129.          END;
  130.  
  131.                                    (* We have starting approximation.    *)
  132.                                    (* Begin improvement loop.            *)
  133.    REPEAT
  134.                                    (* Get probability of current approx. *)
  135.                                    (* to percentage point.               *)
  136.       Q     := Ch;
  137.       P1    := 0.5 * Ch;
  138.       P2    := P - GammaIn( P1, XX, Dprec, MaxIter, Cprec, Iter, Ifault );
  139.  
  140.       IF( Ifault <> 0 ) OR ( Iter > MaxIter ) THEN
  141.          Ifault := 3
  142.       ELSE
  143.          BEGIN
  144.                                       (* Calculate seven-term Taylor series *)
  145.  
  146.             T  := P2 * EXP( XX * LnTwo + G + P1 - C * LN( Ch ) );
  147.             B  := T / Ch;
  148.             A  := 0.5 * T - B * C;
  149.  
  150.             S1 := ( 210.0 + A * ( 140.0 + A * ( 105.0 + A * ( 84.0 + A *
  151.                   ( 70.0 + 60.0 * A ) ) ) ) ) / 420.0;
  152.             S2 := ( 420.0 + A * ( 735.0 + A * ( 966.0 + A * ( 1141.0 +
  153.                     1278.0 * A ) ) ) ) / 2520.0;
  154.             S3 := ( 210.0 + A * ( 462.0 + A * ( 707.0 + 932.0 * A ) ) )
  155.                   / 2520.0;
  156.             S4 := ( 252.0 + A * ( 672.0 + 1182.0 * A ) + C * ( 294.0 + A *
  157.                   ( 889.0 + 1740.0 * A ) ) ) / 5040.0;
  158.             S5 := ( 84.0 + 264.0 * A + C * ( 175.0 + 606.0 * A ) ) / 2520.0;
  159.             S6 := ( 120.0 + C * ( 346.0 + 127.0 * C ) ) / 5040.0;
  160.             Ch := Ch + T * ( 1.0 + 0.5 * T * S1 - B * C * ( S1 - B *
  161.                   ( S2 - B * ( S3 - B * ( S4 - B * ( S5 - B * S6 ) ) ) ) ) );
  162.  
  163.          END;
  164.  
  165.    UNTIL ( ABS( ( Q / Ch ) - 1.0 ) <= E ) OR ( Ifault <> 0 );
  166.  
  167.    IF Ifault = 0 THEN
  168.       Cinv := Ch
  169.    ELSE
  170.       Cinv := -1.0;
  171.  
  172. 9000:  ;
  173.  
  174. END   (* Cinv *);