TCDF Subroutine

public subroutine TCDF(X, Nu, Cdf)

NAME

tcdf(3f) - [M_datapac:CUMULATIVE_DISTRIBUTION] computes the cumulative
distribution function value for student's t distribution with integer
degrees of freedom NU.

SYNOPSIS

   SUBROUTINE TCDF(X,Nu,Cdf)

    REAL(kind=wp) :: X
    INTEGER       :: Nu
    REAL(kind=wp) :: Cdf

DESCRIPTION

TCDF(3f) computes the cumulative distribution function value for
Student's T distribution with integer degrees of freedom parameter
= NU. This distribution is defined for all X.

The probability density function is given in the references below.

Note the mode of internal operations is double precision.

INPUT ARGUMENTS

X      The value at which the cumulative distribution function is to
       be evaluated. X should be non-negative.

NU     The integer number of degrees of freedom. NU should be positive.

OUTPUT ARGUMENTS

CDF    The cumulative distribution function value for the Student's
       T distribution

EXAMPLES

Sample program:

program demo_tcdf
!@(#) line plotter graph of cumulative distribution function
use M_datapac, only : tcdf, plott, label
implicit none
real,allocatable  :: x(:), y(:)
integer           :: nu
integer           :: i
   call label('tcdf')
   x=[(real(i)/20.0,i=0,100,1)]
   if(allocated(y))deallocate(y)
   allocate(y(size(x)))
   nu=12
   do i=1,size(x)
      call tcdf(X(i),Nu,y(i))
   enddo
   call plott(x,y,size(x))
end program demo_tcdf

Results:

 The following is a plot of Y(I) (vertically) versus X(I) (horizontally)
                   I-----------I-----------I-----------I-----------I
  0.5000000E+01 -                                                  X
  0.4791667E+01 I                                                  X
  0.4583333E+01 I                                                  X
  0.4375000E+01 I                                                  X
  0.4166667E+01 I                                                  X
  0.3958333E+01 I                                                  X
  0.3750000E+01 -                                                  X
  0.3541667E+01 I                                                  X
  0.3333333E+01 I                                                  X
  0.3125000E+01 I                                                  X
  0.2916667E+01 I                                                 X
  0.2708333E+01 I                                                 X
  0.2500000E+01 -                                                XX
  0.2291667E+01 I                                                X
  0.2083333E+01 I                                               X
  0.1875000E+01 I                                             XX
  0.1666667E+01 I                                            XX
  0.1458333E+01 I                                         XXX
  0.1250000E+01 -                                     XXXX
  0.1041667E+01 I                                 XXXX
  0.8333335E+00 I                            XXXX
  0.6250000E+00 I                      XX XX
  0.4166670E+00 I               X XX X
  0.2083335E+00 I        XX X X
  0.0000000E+00 -  X X X
                   I-----------I-----------I-----------I-----------I
            0.5000E+00  0.6250E+00  0.7499E+00  0.8749E+00  0.9998E+00

AUTHOR

The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.

MAINTAINER

John Urban, 2022.05.31

LICENSE

CC0-1.0

REFERENCES

  • National Bureau of Standards Applied Mathematics Series 55, 1964, page 948, Formulae 26.7.3 and 26.7.4.
  • Johnson and Kotz, Continuous Univariate Distributions–2, 1970, pages 94-129.
  • Federighi, Extended Tables of the Percentage Points Of Student’S T-Distribution, Journal of the American Statistical Association, 1959, pages 683-688.
  • Owen, Handbook of Statistical Tables, 1962, pages 27-30.
  • Pearson and Hartley, Biometrika Tables for Statisticians, Volume 1, 1954, pages 132-134.

Arguments

Type IntentOptional Attributes Name
real(kind=wp) :: X
integer :: Nu
real(kind=wp) :: Cdf

Source Code

SUBROUTINE TCDF(X,Nu,Cdf)
REAL(kind=wp) :: X
INTEGER       :: Nu
REAL(kind=wp) :: Cdf
REAL(kind=wp) :: anu , cdfn , sd , z
INTEGER :: i , ievodd , imax , imin , nucut
DOUBLE PRECISION dx , dnu , pi , c , csq , s , sum , term , ai
DOUBLE PRECISION DSQRT , DATAN
DOUBLE PRECISION dconst
DOUBLE PRECISION term1 , term2 , term3
DOUBLE PRECISION dcdfn
DOUBLE PRECISION dcdf
DOUBLE PRECISION b11
DOUBLE PRECISION b21 , b22 , b23 , b24 , b25
DOUBLE PRECISION b31 , b32 , b33 , b34 , b35 , b36 , b37
DOUBLE PRECISION d1 , d3 , d5 , d7 , d9 , d11
DATA nucut/1000/
DATA pi/3.14159265358979D0/
DATA dconst/0.3989422804D0/
DATA b11/0.25D0/
DATA b21/0.01041666666667D0/
DATA b22 , b23 , b24 , b25/3.0D0 , -7.0D0 , -5.0D0 , -3.0D0/
DATA b31/0.00260416666667D0/
DATA b32 , b33 , b34 , b35 , b36 , b37/1.0D0 , -11.0D0 , 14.0D0 , &
&     6.0D0 , -3.0D0 , -15.0D0/
!
!     CHECK THE INPUT ARGUMENTS FOR ERRORS
!
      IF ( Nu<=0 ) THEN
         WRITE (G_IO,99001)
         99001 FORMAT (' ***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO TCDF(3f) IS NON-POSITIVE *****')
         WRITE (G_IO,99002) Nu
         99002 FORMAT (' ***** THE VALUE OF THE ARGUMENT IS ',I0,' *****')
         Cdf = 0.0_wp
         RETURN
      ELSE
!
!-----START POINT-----------------------------------------------------
!
         dx = X
         anu = Nu
         dnu = Nu
!
!     IF NU IS 3 THROUGH 9 AND X IS MORE THAN 3000 STANDARD DEVIATIONS BELOW THE MEAN, SET CDF = 0.0 AND RETURN.
!     IF NU IS 10 OR LARGER AND X IS MORE THAN 150 STANDARD DEVIATIONS BELOW THE MEAN, SET CDF = 0.0 AND RETURN.
!     IF NU IS 3 THROUGH 9 AND X IS MORE THAN 3000 STANDARD DEVIATIONS ABOVE THE MEAN, SET CDF = 1.0 AND RETURN.
!     IF NU IS 10 OR LARGER AND X IS MORE THAN 150 STANDARD DEVIATIONS ABOVE THE MEAN, SET CDF = 1.0 AND RETURN.
!
         IF ( Nu<=2 ) GOTO 100
         sd = SQRT(anu/(anu-2.0_wp))
         z = X/sd
         IF ( Nu>=10 .OR. z>=-3000.0_wp ) THEN
            IF ( Nu<10 .OR. z>=-150.0_wp ) THEN
               IF ( Nu<10 .AND. z>3000.0_wp ) GOTO 50
               IF ( Nu<10 .OR. z<=150.0_wp ) GOTO 100
               GOTO 50
            ENDIF
         ENDIF
         Cdf = 0.0_wp
         RETURN
 50      continue
         Cdf = 1.0_wp
         RETURN
      ENDIF
!
!     DISTINGUISH BETWEEN THE SMALL AND MODERATE
!     DEGREES OF FREEDOM CASE VERSUS THE
!     LARGE DEGREES OF FREEDOM CASE
!
 100  continue
      IF ( Nu<nucut ) THEN
!
!     TREAT THE SMALL AND MODERATE DEGREES OF FREEDOM CASE
!     METHOD UTILIZED--EXACT FINITE SUM
!     (SEE AMS 55, page 948, FORMULAE 26.7.3 AND 26.7.4).
!
         c = DSQRT(dnu/(dx*dx+dnu))
         csq = dnu/(dx*dx+dnu)
         s = dx/DSQRT(dx*dx+dnu)
         imax = Nu - 2
         ievodd = Nu - 2*(Nu/2)
         IF ( ievodd==0 ) THEN
!
            sum = 1.0D0
            term = 1.0D0
            imin = 2
         ELSE
!
            sum = c
            IF ( Nu==1 ) sum = 0.0D0
            term = c
            imin = 3
         ENDIF
!
         IF ( imin<=imax ) THEN
            DO i = imin , imax , 2
               ai = i
               term = term*((ai-1.0D0)/ai)*csq
               sum = sum + term
            ENDDO
         ENDIF
!
         sum = sum*s
         IF ( ievodd/=0 ) sum = (2.0D0/pi)*(DATAN(dx/DSQRT(dnu))+sum)
         Cdf = 0.5D0 + sum/2.0D0
         RETURN
      ELSE
!
!     TREAT THE LARGE DEGREES OF FREEDOM CASE.
!     METHOD UTILIZED--TRUNCATED ASYMPTOTIC EXPANSION
!     (SEE JOHNSON AND KOTZ, VOLUME 2, page 102, FORMULA 10;
!     SEE FEDERIGHI, page 687).
!
         CALL NORCDF(X,cdfn)
         dcdfn = cdfn
         d1 = dx
         d3 = dx**3
         d5 = dx**5
         d7 = dx**7
         d9 = dx**9
         d11 = dx**11
         term1 = b11*(d3+d1)/dnu
         term2 = b21*(b22*d7+b23*d5+b24*d3+b25*d1)/(dnu**2)
         term3 = b31*(b32*d11+b33*d9+b34*d7+b35*d5+b36*d3+b37*d1)/(dnu**3)
         dcdf = term1 + term2 + term3
         dcdf = dcdfn - (dconst*(DEXP(-dx*dx/2.0D0)))*dcdf
         Cdf = dcdf
      ENDIF
!
END SUBROUTINE TCDF