[CONTACT]

[ABOUT]

[POLICY]

[ADVERTISE]

C ALGORITHM COLLECTED ALGORITHMS FROM

Found at: ftp.icm.edu.pl:70/packages/netlib/toms-2014-06-10/724

C      ALGORITHM 724, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 19, NO. 4, DECEMBER, 1993, PP. 481-483.
C
C=======================================================================
C Test driver for FINV.
C
       PROGRAM TEST
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION P(5),M(8),N(8),X(8,8)
       OPEN(11,FILE='TESTDATA',STATUS='UNKNOWN')
       P(1) = 9D-1
       P(2) = 99D-2
       P(3) = 999D-3
       P(4) = 9999D-4
       P(5) = 99999D-5
       M(1) = 1
       M(2) = 2
       M(3) = 3
       M(4) = 10
       M(5) = 20
       M(6) = 50
       M(7) = 100
       M(8) = 200
       DO 3 I = 1,8
          N(I) = M(I)
 3     CONTINUE
       DO 5 I = 1,5
          DO 10 J = 1,8
             DO 10 K = 1,8
                X(J,K) = FINV(M(J),N(K),P(I))
                if(X(J,K).LT.0.)write(*,'(A,2I4,F9.6,A)')
     *            'FINV (',M(J),N(K),P(I),' ) fails'
 10       CONTINUE
          WRITE(11,*)
          WRITE(11,'(A5,F8.5)') ' P = ',P(I)
          WRITE(11,*)
          WRITE(11,1) 'N:M',(M(J),J=1,4)
          WRITE(11,*)
          DO 15, K=1,8
             WRITE(11,2)  N(K),(X(J,K),J=1,4)
 15       CONTINUE
          WRITE(11,*)
          WRITE(11,1) 'N:M',(M(J),J=5,8)
          WRITE(11,*)
          DO 20, K=1,8
             WRITE(11,2)  N(K),(X(J,K),J=5,8)
 20       CONTINUE
          WRITE(11,*)
 5     CONTINUE
 1     FORMAT(2X,A5,1X,4(6X,I5,6X))
 2     FORMAT(2X,I5,1X,4(E17.8E3))
 
       M1 = -1
       N1 = 3
       PR = 0.9D0
       WRITE(11,'(1X,A11,I3,I3,F17.9)') 'M,N,P = ',M1,N1,PR
       FI = FINV(M1,N1,PR)
       WRITE(11,'(1X,A11,F17.9)') ' FINV = ',FI
       M1 = 3
       N1 = -1
       PR = 0.9D0
       WRITE(11,'(1X,A11,I3,I3,F17.9)') 'M,N,P = ',M1,N1,PR
       FI = FINV(M1,N1,PR)
       WRITE(11,'(1X,A11,F17.9)') ' FINV = ',FI
       M1 = 3
       N1 = 3
       PR = -0.5D0
       WRITE(11,'(1X,A11,I3,I3,F17.9)') 'M,N,P = ',M1,N1,PR
       FI = FINV(M1,N1,PR)
       WRITE(11,'(1X,A11,F17.9)') ' FINV = ',FI
       M1 = 3
       N1 = 3
       PR = 0D0
       WRITE(11,'(1X,A11,I3,I3,F17.9)') 'M,N,P = ',M1,N1,PR
       FI = FINV(M1,N1,PR)
       WRITE(11,'(1X,A11,F17.9)') ' FINV = ',FI
       CLOSE(11)
       END
C
C=======================================================================
C Software package FINV.
C
      DOUBLE PRECISION FUNCTION DBETAI(X,PIN,QIN)
C      .. Scalar Arguments ..
      DOUBLE PRECISION           X, PIN, QIN
C      .. Local Scalars ..
      DOUBLE PRECISION           C, FINSUM, P, P1, PS, Q, TERM, XB, Y
      INTEGER                    MAX1
      REAL                       SNGL
      DOUBLE PRECISION           ALNEPS, ALNSML, EPS, SML
      SAVE                       ALNEPS, ALNSML, EPS, SML
C      .. External Functions ..
      DOUBLE PRECISION           DLBETA
      EXTERNAL                   DLBETA
C      .. Intrinsic Functions ..
      INTRINSIC                  DBLE, DEXP, DINT, DLOG, DMAX1, DMIN1,
     +                           MAX1, SNGL
C      .. Executable Statements ..
      DATA                       ALNEPS/0.0D0/, ALNSML/0.0D0/,
     +                           EPS/0.0D0/, SML/0.0D0/
C**********************************************************************
C
C PURPOSE: Evaluate the incomplete beta function ratio.
C
C ARGUMENTS:
C     X   -  Upper Limit of integration.
C         X must be in the interval (0.0,1.0) inclusive.
C     PIN -  First beta distribution parameter
C         PIN must be a positive real number.
C     QIN -  Second beta distribution parameter
C         QIN must be a positive real number.
C
C EXTERNAL FUNCTIONS CALLED:
C     DLBETA-DLBETA(A,B) returns the value of the logarithm of the Beta
C         function having parameters A and B.
C
C     Note:  DBETAI is an IMSL library routine.  The authors have been
C         granted special permission to include this source code from
C         the IMSL library.
C         However, anyone wishing to use this code must first
C         purchased the library routines from IMSL.
*
C**********************************************************************
C
      IF (EPS .EQ. 0.0D0) THEN
          EPS = 1.19237D-7
          ALNEPS = DLOG(EPS)
          SML = 100.0D0*2.93941D-39
          ALNSML = DLOG(SML)
      END IF
 
      Y = X
      P = PIN
      Q = QIN
      IF (Q. LE. P .AND. X .LT. 0.8D0) GO TO 10
      IF (X .LT. 0.2D0) GO TO 10
      Y = 1.0D0 - Y
      P = QIN
      Q = PIN
 
 10   IF ((P+Q)*Y/(P+1.0D0) .LT. EPS) GO TO 70
 
      PS = Q - DINT(Q)
      IF (PS .EQ. 0.0D0) PS = 1.0D0
      XB = P*DLOG(Y) - DLBETA(PS,P) - DLOG(P)
      DBETAI = 0.0D0
      IF (XB .LT. ALNSML) GO TO 30
 
      DBETAI = DEXP(XB)
      TERM = DBETAI*P
      IF (PS .EQ. 1.0D0) GO TO 30
 
      N = MAX1(SNGL(ALNEPS/DLOG(Y)),4.0)
 
      DO 20 I=1, N
          TERM = TERM *(DBLE(I) - PS)*Y/DBLE(I)
          DBETAI = DBETAI + TERM/(P+DBLE(I))
  20  CONTINUE
 
 30   IF (Q .LE. 1.0D0) GO TO 60
          XB = P*DLOG(Y) + Q*DLOG(1.0D0-Y) - DLBETA(P,Q) - DLOG(Q)
          IB = MAX1(SNGL(XB/ALNSML),0.0)
          TERM = DEXP(XB-DBLE(IB)*ALNSML)
          C = 1.0D0/(1.0D0-Y)
          P1 = Q*C/(P+Q-1.0D0)
 
          FINSUM = 0.0D0
          N = Q
          IF (Q .EQ. DBLE(N)) N = N - 1
          DO 40 I=1, N
              IF (P1.LE.1.0D0 .AND. TERM/EPS.LE.FINSUM) GO TO 50
              TERM = (Q-DBLE(I)+1.0D0)*C*TERM/(P+Q-DBLE(I))
 
              IF (TERM .GT. 1.0D0) THEN
                  IB = IB - 1
                  TERM = TERM*SML
              END IF
 
          IF (IB .EQ. 0) FINSUM = FINSUM + TERM
 40   CONTINUE
 
 50   DBETAI = DBETAI + FINSUM
 60   IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI
      DBETAI = DMAX1(DMIN1(DBETAI,1.0D0),0.0D0)
      GO TO 9000
 
 70   DBETAI = 0.0D0
      XB = P*DLOG(DMAX1(Y,SML)) - DLOG(P) - DLBETA(P,Q)
      IF (XB.GT.ALNSML .AND. Y.NE.0.0D0) DBETAI = DEXP(XB)
      IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0 - DBETAI
 9000 RETURN
      END
 
       DOUBLE PRECISION FUNCTION DLBETA(A,B)
C      .. Parameters ..
       DOUBLE PRECISION          PI, TOL
       PARAMETER                 (PI = 3.14159265359D0, TOL = 1D-4)
C      .. Scalar Arguments ..
       DOUBLE PRECISION          A, B
C      .. Local Scalars ..
       DOUBLE PRECISION          AT, BT, DSPI, HOLD, TEMP
C      .. Intrinsic Functions ..
       INTRINSIC                 DABS, DBLE, DEXP, DLOG, DSQRT
C      .. Executable Statements ..
C**********************************************************************
C
C PURPOSE: The function returns the logarithm of the ratio:
C                      Gamma(A)*Gamma(B)/Gamma(A+B)
C          This routine was written by the authors to complement the
C          IMSL routine DBETAI which calls this function.
C
C ARGUMENTS:
C     A, B  - These are the parameters of the Beta distribution.
C          A and B are real numbers of the form (integer)/2.
C
C**********************************************************************
 
       DSPI = DSQRT(PI)
       AT = A
       BT = B
       IF (AT.LT.BT) THEN
           HOLD = AT
           AT = BT
           BT = HOLD
       ENDIF
       IF (AT .GT. 60.0D0) THEN
           IF (BT .GT. 60.0D0) THEN
               TEMP = AT*DLOG((AT+BT)/AT) + BT*DLOG((AT+BT)/BT)
               TEMP = DEXP(TEMP)*DSQRT(AT*BT/(DBLE(2)*PI*(AT+BT)))
           ELSE
               TEMP = 1.0D0
 10            IF (BT.LE.(1.0D0+TOL)) GOTO 20
                   BT = BT - 1.0D0
                   TEMP = TEMP*(AT+BT)/BT
                   GOTO 10
 20            CONTINUE
               IF (DABS(BT-0.5D0).LT.TOL)  TEMP = TEMP/DSPI
               TEMP = DLOG(TEMP) + AT*DLOG((AT+BT)/AT) - BT
               TEMP = DEXP(TEMP)*DSQRT(AT)
               IF(DABS(BT-1.0D0).LT.TOL) TEMP=TEMP*DSQRT(AT+BT)
           ENDIF
       ELSE
          TEMP = 1.0D0
 30       IF (AT.LE.(1.0D0+TOL)) GOTO 40
              AT = AT - 1.0D0
              TEMP = TEMP*(AT+BT)/AT
              GOTO 30
 40       IF (BT.LE.(1.0D0+TOL)) GOTO 50
              BT = BT - 1.0D0
              TEMP = TEMP*(AT+BT)/BT
              GOTO 40
 50       IF (DABS(AT+BT-1.5D0).LT.TOL) TEMP = TEMP*0.5D0*DSPI
          IF (DABS(AT-0.5D0).LT.TOL)  TEMP = TEMP/DSPI
          IF (DABS(BT-0.5D0).LT.TOL)  TEMP = TEMP/DSPI
       ENDIF
       DLBETA = -DLOG(TEMP)
       END

		

		
 

		
       DOUBLE PRECISION FUNCTION FINV(M,N,P)
C      .. Parameters ..
       DOUBLE PRECISION          PI
       PARAMETER               ( PI = 3.14159265359D0 )
C      .. Scalar Arguments ..
       INTEGER                   M, N
       DOUBLE PRECISION          P
C      .. Local Scalars ..
       DOUBLE PRECISION          A, B, POWER, T, F
C      .. External Functions ..
       DOUBLE PRECISION          BINV, FINIT, TINIT
       EXTERNAL                  BINV, FINIT, TINIT
C      .. Intrinsic Functions ..
       INTRINSIC                 DTAN, DBLE
C      .. Executable Statements ..
C***********************************************************************
C
C PURPOSE: This function returns percentage points for an F-distribution
C          having (M,N) degrees of freedom --- P{F(M,N) <= FINV} = P.
C
C ARGUMENTS:
C    M,N  -  The degrees of freedom of an F-distribution.
C          M and N must both be positive integers.
C    P    -  The probability for which the inverse of the F-distribution
C          is to be evaluated. 0 <= P < 1.
C    FINV -  Function value  (output)
C          The probability that a random variable with F-distribution
C          takes a value less than or equal to FINV is equal to P.
C          If M or N is negative, or if P >=1 or < 0, then FINV
C          returns with a negative value and an appropriate error
C          message is sent to standard output.
C
C EXTERNAL FUNCTIONS CALLED:
C
C    BINV  -  BINV(A,B,VAPP,P) returns the percentage points for a Beta
C          random variable, B(A,B) --- P{B(A,B) <= BINV} = P.
C          0 <= VAPP <= 1 is an initial approximation to BINV.
C          A and B are elements in {N/2 | N is a positive integer}.
C          If the series used in BINV does not converge, the
C          expansion is terminated, a negative value is returned and
C          a message noting this divergence is sent to standard
C          output.
C
C    FINIT -  FINIT(M,N,P) returns an initial approximation for the
C          point x such that a random variable, F, with F-distribution
C          having (M,N) degrees of freedom will satisfy  P{F <= x} = P.
C          0 < P < 1  and  M,N are positive integers >= 3.
C          This is the approximation due to Carter (1947).
C
C    TINIT -  TINIT(N,P) returns an initial approximation for the point
C          x such that a random variable, T, with t-distribution having
C          N degrees of freedom will satisfy  P{T >= x} = P.
C          0 < P < 1  and  N is a positive integer.
C          This is the approximation due to Goldberg and Levine (1945).
C
C***********************************************************************

		
C                                 ( FIRST, ELIMINATE BAD VALUES FOR P. )
       IF ((P.LE.0.0D0).OR.(P.GE.1.0D0))  THEN
           IF (P.EQ.0.0D0) THEN
               FINV = 0.0D0
           ELSE IF ((P.LT. 0.0D0).OR.(P.GE.1.0D0))  THEN
               FINV = -1.0D0
               PRINT *,'ERROR: PROBABILITY P IS OUT OF RANGE'
           END IF
       ELSE IF ((M.LE.0).OR.(N.LE.0)) THEN
            FINV = -1.0D0
            PRINT *,'ERROR: A PARAMETER VALUES (M OR N) IS NEGATIVE'
C                             (          THEN, CONSIDER SPECIAL CASES. )
C  CASE:  M>2 AND N>2:        (       F-DISTRIBUTION--CONVERT TO BETA. )
       ELSE IF ((M.GT.2).AND.(N.GT.2)) THEN
           F = FINIT(M,N,P)
C          (          TRANSFORM F APPROXIMATION TO BETA APPROXIMATION. )
           B = DBLE(N)/(DBLE(N)+DBLE(M)*F)
C          (                           FIND THE BETA PERCENTAGE POINT. )
           A = BINV( DBLE(N)/2.0D0, DBLE(M)/2.0D0, B, 1.0D0 - P )
C          (                     CONVERT RESULT TO F PERCENTAGE POINT. )
           FINV = ( (1.0D0-A)*DBLE(N))/(A*DBLE(M) )
C  CASE:  M=1, N=1:   (     CAUCHY DISTRIBUTION--DO DIRECT CALCULATION.)
       ELSE IF ((M.EQ.1).AND.(N.EQ.1)) THEN
           A = DTAN(P*PI/2)
           FINV = A*A
C  CASE:  M=1 OR N=1, BUT NOT BOTH: ( T DISTRIBUTION--CONVERT TO BETA. )
       ELSE IF ((M.EQ.1).OR.(N.EQ.1)) THEN
C          (            GET INITIAL APPROXIMATION FOR T, CONVERT TO F. )
           IF (M.EQ.1)  THEN
               T = TINIT(N,(1.0D0-P)/2.0D0)
               F = T*T
           ELSE
               T = TINIT(M,P/2.0D0)
               F = 1.0D0/(T*T)
           ENDIF
C          (          TRANSFORM F APPROXIMATION TO BETA APPROXIMATION. )
           B = DBLE(N)/(DBLE(N)+DBLE(M)*F)
C          (                           FIND THE BETA PERCENTAGE POINT. )
           A = BINV( DBLE(N)/2.0D0, DBLE(M)/2.0D0, B, 1.0D0 - P )
C          (                     CONVERT RESULT TO F PERCENTAGE POINT. )
           FINV = ( (1.0D0-A)*DBLE(N))/(A*DBLE(M) )
C  CASE:  M=2 OR N=2:  (    POLYNOMIAL DENSITY--DO DIRECT CALCULATION. )
       ELSE
           IF(M.EQ.2) THEN
               POWER = (1.0D0-P)**(2.0D0/DBLE(N))
               FINV = DBLE(N)*(1.0D0-POWER)/(2.0D0*POWER)
           ELSE
               POWER = P**(2.0D0/DBLE(M))
               FINV = (2.0D0*POWER)/((1.0D0-POWER)*DBLE(M))
           ENDIF
       ENDIF
       END

		
       DOUBLE PRECISION FUNCTION BINV(A,B,VAPP,P)
C      .. Parameters ..
       DOUBLE PRECISION          ERROR, ERRAPP
       PARAMETER               ( ERROR = 1.0D-8, ERRAPP = 1.0D-3 )
C      .. Scalar Arguments ..
       DOUBLE PRECISION          A, B, P, VAPP
C      .. Local Scalars ..
       DOUBLE PRECISION          BCOEFF, Q, S1, S2, SUM, T, TAIL,
     +                           VHOLD, V
       INTEGER                   I, J, K, LOOPCT
C      .. Local Array ..
       DOUBLE PRECISION          D(2:20,0:18)
C      .. External Functions ..
       DOUBLE PRECISION          BETA, DBETAI
       EXTERNAL                  BETA, DBETAI
C      .. Intrinsic Functions ..
       INTRINSIC                 DABS, DBLE, DMAX1, DMIN1
C      .. Executable Statements ..
C***********************************************************************
C
C PURPOSE: This function uses a series expansion method to compute
C          percentage points for the Beta distribution.
C
C ARGUMENTS:
C     A,B  -  The parameters of the Beta distribution.
C          A and B must both be positive real numbers of the form
C          (integer)/2.
C     VAPP -  The initial approximation to the Beta percentile.
C          VAPP is a real number between 0.0 and 1.0.
C     P    -  The probability for which the inverse Beta percentile is
C          to be evaluated.
C          P is a real number between 0.0 and 1.0.
C
C EXTERNAL FUNCTION CALLED:
C     BETA  -  BETA(A,B,X) returns the value of the density function for
C          a Beta(A,B) random variable.
C          A and B are elements in  {N/2 | N is a positive integer}.
C          0 <= X <= 1.
C
C     DBETAI - DBETAI(X,A,B) returns the probability that a random
C          variable from a Beta distribution having parameters (A,B)
C          will be less than or equal to X --- P{B(A,B) <= X} = DBETAI.
C          A and B are positive real numbers.
C
C***********************************************************************

		
       V = VAPP
       VHOLD = 0.0D0
       LOOPCT = 2
 10    IF ((DABS((V-VHOLD)/V).GE.ERRAPP).AND.(LOOPCT.NE.0)) THEN
          VHOLD = V
          LOOPCT = LOOPCT - 1
C        (          USE DBETAI TO FIND  F(V) = PROB{ BETA(A,B) <= V }. )
C        (          AND THEN COMPUTE    Q = (P - F(V))/f(V).           )
          Q = (P-DBETAI(V,A,B))/BETA(V,A,B)
C        (                       LET D(N,K) = C(N,K)*Q**(N+K-1)/(N-1)! )
          T = 1.0D0 - V
          S1 = Q*(B-1.0D0)/T
          S2 = Q*(1.0D0-A)/V
          D(2,0) = S1 + S2
          TAIL = D(2,0)*Q/2.0D0
          V = V + Q + TAIL
          K = 3
 20       IF ((DABS(TAIL/V).GT.ERROR).AND.(K.LE.20)) THEN
C           (                                   FIRST FIND  D(2,K-2).  )
             S1 = Q*(DBLE(K)-2.0D0)*S1/T
             S2 = Q*(2.0D0-DBLE(K))*S2/V
             D(2,K-2) = S1 + S2
C           (  NOW FIND  D(3,K-3), D(4,K-4), D(5,K-5), ... , D(K-1,1). )
             DO 40 I=3,K-1
                 SUM = D(2,0)*D(I-1,K-I)
                 BCOEFF = 1.0D0
                 DO 30 J = 1,K-I
                     BCOEFF = (BCOEFF*DBLE(K-I-J+1))/DBLE(J)
                     SUM = SUM + BCOEFF*D(2,J)*D(I-1,K-I-J)
 30                  CONTINUE
                 D(I,K-I) = SUM + D(I-1,K-I+1)/DBLE(I-1)
 40              CONTINUE
C           ( AND THEN COMPUTE D(K,0) AND USE IT TO EXPAND THE SERIES. )
             D(K,0) = D(2,0)*D(K-1,0) + D(K-1,1)/DBLE(K-1)
             TAIL = D(K,0)*Q/DBLE(K)
             V = V + TAIL
C           (                           CHECK FOR A DIVERGENT SERIES.  )
             IF ((V .LE. 0.0D0) .OR. (V .GE. 1.0D0))  THEN
                PRINT *,'SERIES IN BINV DIVERGES'
                V = -1.0D0
                GO TO 50
             END IF
             K = K+1
             GO TO 20
          END IF
          GO TO 10
       END IF
 50    BINV = V
       END

		
       DOUBLE PRECISION FUNCTION BETA(X,A,B)
C      .. Parameters ..
       DOUBLE PRECISION                 PI, TOL
       PARAMETER                      ( PI=3.14159265359D0, TOL=1D-4 )
C      .. Scalar Arguments ..
       DOUBLE PRECISION                 A, B, X
C      .. Local Scalars ..
       DOUBLE PRECISION                 AT, BT, HOLD, TEMP, XT, YT
C      .. Intrinsic Functions ..
       INTRINSIC                        DABS, DBLE, DEXP, DLOG, DSQRT
C      .. Executable Statements ..
C***********************************************************************
C
C PURPOSE: This function returns the value of the Beta density function:
C         BETA(X|A,B), where A and B are the parameters of the Beta.
C
C ARGUMENTS:
C     X    - The argument of the Beta density function.
C         X is a real number.
C     A,B  - The parameters of the Beta density function.
C         A and B are real numbers of the form (integer)/2.
C
C***********************************************************************

		
       AT = A
       BT = B
       XT = X
       YT = 1.0D0 - XT
       IF (AT.LT.BT) THEN
           HOLD = AT
           AT = BT
           BT = HOLD
           YT = XT
           XT = 1.0D0 - XT
       ENDIF
       IF (AT .GT. 60.0D0) THEN
           IF (BT .GT. 60.0D0) THEN
               TEMP = (AT-0.5D0)*DLOG(XT*(AT+BT-1.0D0)/(AT-1.0D0))
     +              + (BT-0.5D0)*DLOG(YT*(AT+BT-1.0D0)/(BT-1.0D0))
               TEMP = DEXP(TEMP-1.0D0)*DSQRT((AT+BT-1.0D0)/
     +                  (DBLE(2)*PI*XT*YT))
              ELSE
                  TEMP = 1.0D0
 10               IF (BT .LE. (1.0D0+TOL)) GOTO 20
                      BT = BT - 1.0D0
                      TEMP = TEMP*YT*(AT+BT)/BT
                      GOTO 10
 20               CONTINUE
                  IF (DABS(BT-0.5D0).LT.TOL)  TEMP = TEMP/DSQRT(PI*YT)
                  TEMP = DLOG(TEMP) - BT
     +                   + (AT-0.5D0)*DLOG(XT*(AT+BT-1.0D0)/(AT-1.0D0))
     +                   + BT*DLOG(AT+BT-1.0D0)
                  TEMP = DEXP(TEMP)/DSQRT(XT)
              ENDIF
          ELSE
             TEMP = 1.0D0
 30         IF (AT.LE.(1.0D0+TOL)) GOTO 40
                 AT = AT - 1.0D0
                 TEMP = TEMP*XT*(AT+BT)/AT
                 GOTO 30
 40          IF (BT.LE.(1.0D0+TOL)) GOTO 50
                 BT = BT - 1.0D0
                 TEMP = TEMP*YT*(AT+BT)/BT
                 GOTO 40
 50          IF (DABS(AT+BT-1.5D0).LT.TOL) TEMP = TEMP*0.5D0*DSQRT(PI)
             IF (DABS(AT-0.5D0).LT.TOL)  TEMP = TEMP/DSQRT(PI*XT)
             IF (DABS(BT-0.5D0).LT.TOL)  TEMP = TEMP/DSQRT(PI*YT)
          ENDIF
          BETA = TEMP
       END

		
       DOUBLE PRECISION FUNCTION ZINV(P)
C      .. Scalar Arguments ..
       DOUBLE PRECISION                  P
C      .. Local Scalars ..
       DOUBLE PRECISION                  C0, C1, C2, D1, D2, D3,
     +                                   DENOM, NUMER, PTEMP, Z
C      .. Intrinsic Functions ..
       INTRINSIC                         DLOG, DSQRT
C      .. Executable Statements ..
       DATA         C0,C1,C2 / 2.515517D0, 0.802853D0, 0.010328D0 /
       DATA         D1,D2,D3 / 1.432788D0, 0.189269D0, 0.001308D0 /
C***********************************************************************
C
C PURPOSE: This function returns an approximation for the pth percentile
C          of the standard normal distribution function.
C          This approximation is due to Hastings (1955).
C
C ARGUMENTS:
C     P   - The given percentile.
C        P is a real number between 0.0 and 1.0.
C
C***********************************************************************

		
       PTEMP = P
       IF (P .GT. 0.5D0)  PTEMP = 1.0D0 - P
       Z = DSQRT(-2.0D0*DLOG(PTEMP))
       NUMER = C0 + Z*(C1 + Z*C2)
       DENOM = 1.0D0 + Z*(D1 + Z*(D2 + Z*D3))
       Z = Z - NUMER/DENOM
       IF (P .LE. 0.5D0)  Z = -Z
       ZINV = Z
       END

		
       DOUBLE PRECISION FUNCTION FINIT(M,N,P)
C      .. Scalar Arguments ..
       DOUBLE PRECISION                  P
       INTEGER                           M, N
C      .. Local Scalars ..
       DOUBLE PRECISION                  A, B, C ,D, X
C      .. External Functions ..
       DOUBLE PRECISION                  ZINV
       EXTERNAL                          ZINV
C      .. Intrinsic Functions ..
       INTRINSIC                         DBLE, DEXP, DSQRT
C      .. Executable Statements ..
C***********************************************************************
C
C PURPOSE: This function returns an approximation to the pth percentile
C          for an F distribution with (m,n) degrees of freedom.
C          This approximation is due to Carter (1947).
C
C ARGUMENTS:
C     M,N  - The degrees of freedom of the F distribution.
C         M and N are positive integers, M >=3 and N >= 3.
C
C EXTERNAL FUNCTION CALLED:
C     ZINV -   ZINV returns X, the inverse probability for a Standard
C           Normal distribution function --- P{Z <= X} = P.
C
C***********************************************************************

		
C           (        USE ZINV TO FIND X WHEN  P = PROB{ N(0,1) <= X }. )
       X = ZINV(P)
       A = 1.0D0/(DBLE(M)-1.0D0) + 1.0D0/(DBLE(N)-1.0D0)
       B = 1.0D0/(DBLE(M)-1.0D0) - 1.0D0/(DBLE(N)-1.0D0)
       C = (X*X-3.0D0)/6.0D0
       D = X*A*DSQRT((2.0D0/A)+C) - 2.0D0*B*(C+5.0D0/6.0D0-A/3.0D0)
       FINIT = DEXP(D)
       END

		
       DOUBLE PRECISION FUNCTION TINIT(N,P)
C      .. Scalar Arguments ..
       DOUBLE PRECISION                  P
       INTEGER                           N
C      .. Local Scalars ..
       DOUBLE PRECISION                  X, XSQUAR
C      .. External Functions ..
       DOUBLE PRECISION                  ZINV
       EXTERNAL                          ZINV
C      .. Intrinsic Functions ..
       INTRINSIC                         DBLE
C      .. Executable Statements ..
C***********************************************************************
C
C PURPOSE: This function returns an approximation to the pth percentile
C          for a students-t distribution with n degrees of freedom.
C          This approximation is due to Goldberg and Levine (1945).
C
C ARGUMENTS:
C     N    -  The degree of freedom of the students-t distribution
C          N is a positive integer.
C
C EXTERNAL FUNCTION CALLED:
C     ZINV -   ZINV returns X, the inverse probability for a Standard
C          Normal distribution function --- P{Z <= X} = P.
C
C***********************************************************************

		
C         (          USE ZINV TO FIND X WHEN  P = PROB{N(0,1) <= X}. )
       X =  ZINV(P)
       XSQUAR = X*X
       TINIT = X + X*(XSQUAR+1.0D0)/(4.0D0*DBLE(N)) +  X*(XSQUAR*
     +    (XSQUAR*5.0D0 + 1.6D1) + 3.0D0)/(9.6D1*DBLE(N)*DBLE(N))
       END
C
C=======================================================================
C Output from test driver TEST.
C

		

		

		

		

		
 P =  0.90000
 
    N:M           1                2                3               10
 
      1   0.39863458E+002  0.49500000E+002  0.53593245E+002  0.60194980E+002
      2   0.85263158E+001  0.90000000E+001  0.91617902E+001  0.93915728E+001
      3   0.55383195E+001  0.54623833E+001  0.53907733E+001  0.52304113E+001
     10   0.32850153E+001  0.29244660E+001  0.27276731E+001  0.23226039E+001
     20   0.29746530E+001  0.25892541E+001  0.23800871E+001  0.19367383E+001
     50   0.28086577E+001  0.24119549E+001  0.21967298E+001  0.17291496E+001
    100   0.27563780E+001  0.23564274E+001  0.21393762E+001  0.16632251E+001
    200   0.27307304E+001  0.23292992E+001  0.21113394E+001  0.16307758E+001
 
    N:M          20               50              100              200
 
      1   0.61740292E+002  0.62688052E+002  0.63007277E+002  0.63167977E+002
      2   0.94413094E+001  0.94712356E+001  0.94812251E+001  0.94862225E+001
      3   0.51844817E+001  0.51546171E+001  0.51442595E+001  0.51390194E+001
     10   0.22007439E+001  0.21170725E+001  0.20869169E+001  0.20713517E+001
     20   0.17938433E+001  0.16896163E+001  0.16501336E+001  0.16292174E+001
     50   0.15681071E+001  0.14409422E+001  0.13884651E+001  0.13590479E+001
    100   0.14943480E+001  0.13548104E+001  0.12934390E+001  0.12570616E+001
    200   0.14574916E+001  0.13100201E+001  0.12418229E+001  0.11992434E+001
 
 
 P =  0.99000
 
    N:M           1                2                3               10
 
      1   0.40521807E+004  0.49995000E+004  0.54033520E+004  0.60558467E+004
      2   0.98502513E+002  0.99000000E+002  0.99166201E+002  0.99399196E+002
      3   0.34116222E+002  0.30816520E+002  0.29456695E+002  0.27228734E+002
     10   0.10044289E+002  0.75594322E+001  0.65523126E+001  0.48491468E+001
     20   0.80959581E+001  0.58489319E+001  0.49381934E+001  0.33681864E+001
     50   0.71705768E+001  0.50566109E+001  0.41993434E+001  0.26981394E+001
    100   0.68953010E+001  0.48239098E+001  0.39836953E+001  0.25033111E+001
    200   0.67626803E+001  0.47128548E+001  0.38807134E+001  0.24103654E+001
 
    N:M          20               50              100              200
 
      1   0.62087302E+004  0.63025172E+004  0.63341100E+004  0.63500152E+004
      2   0.99449171E+002  0.99479164E+002  0.99489163E+002  0.99494163E+002
      3   0.26689791E+002  0.26354225E+002  0.26240242E+002  0.26182916E+002
     10   0.44053948E+001  0.41154517E+001  0.40137194E+001  0.39617551E+001
     20   0.29377353E+001  0.26429545E+001  0.25353126E+001  0.24791618E+001
     50   0.22652428E+001  0.19489642E+001  0.18247532E+001  0.17567125E+001
    100   0.20666461E+001  0.17352918E+001  0.15976691E+001  0.15183428E+001
    200   0.19711304E+001  0.16294909E+001  0.14810567E+001  0.13912712E+001
 
 
 P =  0.99900
 
    N:M           1                2                3               10
 
      1   0.40528407E+006  0.49999950E+006  0.54037920E+006  0.60562097E+006
      2   0.99850025E+003  0.99900000E+003  0.99916662E+003  0.99939992E+003
      3   0.16702922E+003  0.14850000E+003  0.14110846E+003  0.12924668E+003
     10   0.21039595E+002  0.14905359E+002  0.12552745E+002  0.87538663E+001
     20   0.14818776E+002  0.99526231E+001  0.80983798E+001  0.50752462E+001
     50   0.12222106E+002  0.79564185E+001  0.63363706E+001  0.36710521E+001
    100   0.11495431E+002  0.74076811E+001  0.58568069E+001  0.32958666E+001
    200   0.11148262E+002  0.71519305E+001  0.56310496E+001  0.31203674E+001
 
    N:M          20               50              100              200
 
      1   0.62090767E+006  0.63028538E+006  0.63344433E+006  0.63503468E+006
      2   0.99944992E+003  0.99947992E+003  0.99948992E+003  0.99949492E+003
      3   0.12641777E+003  0.12466351E+003  0.12406883E+003  0.12376993E+003
     10   0.78037471E+001  0.71926834E+001  0.69801939E+001  0.68720465E+001
     20   0.42899664E+001  0.37650223E+001  0.35761909E+001  0.34783319E+001
     50   0.29506046E+001  0.24413304E+001  0.22458439E+001  0.21399504E+001
    100   0.25908800E+001  0.20755944E+001  0.18674014E+001  0.17484688E+001
    200   0.24227258E+001  0.19022348E+001  0.16824310E+001  0.15517293E+001
 
 
 P =  0.99990
 
    N:M           1                2                3               10
 
      1   0.40528473E+008  0.50000000E+008  0.54037964E+008  0.60562134E+008
      2   0.99985000E+004  0.99990000E+004  0.99991667E+004  0.99994000E+004
      3  -0.60000000E+001  0.69473833E+003 -0.20000000E+001  0.60275497E+003
     10   0.38577153E+002  0.26547867E+002  0.22037622E+002  0.14900724E+002
     20   0.23399482E+002  0.15118864E+002  0.12049800E+002  0.71805392E+001
     50   0.17878604E+002  0.11135994E+002  0.86524418E+001  0.46952294E+001
    100   0.16429538E+002  0.10113222E+002  0.77912271E+001  0.40841334E+001
    200   0.15704014E+002  0.96478196E+001  0.73708346E+001  0.37871583E+001
 
    N:M          20               50              100              200
 
      1   0.62090802E+008  0.63028572E+008  0.63344467E+008  0.63503548E+008
      2   0.99994500E+004  0.99994800E+004  0.99994900E+004  0.99994950E+004
      3   0.58929668E+003  0.58095747E+003  0.57813163E+003  0.57671147E+003
     10   0.13149531E+002  0.12031914E+002  0.11645008E+002  0.11448448E+002
     20   0.59516126E+001  0.51413101E+001  0.48523921E+001  0.47032518E+001
     50   0.36641832E+001  0.29500260E+001  0.26799064E+001  0.25346666E+001
    100   0.31036561E+001  0.24037776E+001  0.21261428E+001  0.19626450E+001
    200   0.28527336E+001  0.21550788E+001  0.18674140E+001  0.16984149E+001
 
 
 P =  0.99999
 
    N:M           1                2                3               10
 
      1   0.40528474E+010  0.50000000E+010  0.54037965E+010  0.60562134E+010
      2   0.99998500E+005  0.99999000E+005  0.99999167E+005  0.99999400E+005
      3  -0.60000000E+001  0.32301520E+004 -0.20000000E+001 -0.60000000E+000
     10   0.66427171E+002  0.45000000E+002 -0.66666667E+001  0.24621400E+002
     20   0.34265428E+002  0.21622777E+002 -0.13333333E+002  0.98057330E+001
     50   0.24146480E+002  0.14622330E+002 -0.33333333E+002  0.57926949E+001
    100   0.21661951E+002  0.12946271E+002 -0.66666667E+002  0.48844467E+001
    200   0.20012255E+002  0.12201845E+002 -0.13333333E+003  0.43125991E+001
 
    N:M          20               50              100              200
 
      1   0.62090802E+010  0.63028572E+010  0.63344467E+010  0.63503548E+010
      2   0.99999450E+005  0.99999480E+005  0.99999490E+005  0.99999495E+005
      3  -0.30000000E+000 -0.12000000E+000 -0.60000000E-001 -0.30000000E-001
     10   0.21601298E+002  0.19682063E+002  0.19019299E+002  0.18682939E+002
     20   0.80199225E+001  0.68528679E+001  0.64391713E+001  0.62261808E+001
     50   0.44237026E+001  0.34888613E+001  0.31390205E+001  0.29519431E+001
    100   0.36185942E+001  0.27301830E+001  0.23825305E+001  0.21887287E+001
    200   0.32697750E+001  0.23977060E+001  0.20437407E+001  0.18377357E+001
 
    M,N,P =  -1  3      0.900000000
     FINV =      -1.000000000
    M,N,P =   3 -1      0.900000000
     FINV =      -1.000000000
    M,N,P =   3  3     -0.500000000
     FINV =      -1.000000000
    M,N,P =   3  3      0.000000000
     FINV =       0.000000000

		
My output:

		
 
 P =  0.90000
 
    N:M           1                2                3               10
 
      1   0.39863458E+002  0.49500000E+002  0.53593245E+002  0.60194980E+002
      2   0.85263158E+001  0.90000000E+001  0.91617902E+001  0.93915728E+001
      3   0.55383195E+001  0.54623833E+001  0.53907733E+001  0.52304113E+001
     10   0.32850153E+001  0.29244660E+001  0.27276731E+001  0.23226039E+001
     20   0.29746530E+001  0.25892541E+001  0.23800871E+001  0.19367383E+001
     50   0.28086577E+001  0.24119549E+001  0.21967298E+001  0.17291496E+001
    100   0.27563780E+001  0.23564274E+001  0.21393762E+001  0.16632251E+001
    200   0.27307304E+001  0.23292992E+001  0.21113394E+001  0.16307758E+001
 
    N:M          20               50              100              200
 
      1   0.61740292E+002  0.62688052E+002  0.63007277E+002  0.63167977E+002
      2   0.94413094E+001  0.94712356E+001  0.94812251E+001  0.94862225E+001
      3   0.51844817E+001  0.51546171E+001  0.51442595E+001  0.51390194E+001
     10   0.22007439E+001  0.21170725E+001  0.20869169E+001  0.20713517E+001
     20   0.17938433E+001  0.16896163E+001  0.16501336E+001  0.16292174E+001
     50   0.15681071E+001  0.14409422E+001  0.13884651E+001  0.13590479E+001
    100   0.14943480E+001  0.13548104E+001  0.12934390E+001  0.12570616E+001
    200   0.14574916E+001  0.13100201E+001  0.12418229E+001  0.11992434E+001
 
 
 P =  0.99000
 
    N:M           1                2                3               10
 
      1   0.40521807E+004  0.49995000E+004  0.54033520E+004  0.60558467E+004
      2   0.98502513E+002  0.99000000E+002  0.99166201E+002  0.99399196E+002
      3   0.34116222E+002  0.30816520E+002  0.29456695E+002  0.27228734E+002
     10   0.10044289E+002  0.75594322E+001  0.65523126E+001  0.48491468E+001
     20   0.80959581E+001  0.58489319E+001  0.49381934E+001  0.33681864E+001
     50   0.71705768E+001  0.50566109E+001  0.41993434E+001  0.26981394E+001
    100   0.68953010E+001  0.48239098E+001  0.39836953E+001  0.25033111E+001
    200   0.67626803E+001  0.47128548E+001  0.38807134E+001  0.24103654E+001
 
    N:M          20               50              100              200
 
      1   0.62087302E+004  0.63025172E+004  0.63341100E+004  0.63500152E+004
      2   0.99449171E+002  0.99479164E+002  0.99489163E+002  0.99494163E+002
      3   0.26689791E+002  0.26354225E+002  0.26240242E+002  0.26182916E+002
     10   0.44053948E+001  0.41154517E+001  0.40137194E+001  0.39617551E+001
     20   0.29377353E+001  0.26429545E+001  0.25353126E+001  0.24791618E+001
     50   0.22652428E+001  0.19489642E+001  0.18247532E+001  0.17567125E+001
    100   0.20666461E+001  0.17352918E+001  0.15976691E+001  0.15183428E+001
    200   0.19711304E+001  0.16294909E+001  0.14810567E+001  0.13912712E+001
 
 
 P =  0.99900
 
    N:M           1                2                3               10
 
      1   0.40528407E+006  0.49999950E+006  0.54037920E+006  0.60562097E+006
      2   0.99850025E+003  0.99900000E+003  0.99916662E+003  0.99939992E+003
      3   0.16702927E+003  0.14850000E+003  0.14110846E+003  0.12924668E+003
     10   0.21039595E+002  0.14905359E+002  0.12552745E+002  0.87538663E+001
     20   0.14818776E+002  0.99526231E+001  0.80983798E+001  0.50752462E+001
     50   0.12222106E+002  0.79564185E+001  0.63363706E+001  0.36710521E+001
    100   0.11495431E+002  0.74076811E+001  0.58568069E+001  0.32958666E+001
    200   0.11148262E+002  0.71519305E+001  0.56310496E+001  0.31203674E+001
 
    N:M          20               50              100              200
 
      1   0.62090767E+006  0.63028538E+006  0.63344433E+006  0.63503468E+006
      2   0.99944992E+003  0.99947992E+003  0.99948992E+003  0.99949492E+003
      3   0.12641777E+003  0.12466351E+003  0.12406883E+003  0.12376993E+003
     10   0.78037471E+001  0.71926834E+001  0.69801939E+001  0.68720465E+001
     20   0.42899664E+001  0.37650223E+001  0.35761909E+001  0.34783319E+001
     50   0.29506046E+001  0.24413304E+001  0.22458439E+001  0.21399504E+001
    100   0.25908800E+001  0.20755944E+001  0.18674014E+001  0.17484688E+001
    200   0.24227258E+001  0.19022348E+001  0.16824310E+001  0.15517293E+001
 
 
 P =  0.99990
 
    N:M           1                2                3               10
 
      1   0.40528473E+008  0.49999999E+008  0.54037964E+008  0.60562134E+008
      2   0.99985000E+004  0.99990000E+004  0.99991667E+004  0.99994000E+004
      3  -0.60000000E+001  0.69473833E+003 -0.20000000E+001  0.60275497E+003
     10   0.38577153E+002  0.26547867E+002  0.22037622E+002  0.14900724E+002
     20   0.23399482E+002  0.15118864E+002  0.12049800E+002  0.71805392E+001
     50   0.17878604E+002  0.11135994E+002  0.86524418E+001  0.46952294E+001
    100   0.16429538E+002  0.10113222E+002  0.77912271E+001  0.40841334E+001
    200   0.15704014E+002  0.96478196E+001  0.73708353E+001  0.37871583E+001
 
    N:M          20               50              100              200
 
      1   0.62090802E+008  0.63028572E+008  0.63344467E+008  0.63503548E+008
      2   0.99994500E+004  0.99994800E+004  0.99994900E+004  0.99994950E+004
      3   0.58929668E+003  0.58095747E+003  0.57813163E+003  0.57671147E+003
     10   0.13149531E+002  0.12031914E+002  0.11645008E+002  0.11448448E+002
     20   0.59516126E+001  0.51413101E+001  0.48523921E+001  0.47032518E+001
     50   0.36641832E+001  0.29500260E+001  0.26799064E+001  0.25346666E+001
    100   0.31036561E+001  0.24037776E+001  0.21261428E+001  0.19626450E+001
    200   0.28527336E+001  0.21550788E+001  0.18674140E+001  0.16984149E+001
 
 
 P =  0.99999
 
    N:M           1                2                3               10
 
      1   0.40528474E+010  0.50000000E+010  0.54037965E+010  0.60562134E+010
      2   0.99998500E+005  0.99999000E+005  0.99999167E+005  0.99999400E+005
      3  -0.60000000E+001  0.32301520E+004 -0.20000000E+001 -0.60000000E+000
     10   0.66427171E+002  0.45000000E+002 -0.66666667E+001  0.24621400E+002
     20   0.34265428E+002  0.21622777E+002 -0.13333333E+002  0.98057330E+001
     50   0.24146480E+002  0.14622330E+002 -0.33333333E+002  0.57926949E+001
    100   0.21661951E+002  0.12946271E+002 -0.66666667E+002  0.48844467E+001
    200   0.20012255E+002  0.12201845E+002 -0.13333333E+003  0.43125991E+001
 
    N:M          20               50              100              200
 
      1   0.62090802E+010  0.63028572E+010  0.63344467E+010  0.63503548E+010
      2   0.99999450E+005  0.99999480E+005  0.99999490E+005  0.99999495E+005
      3  -0.30000000E+000 -0.12000000E+000 -0.60000000E-001 -0.30000000E-001
     10   0.21601298E+002  0.19682063E+002  0.19019299E+002  0.18682939E+002
     20   0.80199225E+001  0.68528679E+001  0.64391713E+001  0.62261808E+001
     50   0.44237026E+001  0.34888613E+001  0.31390205E+001  0.29519431E+001
    100   0.36185942E+001  0.27301830E+001  0.23825305E+001  0.21887287E+001
    200   0.32697750E+001  0.23977060E+001  0.20437407E+001  0.18377357E+001
 
    M,N,P =  -1  3      0.900000000
    M,N,P =   3 -1      0.900000000
    M,N,P =   3  3     -0.500000000
    M,N,P =   3  3      0.000000000
     FINV =       0.000000000

		
.


AD:

NEW PAGES:

[ODDNUGGET]

[GOPHER]