GAMPPF Subroutine

public subroutine GAMPPF(P, Gamma, Ppf)

NAME

gamppf(3f) - [M_datapac:PERCENT_POINT] compute the gamma percent
point function

SYNOPSIS

   SUBROUTINE GAMPPF(P,Gamma,Ppf)

    REAL(kind=wp),intent(in)  :: P
    REAL(kind=wp),intent(in)  :: Gamma
    REAL(kind=wp),intent(out) :: Ppf

DESCRIPTION

GAMPPF(3f) computes the percent point function value for the gamma
distribution with REAL tail length parameter = GAMMA.

The gamma distribution used herein has mean = GAMMA and standard
deviation = sqrt(GAMMA). This distribution is defined for all positive
X, and has the probability density function

    f(X) = (1/constant) * (X**(GAMMA-1)) * exp(-X)

where the constant = the gamma function evaluated at the value GAMMA.

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 (exclusively) and 1.0 (exclusively))
        at which the percent point function is to be evaluated.

GAMMA   The value of the tail length parameter. GAMMA should be positive.

OUTPUT ARGUMENTS

PPF     The percent point function value for the gamma distribution

EXAMPLES

Sample program:

program demo_gamppf
use M_datapac, only : gamppf
implicit none
! call gamppf(x,y)
end program demo_gamppf

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.
  • Johnson and Kotz, Continuous Univariate Distributions–1, 1970, pages 166-206.
  • Hastings and Peacock, Statistical Distributions–A Handbook for Students and Practitioners, 1975, pages 68-73.
 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), intent(in) :: P
real(kind=wp), intent(in) :: Gamma
real(kind=wp), intent(out) :: Ppf

Source Code

SUBROUTINE GAMPPF(P,Gamma,Ppf)
REAL(kind=wp),intent(in)  :: P
REAL(kind=wp),intent(in)  :: Gamma
REAL(kind=wp),intent(out) :: Ppf
INTEGER :: icount , iloop , j , maxit
!
!     ACCURACY--(ON THE UNIVAC 1108, EXEC 8 SYSTEM AT NBS)
!               COMPARED TO THE KNOWN GAMMA = 1 (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.)
!
!---------------------------------------------------------------------
!
DOUBLE PRECISION dp , dgamma
DOUBLE PRECISION z , z2 , z3 , z4 , z5 , den , a , b , c , d , 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
DIMENSION d(10)
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 GAMPPF(3f) is outside the allowable (0,1) interval *****')
         WRITE (G_IO,99007) P
         Ppf = 0.0_wp
         RETURN
      ELSEIF ( Gamma<=0.0_wp ) THEN
         WRITE (G_IO,99002)
         99002 FORMAT (' ***** FATAL ERROR--The second input argument to GAMPPF(3f) is non-positive *****')
         WRITE (G_IO,99007) Gamma
         Ppf = 0.0_wp
         RETURN
      ELSE
!
!-----START POINT-----------------------------------------------------
!
         dp = P
         dgamma = Gamma
         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 = 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,99003) maxit
!
      99003 FORMAT (' ','*****ERROR IN INTERNAL OPERATIONS in GAMPPF(3f) --The number of iterations exceeds ',I0)
      WRITE (G_IO,99004) P
      99004 FORMAT (' ','     The input value of P     is ',E15.8)
      WRITE (G_IO,99005) Gamma
      99005 FORMAT (' ','     The input value of GAMMA is ',E15.8)
      WRITE (G_IO,99006)
      99006 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
99007 FORMAT (' ***** The value of the argument is ',E15.8,' *****')
!
END SUBROUTINE GAMPPF