* CDFNML c c Latest revision - 03/15/90 (JRD) c c This routine computes the normal cumulative distribution c function from the error function as described in chapter 13 c of distributions in statistics - continuous univariate c distributions - 1, by Johnson and Kotz. c c Written by - Janet R. Donaldson c Statistical Engineering Division c National Bureau of Standards, Boulder, Colorado c c Creation date - December 7, 1981 c c----------------------------------------------------------------------- REAL FUNCTION CDFNML(X) c variable declarations IMPLICIT NONE c scalar arguments REAL X c external functions REAL ERF EXTERNAL ERF c variable definitions (alphabetically) c c REAL X c the percent point from the chi squared distribution. c c commence body of routine c CDFNML = 0.5E0 * (1.0E0 + ERF(X/SQRT(2.0E0))) C RETURN END * CDFT c c Latest revision - 03/15/90 (JRD) c c Purpose--This subroutine computes the cumulative distribution c function value for Student's t distribution c with integer degrees of freedom parameter = idf. c This distribution is defined for all x. c The probability density function is given c in the references below. c Input arguments--X = The value at c which the cumulative distribution c function is to be evaluated. c X should be non-negative. c --IDF = the integer number of degrees c of freedom. c IDF should be positive. c Output arguments--CDF = the single precision cumulative c distribution function value. c Output--The single precision cumulative distribution c function value (cdf) for the Student's t distribution c with degrees of freedom parameter = IDF. c Printing--None unless an input argument error condition exists. c Restrictions--IDF should be a positive integer variable. c Other datapac subroutines needed--NORCDFC [? is now CDFNML (AJA, 2016)] C Language--ANSI FORTRAN c References--National Bureau of Standards Applied Mathmatics c Series 55, 1964, page 948, formulae 26.7.3 and 26.7.4. c --Johnson and Kotz, Continuous Univariate c Distributions--2, 1970, pages 94-129. c --Federighi, Extended Tables of the c Percentage Points of Student's c t-distribution, Journal of the c American Statistical Association, c 1959, pages 683-688. c --Owen, Handbook of Statistical Tables, c 1962, pages 27-30. c --Pearson and Hartley, Biometrika Tables c for Statisticians, volume 1, 1954, c pages 132-134. c Written by--James J. Filliben c Statistical Engineering Laboratory (205.03) c National Bureau of Standards c Washington, D. C. 20234 c Original version--June 1972. c Updated --May 1974. c Updated --September 1975. c Updated --November 1975. c Updated --October 1976. c c--------------------------------------------------------------------- REAL FUNCTION CDFT (X, IDF) c Variable declarations IMPLICIT NONE c Scalar arguments REAL X INTEGER IDF c Local scalars REAL $ B11,B21,B22,B23,B24,B25,B31,B32,B33,B34,B35,B36,B37,C,CSQ,D1, $ D11,D3,D5,D7,D9,DCONST,DF,FPSPM,PI,SD,SUM,TERM,TERM1,TERM2, $ TERM3,Z INTEGER I,IDFCUT,IEVODD,IMAX,IMIN c External functions REAL CDFNML, RMACHINE c Variable definitions (alphabetically) C C REAL B11, B21, B22, B23, B24, B25 c constants used in the computations. C REAL B31, B32, B33, B34, B35, B36, B37 c constants used in the computations. C REAL DF the degrees of freedom. C REAL D1, D11, D3, D5, D7, D9 c constants used in the computations. C REAL FPSPM the floating point smallest positive magnitude. C INTEGER I an index. C INTEGER IDF the degrees of freedom. C REAL X the t statistic. DATA IDFCUT /1000/ DATA DCONST /0.3989422804E0/ DATA B11 /0.25E0/ DATA B21 /96.0E0/ C DATA B21 /0.01041666666667E0/ DATA B22, B23, B24, B25 /3.0E0,-7.0E0,-5.0E0,-3.0E0/ DATA B31 /0.00260416666667E0/ DATA B32, B33, B34, B35, B36, B37 $ /1.0E0,-11.0E0,14.0E0,6.0E0,-3.0E0,-15.0E0/ PI = ACOS(-1.) FPSPM = rmachine(2) c Check the input arguments for errors. IF (IDF.LE.0) GO TO 10 GO TO 20 10 PRINT 1000 PRINT 1010, IDF CDFT = 0.0E0 RETURN 20 CONTINUE c-----start point----------------------------------------------------- DF = IDF c If IDF is 3 through 9 and X is more than 3000 c standard deviations below the mean, c set CDFT = 0.0e0 and return. c If IDF is 10 or larger and x is more than 150 c standard deviations below the mean, c set CDFT = 0.0e0 and return. c If IDF is 3 through 9 and x is more than 3000 c standard deviations above the mean, c set CDFT = 1.0e0 and return. c If IDF is 10 or larger and x is more than 150 c standard deviations above the mean, c set CDFT = 1.0e0 and return. IF (IDF .LE. 2) GO TO 50 SD = SQRT(DF/(DF-2.0E0)) Z = X/SD IF (IDF.LT.10 .AND. Z.LT.(-3000.0E0)) GO TO 30 IF (IDF.GE.10 .AND. Z.LT.(-150.0E0)) GO TO 30 IF (IDF.LT.10 .AND. Z.GT.3000.0E0) GO TO 40 IF (IDF.GE.10 .AND. Z.GT.150.0E0) GO TO 40 GO TO 50 30 CDFT = 0.0E0 RETURN 40 CDFT = 1.0E0 RETURN 50 CONTINUE c Distinguish between the small and moderate c degrees of freedom case versus the c large degrees of freedom case IF (IDF.LT.IDFCUT) GO TO 60 GO TO 120 c Treat the small and moderate degrees of freedom case c method utilized--exact finite sum c (see AMS 55, page 948, formulae 26.7.3 and 26.7.4). 60 CONTINUE CSQ = DF/(X*X+DF) C = SQRT(CSQ) IMAX = IDF - 2 IEVODD = IDF - 2*(IDF/2) IF (IEVODD.NE.0) THEN IF (IDF.EQ.1) THEN SUM = 0.0E0 ELSE SUM = C END IF TERM = C IMIN = 3 ELSE SUM = 1.0E0 TERM = 1.0E0 IMIN = 2 END IF DO 90 I=IMIN,IMAX,2 IF (TERM.NE.0.0E0) THEN IF (LOG(TERM)+LOG((I-1.0E0)/I)+LOG(CSQ).GE.LOG(FPSPM)) THEN TERM = TERM*((I-1.0E0)/I)*CSQ SUM = SUM + TERM ELSE TERM = 0.0E0 END IF END IF 90 CONTINUE IF (SUM.EQ.0.0E0 .OR. X.EQ.0.0E0) THEN SUM = 0.0E0 ELSE IF (LOG(SUM)+LOG(ABS(X))-0.5*LOG(X*X+DF) .LT. LOG(FPSPM)) THEN SUM = 0.0E0 ELSE SUM = SUM*X/SQRT(X*X+DF) END IF END IF IF (IEVODD.EQ.0) GO TO 110 SUM = (2.0E0/PI)*(ATAN(X/SQRT(DF))+SUM) 110 CDFT = 0.5E0 + SUM/2.0E0 RETURN c Treat the large degrees of freedom case. c Method utilized--truncated asymptotic expansion c (see Johnson and Kotz, volume 2, page 102, formula 10? c see Federighi, page 687). 120 CONTINUE D1 = X D3 = X**3 D5 = X**5 D7 = X**7 D9 = X**9 D11 = X**11 TERM1 = B11*(D3+D1)/DF C TERM2 = B21*(B22*D7+B23*D5+B24*D3+B25*D1)/(DF**2) TERM2 = (B22*D7+B23*D5+B24*D3+B25*D1)/(DF**2) / B21 TERM3 = B31*(B32*D11+B33*D9+B34*D7+B35*D5+B36*D3+B37*D1)/(DF**3) CDFT = TERM1 + TERM2 + TERM3 CDFT = CDFNML(X) - (DCONST*(EXP(-X*X/2.0E0)))*CDFT RETURN C Format statements 1000 FORMAT (' ', '***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO ', $ 'THE CDFT SUBROUTINE IS NON-POSITIVE *****') 1010 FORMAT (' ', '***** THE VALUE OF THE ARGUMENT IS ', I8, ' *****') END