TPPF Subroutine

public subroutine TPPF(P, Nu, Ppf)

NAME

tppf(3f) - [M_datapac:PERCENT_POINT] computes the percent
point function value for the student's T distribution

SYNOPSIS

   SUBROUTINE TPPF(P,Nu,Ppf)

DESCRIPTION

tppf(3f) computes the percent point function value for the student's
t distribution with integer degrees of freedom parameter = nu.
the student's t distribution used herein is defined for all x, and
its probability density function is given in the references below.

note that the percent point function of a distribution is identically
the same as the inverse cumulative distribution function of the
distribution.

OPTIONS

 X   description of parameter
 Y   description of parameter

EXAMPLES

Sample program:

program demo_tppf
use M_datapac, only : tppf
implicit none
! call tppf(x,y)
end program demo_tppf

Results:

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 949, Formula 26.7.5.
  • Johnson and Kotz, Continuous Univariate Distributions–2, 1970, page 102, Formula 11.
  • Federighi, ‘Extended Tables of the Percentage Points of Student’s T Distribution, Journal of the American Statistical Association, 1969, pages 683-688.
  • Hastings and Peacock, Statistical Distributions–A Handbook for Students and Practitioners, 1975, pages 120-123.

NAME

tran(3f) - [M_datapac:RANDOM] a random sample of size n from the
Student's t distribution with integer degrees of freedom parameter NU.

SYNOPSIS

   SUBROUTINE TRAN(N,Nu,Iseed,X)

    INTEGER,intent(in)        :: N
    INTEGER,intent(in)        :: Nu
    INTEGER,intent(inout)     :: Iseed
    REAL(kind=wp),intent(out) :: X(:)

DESCRIPTION

This subroutine generates a random sample of size N from the Student's
T distribution with integer degrees of freedom parameter = NU.

INPUT ARGUMENTS

N     The desired integer number of random numbers to be generated.

NU    The integer degrees of freedom (parameter) for the T
      distribution. NU should be a positive integer variable.

ISEED An integer seed value. Should be set to a non-negative value to start a new sequence of values. Will be set to -1 on return to indicate the next call should continue the current random

OUTPUT ARGUMENTS

X     A vector (of dimension at least N) into which the generated
      random sample of size N from the Student's T distribution
      will be placed.

EXAMPLES

Sample program:

program demo_tran
use m_datapac, only : tran, plott, label, plotxt, sort
implicit none
integer,parameter :: n=400
real :: x(n)
integer :: iseed
integer :: nu
   call label('tran')
   nu=3
   iseed=12345
   call tran(N,Nu,Iseed,X)
   call plotxt(x,n)
   call sort(x,n,x) ! sort to show distribution
   call plotxt(x,n)
end program demo_tran

Results:

 THE FOLLOWING IS A PLOT OF X(I) (VERTICALLY) VERSUS I (HORIZONTALLY
                   I-----------I-----------I-----------I-----------I
  0.1029407E+02 -                                          X
  0.9486016E+01 I
  0.8677961E+01 I
  0.7869905E+01 I
  0.7061851E+01 I                          X
  0.6253795E+01 I                               X
  0.5445739E+01 -
  0.4637684E+01 I
  0.3829628E+01 I          X        X         X      X
  0.3021573E+01 I      X X  X XX       X X     X                 X
  0.2213517E+01 I  X    XX       XX  X    X X      X X   X
  0.1405462E+01 I   X   XXXX XXXXXXXXXXXX XX    XXX  XXX  XXXXXX
  0.5974064E+00 -  XX XXXXXXXXXXXXXXXXXXXXXXXX XXX XXXXXXXXXXXX XXX
 -0.2106485E+00 I   XXXXXXXXXXXXXX XXXXXXXXXXXXXXXXXXXXXX XXXXXXXXXX
 -0.1018704E+01 I  XX XX XX X XXX XX XXXXXXXXXXXX  XX  XXXX XX XXX X
 -0.1826759E+01 I    XXXXXX XX  XX  XX  X   XX  X     X X    XX   XX
 -0.2634815E+01 I     X           XX        X         X   X    XX
 -0.3442871E+01 I          X                   X
 -0.4250926E+01 -                        X          X          X  X
 -0.5058982E+01 I                                X              X
 -0.5867038E+01 I                     X
 -0.6675092E+01 I
 -0.7483148E+01 I
 -0.8291203E+01 I
 -0.9099259E+01 -                          X                      X
                   I-----------I-----------I-----------I-----------I
            0.1000E+01  0.1008E+03  0.2005E+03  0.3002E+03  0.4000E+03

 THE FOLLOWING IS A PLOT OF X(I) (VERTICALLY) VERSUS I (HORIZONTALLY
                   I-----------I-----------I-----------I-----------I
  0.1029407E+02 -                                                  X
  0.9486016E+01 I
  0.8677961E+01 I
  0.7869905E+01 I
  0.7061851E+01 I                                                  X
  0.6253795E+01 I                                                  X
  0.5445739E+01 -
  0.4637684E+01 I
  0.3829628E+01 I                                                 XX
  0.3021573E+01 I                                                XX
  0.2213517E+01 I                                              XXX
  0.1405462E+01 I                                         XXXXXX
  0.5974064E+00 -                              XXXXXXXXXXXX
 -0.2106485E+00 I                XXXXXXXXXXXXXXX
 -0.1018704E+01 I        XXXXXXXXX
 -0.1826759E+01 I    XXXXX
 -0.2634815E+01 I   XX
 -0.3442871E+01 I   X
 -0.4250926E+01 -   X
 -0.5058982E+01 I  X
 -0.5867038E+01 I  X
 -0.6675092E+01 I
 -0.7483148E+01 I
 -0.8291203E+01 I
 -0.9099259E+01 -  X
                   I-----------I-----------I-----------I-----------I
            0.1000E+01  0.1008E+03  0.2005E+03  0.3002E+03  0.4000E+03

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

  • Mood and Grable, Introduction to the Theory of Statistics, 1963, page 233.
  • Johnson and Kotz, Continuous Univariate Distributions–2, 1970, page 94.
  • Hastings and Peacock, Statistical Distributions–A Handbook for Students and Practitioners, 1975, page 121.

Arguments

Type IntentOptional Attributes Name
real(kind=wp) :: P
integer :: Nu
real(kind=wp) :: Ppf

Source Code

SUBROUTINE TPPF(P,Nu,Ppf)
INTEGER ipass , maxit , Nu
REAL(kind=wp) :: P , Ppf , ppfn
!
!     INPUT ARGUMENTS--P      = THE  VALUE
!                                (BETWEEN 0.0 (EXCLUSIVELY)
!                                AND 1.0 (EXCLUSIVELY))
!                                AT WHICH THE PERCENT POINT
!                                FUNCTION IS TO BE EVALUATED.
!                     --NU     = THE INTEGER NUMBER OF DEGREES
!                                OF FREEDOM.
!                                NU SHOULD BE POSITIVE.
!     OUTPUT ARGUMENTS--PPF    = THE  PERCENT
!                                POINT FUNCTION VALUE.
!     OUTPUT--THE  PERCENT POINT FUNCTION .
!             VALUE PPF FOR THE STUDENT'S T DISTRIBUTION
!             WITH DEGREES OF FREEDOM PARAMETER = NU.
!     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
!     RESTRICTIONS--NU SHOULD BE A POSITIVE INTEGER VARIABLE.
!                 --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
!                   AND 1.0 (EXCLUSIVELY).
!     OTHER DATAPAC   SUBROUTINES NEEDED--NORPPF.
!     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
!     COMMENT--FOR NU = 1 AND NU = 2, THE PERCENT POINT FUNCTION
!              FOR THE T DISTRIBUTION EXISTS IN SIMPLE CLOSED FORM
!              AND SO THE COMPUTED PERCENT POINTS ARE EXACT.
!            --FOR OTHER SMALL VALUES OF NU (NU BETWEEN 3 AND 6,
!              INCLUSIVELY), THE APPROXIMATION
!              OF THE T PERCENT POINT BY THE FORMULA
!              GIVEN IN THE REFERENCE BELOW IS AUGMENTED
!              BY 3 ITERATIONS OF NEWTON'S METHOD FOR
!              ROOT DETERMINATION.
!              THIS IMPROVES THE ACCURACY--ESPECIALLY FOR
!              VALUES OF P NEAR 0 OR 1.
!     ORIGINAL VERSION--OCTOBER   1975.
!     UPDATED         --NOVEMBER  1975.
!
!---------------------------------------------------------------------
!
      DOUBLE PRECISION pi
      DOUBLE PRECISION sqrt2
      DOUBLE PRECISION dp
      DOUBLE PRECISION dnu
      DOUBLE PRECISION term1 , term2 , term3 , term4 , term5
      DOUBLE PRECISION dppfn
      DOUBLE PRECISION dppf , dcon , darg , z , s , c
      DOUBLE PRECISION b21
      DOUBLE PRECISION b31 , b32 , b33 , b34
      DOUBLE PRECISION b41 , b42 , b43 , b44 , b45
      DOUBLE PRECISION b51 , b52 , b53 , b54 , b55 , b56
      DOUBLE PRECISION d1 , d3 , d5 , d7 , d9
      DATA pi/3.14159265358979D0/
      DATA sqrt2/1.414213562D0/
      DATA b21/0.25D0/
      DATA b31 , b32 , b33 , b34/0.01041666666667D0 , 5.0D0 , 16.0D0 ,  &
     &     3.0D0/
      DATA b41 , b42 , b43 , b44 , b45/0.00260416666667D0 , 3.0D0 ,     &
     &     19.0D0 , 17.0D0 , -15.0D0/
      DATA b51 , b52 , b53 , b54 , b55 , b56/0.00001085069444D0 ,       &
     &     79.0D0 , 776.0D0 , 1482.0D0 , -1920.0D0 , -945.0D0/
!
!     CHECK THE INPUT ARGUMENTS FOR ERRORS
!
      IF ( P<=0.0_wp .OR. P>=1.0_wp ) THEN
         WRITE (G_IO,99001)
99001    FORMAT (' ',                                                   &
     &'***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE TPPF   SUBROU&
     &TINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****')
         WRITE (G_IO,99002) P
99002    FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,       &
     &           ' *****')
         RETURN
      ELSE
!
!-----START POINT-----------------------------------------------------
!
         dnu = Nu
         dp = P
         maxit = 5
!
         IF ( Nu>=3 ) THEN
!
!     TREAT THE NU GREATER THAN OR EQUAL TO 3 CASE
!
            CALL NORPPF(P,ppfn)
            dppfn = ppfn
            d1 = dppfn
            d3 = dppfn**3
            d5 = dppfn**5
            d7 = dppfn**7
            d9 = dppfn**9
            term1 = d1
            term2 = b21*(d3+d1)/dnu
            term3 = b31*(b32*d5+b33*d3+b34*d1)/(dnu**2)
            term4 = b41*(b42*d7+b43*d5+b44*d3+b45*d1)/(dnu**3)
            term5 = b51*(b52*d9+b53*d7+b54*d5+b55*d3+b56*d1)/(dnu**4)
            dppf = term1 + term2 + term3 + term4 + term5
            Ppf = dppf
            IF ( Nu>=7 ) RETURN
            IF ( Nu==3 ) THEN
!
!     AUGMENT THE RESULTS FOR THE NU = 3 CASE
!
               dcon = pi*(dp-0.5D0)
               darg = dppf/DSQRT(dnu)
               z = DATAN(darg)
               DO ipass = 1 , maxit
                  s = DSIN(z)
                  c = DCOS(z)
                  z = z - (z+s*c-dcon)/(2.0D0*c*c)
               ENDDO
               Ppf = DSQRT(dnu)*s/c
               RETURN
            ELSEIF ( Nu==4 ) THEN
!
!     AUGMENT THE RESULTS FOR THE NU = 4 CASE
!
               dcon = 2.0D0*(dp-0.5D0)
               darg = dppf/DSQRT(dnu)
               z = DATAN(darg)
               DO ipass = 1 , maxit
                  s = DSIN(z)
                  c = DCOS(z)
                  z = z - ((1.0D0+0.5D0*c*c)*s-dcon)/(1.5D0*c*c*c)
               ENDDO
               Ppf = DSQRT(dnu)*s/c
               RETURN
            ELSEIF ( Nu==5 ) THEN
!
!     AUGMENT THE RESULTS FOR THE NU = 5 CASE
!
               dcon = pi*(dp-0.5D0)
               darg = dppf/DSQRT(dnu)
               z = DATAN(darg)
               DO ipass = 1 , maxit
                  s = DSIN(z)
                  c = DCOS(z)
                  z = z - (z+(c+(2.0D0/3.0D0)*c*c*c)*s-dcon)            &
     &                /((8.0D0/3.0D0)*c**4)
               ENDDO
               Ppf = DSQRT(dnu)*s/c
               RETURN
            ELSEIF ( Nu==6 ) THEN
!
!     AUGMENT THE RESULTS FOR THE NU = 6 CASE
!
               dcon = 2.0D0*(dp-0.5D0)
               darg = dppf/DSQRT(dnu)
               z = DATAN(darg)
               DO ipass = 1 , maxit
                  s = DSIN(z)
                  c = DCOS(z)
                  z = z - ((1.0D0+0.5D0*c*c+0.375D0*c**4)*s-dcon)       &
     &                /((15.0D0/8.0D0)*c**5)
               ENDDO
               Ppf = DSQRT(dnu)*s/c
               GOTO 99999
            ENDIF
         ELSEIF ( Nu==1 ) THEN
!
!     TREAT THE NU = 1 (CAUCHY) CASE
!
            darg = pi*dp
            Ppf = -DCOS(darg)/DSIN(darg)
            RETURN
         ELSEIF ( Nu==2 ) THEN
!
!     TREAT THE NU = 2 CASE
!
            term1 = sqrt2/2.0D0
            term2 = 2.0D0*dp - 1.0D0
            term3 = DSQRT(dp*(1.0D0-dp))
            Ppf = term1*term2/term3
            RETURN
         ELSE
            WRITE (G_IO,99003)
99003       FORMAT (' ','INTERNAL ERROR IN TPPF SUBROUTINE')
            Ppf = 0.0_wp
            RETURN
         ENDIF
      ENDIF
      RETURN
!
99999 END SUBROUTINE TPPF