GAMCDF Subroutine

public subroutine GAMCDF(X, Gamma, Cdf)

NAME

gamcdf(3f) - [M_datapac:CUMULATIVE_DISTRIBUTION] compute the gamma cumulative
distribution function

SYNOPSIS

   SUBROUTINE GAMCDF(X,Gamma,Cdf)

    REAL(kind=wp),intent(in)  :: Gamma
    REAL(kind=wp),intent(in)  :: X
    REAL(kind=wp),intent(out) :: Cdf

DESCRIPTION

GAMCDF(3f) computes the cumulative distribution 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 the mode of internal operations is DOUBLE PRECISION.

ACCURACY

(On the UNIVAC 1108, EXEC 8 system at NBS)

Compared to the known GAMMA = 1 (exponential) results, agreement
was had out to 7 significant digits for all tested X.  The tested X
values covered the entire range of the distribution--from the 0.00001
percent point up to the 99.99999 percent point of the distribution.

INPUT ARGUMENTS

X      The value at which the cumulative distribution function is
       to be evaluated. X should be positive.
GAMMA  The value of the tail length parameter. GAMMA should be positive.

OUTPUT ARGUMENTS

CDF    The cumulative distribution function value for the gamma
       distribution

EXAMPLES

Sample program:

program demo_gamcdf
use M_datapac, only : gamcdf
implicit none
! call gamcdf(x,y)
end program demo_gamcdf

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.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: X
real(kind=wp), intent(in) :: Gamma
real(kind=wp), intent(out) :: Cdf

Source Code

SUBROUTINE GAMCDF(X,Gamma,Cdf)
REAL(kind=wp),intent(in)  :: Gamma
REAL(kind=wp),intent(in)  :: X
REAL(kind=wp),intent(out) :: Cdf
INTEGER :: i , maxit
DOUBLE PRECISION dx , dgamma , ai , term , sum , cut1 , cut2 , cutoff , t
DOUBLE PRECISION z , z2 , z3 , z4 , z5 , den , a , b , c , d , g
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 ( X<=0.0_wp ) THEN
         WRITE (G_IO,99001)
99001    FORMAT (' ',                                                   &
     &'***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUMENT TO THE GAMC&
     &DF SUBROUTINE IS NON-POSITIVE *****')
         WRITE (G_IO,99007) X
         Cdf = 0.0_wp
         RETURN
      ELSEIF ( Gamma<=0.0_wp ) THEN
         WRITE (G_IO,99002)
99002    FORMAT (' ',                                                   &
     &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE GAMCDF SUBROU&
     &TINE IS NON-POSITIVE *****')
         WRITE (G_IO,99007) Gamma
         Cdf = 0.0_wp
         RETURN
      ELSE
!
!-----START POINT-----------------------------------------------------
!
         dx = X
         dgamma = Gamma
         maxit = 10000
!
!     COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE
!     NBS APPLIED MATHEMATICS SERIES REFERENCE.
!
         z = dgamma
         den = 1.0D0
         DO WHILE ( z<10.0D0 )
            den = den*z
            z = z + 1
         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
!
!     COMPUTE T-SUB-Q AS DEFINED ON page 4 OF THE WILK, GNANADESIKAN,
!     AND HUYETT REFERENCE
!
         sum = 1.0D0/dgamma
         term = 1.0D0/dgamma
         cut1 = dx - dgamma
         cut2 = dx*10000000000.0D0
         DO i = 1 , maxit
            ai = i
            term = dx*term/(dgamma+ai)
            sum = sum + term
            cutoff = cut1 + (cut2*term/sum)
            IF ( ai>cutoff ) GOTO 50
         ENDDO
         WRITE (G_IO,99003) maxit
!
99003    FORMAT (' ','*****ERROR IN INTERNAL OPERATIONS IN THE GAMCDF ',&
     &           'SUBROUTINE--THE NUMBER OF ITERATIONS EXCEEDS ',I0)
         WRITE (G_IO,99004) X
99004    FORMAT (' ','     THE INPUT VALUE OF X     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 CDF HAS BEEN SET TO 1.0')
         Cdf = 1.0_wp
         RETURN
!
 50      t = sum
         Cdf = (dx**dgamma)*(DEXP(-dx))*t/g
      ENDIF
99007 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,' *****')
!
END SUBROUTINE GAMCDF