BINPPF – M

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

NAME

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

SYNOPSIS

   SUBROUTINE BINPPF(P,Ppar,N,Ppf)

    REAL(kind=wp) :: P
    REAL(kind=wp) :: Ppar
    REAL(kind=wp) :: Ppf
    INTEGER :: N

DESCRIPTION

BINPPF(3f) computes the percent point function value at the precision
precision value P for the binomial distribution with REAL
'Bernoulli probability' parameter = PPAR, and integer 'number of
Bernoulli trials' parameter = N.

The binomial distribution used herein has mean = N*PPAR and standard
deviation = sqrt(N*PPAR*(1-PPAR)).

This distribution is defined for all discrete integer X between 0
(inclusively) and N (inclusively).

This distribution has the probability function

   f(X) = c(N,X) * PPAR**X * (1-PPAR)**(N-X).

where c(N,X) is the combinatorial function equaling the number of
combinations of N items taken X at a time.

The binomial distribution is the distribution of the number of
successes in N 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.

INPUT ARGUMENTS

P     The value (between 0.0 (inclusively) and 1.0 (inclusively))
      at which the percent point function is to be evaluated.
PPAR  The value of the 'Bernoulli probability' parameter for the binomial
      distribution. PPAR should be between 0.0 (exclusively) and
      1.0 (exclusively).
N     The integer value of the 'number of Bernoulli trials' parameter.
      N should be a positive integer.

OUTPUT ARGUMENTS

PPF The percent point function value.

EXAMPLES

Sample program:

program demo_binppf
use M_datapac, only : binppf
implicit none
! call binppf(x,y)
end program demo_binppf

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 50-86, especially page 64, Formula 36.
Hastings and Peacock, Statistical Distributions–A Handbook for Students and Practitioners, 1975, pages 36-41.
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 135-142.
Kendall and Stuart, The Advanced Theory of Statistics, Volume 1, Edition 2, 1963, pages 120-125.
Mood and Grable, Introduction to the Theory of Statistics, Edition 2, 1963, pages 64-69.
Owen, Handbook of Statistical Tables, 1962, pages 264-272.

Arguments

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

Source Code

SUBROUTINE BINPPF(P,Ppar,N,Ppf)
REAL(kind=wp) :: P
REAL(kind=wp) :: Ppar
REAL(kind=wp) :: Ppf
INTEGER :: N
REAL(kind=wp) :: amean , an, p0 , p1 , p2 , pf0 , qfn , sd , x0 , x1 , x2 , zppf
INTEGER :: i , isd , ix0 , ix0p1 , ix1 , ix2
DOUBLE PRECISION :: dppar


!     MODE OF INTERNAL OPERATIONS -- SINGLE 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.
!---------------------------------------------------------------------
!
!     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 BINPPF(3f) 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 BINPPF(3f) 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 BINPPF(3f) 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 OR 1.0
!     2) P = 0.5 AND PPAR = 0.5
!     3) PPF = 0 OR N
!
            IF ( P/=0.0_wp ) THEN
               IF ( P==1.0_wp ) GOTO 20
               IF ( P==0.5_wp .AND. Ppar==0.5_wp ) THEN
                  Ppf = N/2
                  RETURN
               ELSE
                  pf0 = (1.0D0-dppar)**N
                  qfn = 1.0D0 - (dppar**N)
                  IF ( P>pf0 ) THEN
                     IF ( P>qfn ) GOTO 20
!
!     DETERMINE AN INITIAL APPROXIMATION TO THE BINOMIAL
!     PERCENT POINT BY USE OF THE NORMAL APPROXIMATION
!     TO THE BINOMIAL.
!     (SEE JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS,
!     page 64, FORMULA 36).
!
                     amean = an*Ppar
                     sd = SQRT(an*Ppar*(1.0_wp-Ppar))
                     CALL NORPPF(P,zppf)
                     x2 = amean - 0.5_wp + zppf*sd
                     ix2 = x2
!
!     CHECK AND MODIFY (IF NECESSARY) THIS INITIAL
!     ESTIMATE OF THE PERCENT POINT
!     TO ASSURE THAT IT BE IN THE CLOSED INTERVAL 0 TO N.
!
                     IF ( ix2<0 ) ix2 = 0
                     IF ( ix2>N ) ix2 = N
!
!     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 = N
                     isd = sd + 1.0_wp
                     x2 = ix2
                     CALL BINCDF(x2,Ppar,N,p2)
!
                     IF ( p2<P ) THEN
!
                        ix0 = ix2
                        DO i = 1 , 100000
                           ix2 = ix0 + isd
                           IF ( ix2>=ix1 ) GOTO 200
                           x2 = ix2
                           CALL BINCDF(x2,Ppar,N,p2)
                           IF ( p2>=P ) GOTO 50
                           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 200
                           x2 = ix2
                           CALL BINCDF(x2,Ppar,N,p2)
                           IF ( p2<P ) GOTO 100
                           ix1 = ix2
                        ENDDO
                        WRITE (G_IO,99020)
                        WRITE (G_IO,99006)
                        99006 FORMAT (' NO LOWER BOUND FOUND AFTER 10**7 ITERATIONS')
                     ENDIF
                     GOTO 300
                  ENDIF
               ENDIF
            ENDIF
            Ppf = 0.0_wp
            RETURN
         ENDIF
 20      Ppf = N
         RETURN
      ENDIF
 50   ix1 = ix2
      GOTO 200
   ENDIF
 100  ix0 = ix2
!
 200  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 300
      ENDIF
   ENDIF
!
!     COMPUTE BINOMIAL PROBABILITIES FOR THE
!     DERIVED LOWER AND UPPER BOUNDS.
!
   x0 = ix0
   x1 = ix1
   CALL BINCDF(x0,Ppar,N,p0)
   CALL BINCDF(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 BINCDF(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
!
 300  continue
   WRITE (G_IO,99014) ix0 , p0
   99014 FORMAT (' ','IX0  = ',I0,10X,'P0 = ',F14.7)
   WRITE (G_IO,99015) ix1 , p1
   99015 FORMAT (' ','IX1  = ',I0,10X,'P1 = ',F14.7)
   WRITE (G_IO,99016) ix2 , p2
   99016 FORMAT (' ','IX2  = ',I0,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 BINPPF SUBROUTINE *****')
   99021 FORMAT (' ','BISECTION VALUE (X2) = LOWER BOUND (X0)')
   99022 FORMAT (' ','BISECTION VALUE (X2) = UPPER BOUND (X1)')
!
END SUBROUTINE BINPPF

BINPPF Subroutine

Contents

Source Code

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

NAME

SYNOPSIS

DESCRIPTION

INPUT ARGUMENTS

OUTPUT ARGUMENTS

EXAMPLES

AUTHOR

MAINTAINER

LICENSE

REFERENCES

Arguments

Source Code