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
.