NBPPF Subroutine

public subroutine NBPPF(P, Ppar, N, Ppf)

NAME

nbppf(3f) - [M_datapac:PERCENT_POINT] compute the negative binomial
percent point function

SYNOPSIS

   SUBROUTINE NBPPF(P,Ppar,N,Ppf)

DESCRIPTION

nbppf(3f) computes the percent point function value at the precision
precision value p for the negative binomial distribution with precision
precision 'bernoulli probability' parameter = ppar, and integer
'number of successes in bernoulli trials' parameter = n.

the negative binomial distribution used herein has mean =
n*(1-ppar)/ppar and standard deviation = sqrt(n*(1-ppar)/(ppar*ppar))).
this distribution is defined for all non-negative integer x--x = 0,
1, 2, ... .

this distribution has the probability function

    f(x) = c(n+x-1,n) * ppar**n * (1-ppar)**x.

where c(n+x-1,n) is the combinatorial function
equaling the number of combinations of n+x-1 items
taken n at a time.

the negative binomial distribution is the distribution of the number
of failures before obtaining n successes in an indefinite sequence of
bernoulli (0,1) trials where the probability of success in a precision
trial = ppar.

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_nbppf
use M_datapac, only : nbppf
implicit none
! call nbppf(x,y)
end program demo_nbppf

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

  • JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS, 1969, pages 122-142, ESPECIALLY page 127, FORMULA 22.
  • HASTINGS AND PEACOCK, STATISTICAL DISTRIBUTIONS–A HANDBOOK FOR STUDENTS AND PRACTITIONERS, 1975, pages 92-95.
  • NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS SERIES 55, 1964, page 929.
  • FELLER, AN INTRODUCTION TO PROBABILITY THEORY AND ITS APPLICATIONS, VOLUME 1, EDITION 2, 1957, pages 155-157, 210.
  • KENDALL AND STUART, THE ADVANCED THEORY OF STATISTICS, VOLUME 1, EDITION 2, 1963, pages 130-131.
  • WILLIAMSON AND BRETHERTON, TABLES OF THE NEGATIVE BINOMIAL PROBABILITY DISTRIBUTION, 1963.
  • OWEN, HANDBOOK OF STATISTICAL TABLES, 1962, page 304.

Arguments

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

Source Code

      SUBROUTINE NBPPF(P,Ppar,N,Ppf)
REAL(kind=wp) :: amean , an , arcsh , arg , e , P , p0 , p1 , p2 , pf0 ,      &
     &     Ppar , Ppf , sd , sinh , x0 , x1 , x2 , ymean , yppf , ysd
REAL(kind=wp) :: zppf
INTEGER :: i , isd , ix0 , ix0p1 , ix1 , ix2 , N
!
!     INPUT ARGUMENTS--P      = THE  VALUE
!                                (BETWEEN 0.0 (INCLUSIVELY)
!                                AND 1.0 (EXCLUSIVELY))
!                                AT WHICH THE PERCENT POINT
!                                FUNCTION IS TO BE EVALUATED.
!                     --PPAR   = THE  VALUE
!                                OF THE 'BERNOULLI PROBABILITY'
!                                PARAMETER FOR THE NEGATIVE BINOMIAL
!                                DISTRIBUTION.
!                                PPAR SHOULD BE BETWEEN
!                                0.0 (EXCLUSIVELY) AND
!                                1.0 (EXCLUSIVELY).
!                     --N      = THE INTEGER VALUE
!                                OF THE 'NUMBER OF SUCCESSES
!                                IN BERNOULLI TRIALS' PARAMETER.
!                                N SHOULD BE A POSITIVE INTEGER.
!     OUTPUT ARGUMENTS--PPF    = THE  PERCENT
!                                POINT FUNCTION VALUE.
!     OUTPUT--THE  PERCENT POINT  .
!             FUNCTION VALUE PPF
!             FOR THE NEGATIVE BINOMIAL DISTRIBUTION
!             WITH 'BERNOULLI PROBABILITY' PARAMETER = PPAR
!             AND 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS'
!             PARAMETER = N.
!     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
!     RESTRICTIONS--PPAR SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
!                   AND 1.0 (EXCLUSIVELY).
!                 --N SHOULD BE A POSITIVE INTEGER.
!                 --P SHOULD BE BETWEEN 0.0 (INCLUSIVELY)
!                   AND 1.0 (EXCLUSIVELY).
!     OTHER DATAPAC   SUBROUTINES NEEDED--NORPPF, NBCDF.
!     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, EXP, LOG.
!     MODE OF INTERNAL OPERATIONS-- AND DOUBLE PRECISION.
!     COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT
!              FROM THIS DISCRETE DISTRIBUTION
!              PERCENT POINT FUNCTION
!              SUBROUTINE MUST NECESSARILY BE A
!              DISCRETE INTEGER VALUE,
!              THE OUTPUT VARIABLE PPF IS SINGLE
!              PRECISION IN MODE.
!              PPF HAS BEEN SPECIFIED AS SINGLE
!              PRECISION SO AS TO CONFORM WITH THE DATAPAC
!              CONVENTION THAT ALL OUTPUT VARIABLES FROM ALL
!              DATAPAC SUBROUTINES ARE .
!              THIS CONVENTION IS BASED ON THE BELIEF THAT
!              1) A MIXTURE OF MODES (FLOATING POINT
!              VERSUS INTEGER) IS INCONSISTENT AND
!              AN UNNECESSARY COMPLICATION
!              IN A DATA ANALYSIS; AND
!              2) FLOATING POINT MACHINE ARITHMETIC
!              (AS OPPOSED TO INTEGER ARITHMETIC)
!              IS THE MORE NATURAL MODE FOR DOING
!              DATA ANALYSIS.
!     ORIGINAL VERSION--NOVEMBER  1975.
!
!---------------------------------------------------------------------
!
      DOUBLE PRECISION dppar
!
!     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 NBPPF  SUBROU&
     &TINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****')
         WRITE (G_IO,99019) P
         Ppf = 0.0_wp
         RETURN
      ELSE
         IF ( Ppar<=0.0_wp .OR. Ppar>=1.0_wp ) THEN
            WRITE (G_IO,99002)
99002       FORMAT (' ',                                                &
     &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE NBPPF  SUBROU&
     &TINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****')
            WRITE (G_IO,99019) Ppar
            Ppf = 0.0_wp
            RETURN
         ELSE
            IF ( N<1 ) THEN
               WRITE (G_IO,99003)
99003          FORMAT (' ',                                             &
     &'***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE NBPPF  SUBROU&
     &TINE IS NON-POSITIVE *****')
               WRITE (G_IO,99004) N
99004          FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,    &
     &                 ' *****')
               Ppf = 0.0_wp
               RETURN
            ELSE
!
!-----START POINT-----------------------------------------------------
!
               an = N
               dppar = Ppar
               Ppf = 0.0_wp
               ix0 = 0
               ix1 = 0
               ix2 = 0
               p0 = 0.0_wp
               p1 = 0.0_wp
               p2 = 0.0_wp
!
!     TREAT CERTAIN SPECIAL CASES IMMEDIATELY--
!     1) P = 0.0
!     2) P = 0.5 AND PPAR = 0.5
!     3) PPF = 0
!
               IF ( P/=0.0_wp ) THEN
                  IF ( P==0.5_wp .AND. Ppar==0.5_wp ) THEN
                     Ppf = N - 1
                     RETURN
                  ELSE
                     pf0 = dppar**N
                     IF ( P>pf0 ) THEN
!
!     DETERMINE AN INITIAL APPROXIMATION TO THE NEGATIVE BINOMIAL
!     PERCENT POINT BY USE OF THE HYPERBOLIC ARCSIN
!     TRANSFORMATION OF THE NEGATIVE BINOMIAL
!     TO APPROXIMATE NORMALITY.
!     (SEE JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS,
!     page 127, FORMULA 22).
!
                        amean = an*(1.0_wp-Ppar)/Ppar
                        sd = SQRT(an*(1.0_wp-Ppar)/(Ppar*Ppar))
                        arg = SQRT((amean+0.375_wp)/(an-0.75_wp))
                        arcsh = LOG(arg+SQRT(arg*arg+1.0_wp))
                        ymean = (SQRT(an-0.5_wp))*arcsh
                        ysd = 0.5_wp
                        CALL NORPPF(P,zppf)
                        yppf = ymean + zppf*ysd
                        arg = yppf/SQRT(an-0.5_wp)
                        e = EXP(arg)
                        sinh = (e-1.0_wp/e)/2.0_wp
                        x2 = -0.375_wp + (an-0.75_wp)*sinh*sinh
                        x2 = x2 + 0.5_wp
                        ix2 = x2
!
!     CHECK AND MODIFY (IF NECESSARY) THIS INITIAL
!     ESTIMATE OF THE PERCENT POINT
!     TO ASSURE THAT IT BE NON-NEGATIVE.
!
                        IF ( ix2<0 ) ix2 = 0
!
!     DETERMINE UPPER AND LOWER BOUNDS ON THE DESIRED
!     PERCENT POINT BY ITERATING OUT (BOTH BELOW AND ABOVE)
!     FROM THE ORIGINAL APPROXIMATION AT STEPS
!     OF 1 STANDARD DEVIATION.
!     THE RESULTING BOUNDS WILL BE AT MOST
!     1 STANDARD DEVIATION APART.
!
                        ix0 = 0
                        ix1 = huge(0) ! = 10**10
                        isd = sd + 1.0_wp
                        x2 = ix2
                        CALL NBCDF(x2,Ppar,N,p2)
!
                        IF ( p2<P ) THEN
!
                           ix0 = ix2
                           DO i = 1 , 100000
                              ix2 = ix0 + isd
                              IF ( ix2>=ix1 ) GOTO 100
                              x2 = ix2
                              CALL NBCDF(x2,Ppar,N,p2)
                              IF ( p2>=P ) GOTO 20
                              ix0 = ix2
                           ENDDO
                           WRITE (G_IO,99020)
                           WRITE (G_IO,99005)
!
99005                      FORMAT (' ',                                 &
     &                     'NO UPPER BOUND FOUND AFTER 10**7 ITERATIONS'&
     &                     )
                        ELSE
!
                           ix1 = ix2
                           DO i = 1 , 100000
                              ix2 = ix1 - isd
                              IF ( ix2<=ix0 ) GOTO 100
                              x2 = ix2
                              CALL NBCDF(x2,Ppar,N,p2)
                              IF ( p2<P ) GOTO 50
                              ix1 = ix2
                           ENDDO
                           WRITE (G_IO,99020)
                           WRITE (G_IO,99006)
99006                      FORMAT (' ',                                 &
     &                     'NO LOWER BOUND FOUND AFTER 10**7 ITERATIONS'&
     &                     )
                        ENDIF
                        GOTO 200
                     ENDIF
                  ENDIF
               ENDIF
               Ppf = 0.0_wp
               RETURN
            ENDIF
 20         ix1 = ix2
            GOTO 100
         ENDIF
 50      ix0 = ix2
      ENDIF
!
 100  IF ( ix0==ix1 ) THEN
         IF ( ix0==0 ) THEN
            ix1 = ix1 + 1
         ELSEIF ( ix0==N ) THEN
            ix0 = ix0 - 1
         ELSE
            WRITE (G_IO,99020)
            WRITE (G_IO,99007)
99007       FORMAT (' ','LOWER AND UPPER BOUND IDENTICAL')
            GOTO 200
         ENDIF
      ENDIF
!
!     COMPUTE NEGATIVE BINOMIAL PROBABILITIES FOR THE
!     DERIVED LOWER AND UPPER BOUNDS.
!
      x0 = ix0
      x1 = ix1
      CALL NBCDF(x0,Ppar,N,p0)
      CALL NBCDF(x1,Ppar,N,p1)
!
!     CHECK THE PROBABILITIES FOR PROPER ORDERING
!
      IF ( p0<P .AND. P<=p1 ) THEN
         DO
!
!     THE STOPPING CRITERION IS THAT THE LOWER BOUND
!     AND UPPER BOUND ARE EXACTLY 1 UNIT APART.
!     CHECK TO SEE IF IX1 = IX0 + 1;
!     IF SO, THE ITERATIONS ARE COMPLETE;
!     IF NOT, THEN BISECT, COMPUTE PROBABILIIES,
!     CHECK PROBABILITIES, AND CONTINUE ITERATING
!     UNTIL IX1 = IX0 + 1.
!
            ix0p1 = ix0 + 1
            IF ( ix1==ix0p1 ) THEN
               Ppf = ix1
               IF ( p0==P ) Ppf = ix0
               RETURN
            ELSE
               ix2 = (ix0+ix1)/2
               IF ( ix2/=ix0 ) THEN
                  IF ( ix2==ix1 ) THEN
                     WRITE (G_IO,99020)
                     WRITE (G_IO,99021)
                     EXIT
                  ELSE
                     x2 = ix2
                     CALL NBCDF(x2,Ppar,N,p2)
                     IF ( p0<p2 .AND. p2<p1 ) THEN
                        IF ( p2<=P ) THEN
                           ix0 = ix2
                           p0 = p2
                        ELSE
                           ix1 = ix2
                           p1 = p2
                        ENDIF
                        CYCLE
                     ELSEIF ( p2<=p0 ) THEN
                        WRITE (G_IO,99020)
                        WRITE (G_IO,99008)
99008                   FORMAT (' ','BISECTION VALUE PROBABILITY (P2) ',&
     &                          'LESS THAN LOWER BOUND PROBABILITY (P0)'&
     &                          )
                        EXIT
                     ELSEIF ( p2>=p1 ) THEN
                        WRITE (G_IO,99020)
                        WRITE (G_IO,99009)
99009                   FORMAT (' ','BISECTION VALUE PROBABILITY (P2) ',&
     &                       'GREATER THAN UPPER BOUND PROBABILITY (P1)'&
     &                       )
                        EXIT
                     ENDIF
                  ENDIF
               ENDIF
               WRITE (G_IO,99020)
               WRITE (G_IO,99021)
               EXIT
            ENDIF
         ENDDO
      ELSEIF ( p0==P ) THEN
         Ppf = ix0
         RETURN
      ELSEIF ( p1==P ) THEN
         Ppf = ix1
         RETURN
      ELSEIF ( p0>p1 ) THEN
         WRITE (G_IO,99020)
         WRITE (G_IO,99010)
99010    FORMAT (' ','LOWER BOUND PROBABILITY (P0) GREATER THAN ',      &
     &           'UPPER BOUND PROBABILITY (P1)')
      ELSEIF ( p0>P ) THEN
         WRITE (G_IO,99020)
         WRITE (G_IO,99011)
99011    FORMAT (' ','LOWER BOUND PROBABILITY (P0) GREATER THAN ',      &
     &           'INPUT PROBABILITY (P)')
      ELSEIF ( p1<P ) THEN
         WRITE (G_IO,99020)
         WRITE (G_IO,99012)
99012    FORMAT (' ','UPPER BOUND PROBABILITY (P1) LESS    THAN ',      &
     &           'INPUT PROBABILITY (P)')
      ELSE
         WRITE (G_IO,99020)
         WRITE (G_IO,99013)
99013    FORMAT (' ','IMPOSSIBLE BRANCH CONDITION ENCOUNTERED')
      ENDIF
!
 200  WRITE (G_IO,99014) ix0 , p0
99014 FORMAT (' ','IX0  = ',I8,10X,'P0 = ',F14.7)
      WRITE (G_IO,99015) ix1 , p1
99015 FORMAT (' ','IX1  = ',I8,10X,'P1 = ',F14.7)
      WRITE (G_IO,99016) ix2 , p2
99016 FORMAT (' ','IX2  = ',I8,10X,'P2 = ',F14.7)
      WRITE (G_IO,99017) P
99017 FORMAT (' ','P    = ',F14.7)
      WRITE (G_IO,99018) Ppar , N
99018 FORMAT (' ','PPAR = ',F14.7,10X,'N  = ',I0)
      RETURN
99019 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,' *****')
99020 FORMAT (' ','***** INTERNAL ERROR IN NBPPF  SUBROUTINE *****')
99021 FORMAT (' ','BISECTION VALUE (X2) = LOWER BOUND (X0)')
99022 FORMAT (' ','BISECTION VALUE (X2) = UPPER BOUND (X1)')
!
END SUBROUTINE NBPPF