CHSPPF Subroutine

public subroutine CHSPPF(P, Nu, Ppf)

NAME

chsppf(3f) - [M_datapac:PERCENT_POINT] compute the chi-square percent
point function

SYNOPSIS

   SUBROUTINE CHSPPF(P,Nu,Ppf)

    REAL(kind=wp) :: P
    REAL(kind=wp) :: Ppf
    INTEGER       :: Nu

DESCRIPTION

CHSPPF(3f) computes the percent point function value for the
chi-squared distribution with integer degrees of freedom parameter
= nu.

The chi-squared distribution used herein is defined for all
non-negative x, and its probability density function is given in
references 2, 3, and 4 below.

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

INPUT ARGUMENTS

P      The value (between 0.0 (inclusively) 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 for the chi-squared
       distribution

ACCURACY

(On the UNIVAC 1108, EXEC 8 System at NBS) Compared to the known NU
= 2 (exponential) results, agreement was had out to 6 significant
digits for all tested P in the range P = .001 to P = .999.  for P =
.95 And smaller, The agreement was even better--7 significant digits.
(Note that the tabulated values given in the Wilk, Gnanadesikan,
and Huyett reference below, page 20, are in error for at least the
GAMMA = 1 case-- The worst detected error was agreement to only 3
significant digits (in their 8 significant digit table) for P = .999.)

EXAMPLES

Sample program:

program demo_chsppf
use M_datapac, only : chsppf
implicit none
! call chsppf(x,y)
end program demo_chsppf

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

  • Wilk, gnanadesikan, and huyett, ‘probability plots for the gamma Distribution’, technometrics, 1962, pages 1-15, especially pages 3-5.
  • National bureau of standards applied mathematics series 55, 1964, page 257, formula 6.1.41, and pages 940-943.
  • Johnson and kotz, continuous univariate distributions–1, 1970, pages 166-206.
  • Hastings and peacock, statistical distributions–a handbook for Students and practitioners, 1975, pages 46-51.
 THIS SECTION BELOW IS LOGICALLY SEPARATE FROM THE ABOVE.
 THIS SECTION COMPUTES A CDF VALUE FOR ANY GIVEN TENTATIVE
 PERCENT POINT X VALUE AS DEFINED IN EITHER OF THE 2
 ITERATION LOOPS IN THE ABOVE CODE.

 COMPUTE T-SUB-Q AS DEFINED ON page 4 OF THE WILK, GNANADESIKAN,
 AND HUYETT REFERENCE

Arguments

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

Source Code

SUBROUTINE CHSPPF(P,Nu,Ppf)
REAL(kind=wp) :: P
REAL(kind=wp) :: Ppf
INTEGER       :: Nu
REAL(kind=wp) :: anu , dnu , gamma
INTEGER       :: icount , iloop , j , maxit

DOUBLE PRECISION :: dp , dgamma
DOUBLE PRECISION :: z , z2 , z3 , z4 , z5 , den , a , b , c , d(10) , g
DOUBLE PRECISION :: xmin0 , xmin , ai , xmax , dx , pcalc , xmid
DOUBLE PRECISION :: xlower , xupper , xdel
DOUBLE PRECISION :: sum , term , cut1 , cut2 , aj , cutoff , t
DOUBLE PRECISION :: DEXP , DLOG

DATA c/.918938533204672741D0/
DATA d(1) , d(2) , d(3) , d(4) , d(5)/ + .833333333333333333D-1 , &
     &     -.277777777777777778D-2 , +.793650793650793651D-3 ,          &
     &     -.595238095238095238D-3 , +.841750841750841751D-3/
DATA d(6) , d(7) , d(8) , d(9) , d(10)/ - .191752691752691753D-2 ,&
     &     +.641025641025641025D-2 , -.295506535947712418D-1 ,          &
     &     +.179644372368830573D0 , -.139243221690590111D1/
!
!     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 CHSPPF SUBROU&
     &TINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****')
         WRITE (G_IO,99002) P
99002    FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,       &
     &           ' *****')
         Ppf = 0.0_wp
         RETURN
      ELSEIF ( Nu<1 ) THEN
         WRITE (G_IO,99003)
99003    FORMAT (' ',                                                   &
     &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE CHSPPF SUBROU&
     &TINE IS NON-POSITIVE *****')
         WRITE (G_IO,99004) Nu
99004    FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****')
         Ppf = 0.0_wp
         RETURN
      ELSE
!
!-----START POINT-----------------------------------------------------
!
!     EXPRESS THE CHI-SQUARED DISTRIBUTION PERCENT POINT
!     FUNCTION IN TERMS OF THE EQUIVALENT GAMMA
!     DISTRIBUTION PERCENT POINT FUNCTION,
!     AND THEN EVALUATE THE LATTER.
!
         anu = Nu
         gamma = anu/2.0_wp
         dp = P
         dnu = Nu
         dgamma = dnu/2.0D0
         maxit = 10000
!
!     COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE
!     NBS APPLIED MATHEMATICS SERIES REFERENCE.
!     THIS GAMMA FUNCTION NEED BE CALCULATED ONLY ONCE.
!     IT IS USED IN THE CALCULATION OF THE CDF BASED ON
!     THE TENTATIVE VALUE OF THE PPF IN THE ITERATION.
!
         z = dgamma
         den = 1.0D0
         DO WHILE ( z<10.0D0 )
            den = den*z
            z = z + 1.0D0
         ENDDO
         z2 = z*z
         z3 = z*z2
         z4 = z2*z2
         z5 = z2*z3
         a = (z-0.5D0)*DLOG(z) - z + c
         b = d(1)/z + d(2)/z3 + d(3)/z5 + d(4)/(z2*z5) + d(5)/(z4*z5)   &
     &       + d(6)/(z*z5*z5) + d(7)/(z3*z5*z5) + d(8)/(z5*z5*z5) + d(9)&
     &       /(z2*z5*z5*z5)
         g = DEXP(a+b)/den
!
!     DETERMINE LOWER AND UPPER LIMITS ON THE DESIRED 100P
!     PERCENT POINT.
!
         iloop = 1
         xmin0 = (dp*dgamma*g)**(1.0D0/dgamma)
         xmin = xmin0
         icount = 1
      ENDIF
 100  ai = icount
      xmax = ai*xmin0
      dx = xmax
      GOTO 500
 200  xmid = (xmin+xmax)/2.0D0
!
!     NOW ITERATE BY BISECTION UNTIL THE DESIRED ACCURACY IS ACHIEVED.
!
      iloop = 2
      xlower = xmin
      xupper = xmax
      icount = 0
 300  dx = xmid
      GOTO 500
 400  Ppf = 2.0D0*xmid
      RETURN
!
!********************************************************************
!     THIS SECTION BELOW IS LOGICALLY SEPARATE FROM THE ABOVE.
!     THIS SECTION COMPUTES A CDF VALUE FOR ANY GIVEN TENTATIVE
!     PERCENT POINT X VALUE AS DEFINED IN EITHER OF THE 2
!     ITERATION LOOPS IN THE ABOVE CODE.
!
!     COMPUTE T-SUB-Q AS DEFINED ON page 4 OF THE WILK, GNANADESIKAN,
!     AND HUYETT REFERENCE
!
 500  sum = 1.0D0/dgamma
      term = 1.0D0/dgamma
      cut1 = dx - dgamma
      cut2 = dx*10000000000.0D0
      DO j = 1 , maxit
         aj = j
         term = dx*term/(dgamma+aj)
         sum = sum + term
         cutoff = cut1 + (cut2*term/sum)
         IF ( aj>cutoff ) GOTO 600
      ENDDO
      WRITE (G_IO,99005) maxit
!
99005 FORMAT (' ','*****ERROR IN INTERNAL OPERATIONS IN THE CHSPPF ',   &
     &        'SUBROUTINE--THE NUMBER OF ITERATIONS EXCEEDS ',I0)
      WRITE (G_IO,99006) P
99006 FORMAT (' ','     THE INPUT VALUE OF P     IS ',E15.8)
      WRITE (G_IO,99007) Nu
99007 FORMAT (' ','     THE INPUT VALUE OF NU    IS ',I0)
      WRITE (G_IO,99008)
99008 FORMAT (' ','     THE OUTPUT VALUE OF PPF HAS BEEN SET TO 0.0')
      Ppf = 0.0_wp
      RETURN
!
 600  t = sum
      pcalc = (dx**dgamma)*(DEXP(-dx))*t/g
      IF ( iloop==1 ) THEN
         IF ( pcalc>=dp ) GOTO 200
         xmin = xmax
         icount = icount + 1
         IF ( icount>30000 ) GOTO 200
         GOTO 100
      ELSE
         IF ( pcalc==dp ) GOTO 400
         IF ( pcalc>dp ) THEN
            xupper = xmid
            xmid = (xmid+xlower)/2.0D0
         ELSE
            xlower = xmid
            xmid = (xmid+xupper)/2.0D0
         ENDIF
         xdel = xmid - xlower
         IF ( xdel<0.0D0 ) xdel = -xdel
         icount = icount + 1
         IF ( xdel>=0.0000000001D0 .AND. icount<=100 ) GOTO 300
         GOTO 400
      ENDIF
!
END SUBROUTINE CHSPPF