POIPPF – M

public subroutine POIPPF(P, Alamba, Ppf)

NAME

poippf(3f) - [M_datapac:PERCENT_POINT] compute the Poisson percent
point function

SYNOPSIS

   SUBROUTINE POIPPF(P,Alamba,Ppf)

DESCRIPTION

POIPPF(3f) computes the percent point function value at the precision
precision value P for the Poisson distribution with REAL
tail length parameter = alamba.

the poisson distribution used herein has mean = alamba and standard
deviation = sqrt(alamba). this distribution is defined for all
discrete non-negative integer x--x = 0, 1, 2, ... .

this distribution has the probability function

    f(x) = exp(-alamba) * alamba**x / x!.

the poisson distribution is the distribution of the number of events
in the interval (0,alamba) when the waiting time between events is
exponentially distributed with mean = 1 and standard deviation = 1.

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_poippf
use M_datapac, only : poippf
implicit none
! call poippf(x,y)
end program demo_poippf

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 87-121, especially page 102, Formula 36.1. –Hastings and Peacock, Statistical Distributions–A Handbook for Students and Practitioners, 1975, pages 108-113.
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 146-154.
Cox and Miller, The Theory of Stochastic Processes, 1965, page 7.
General Electric Company, Tables of the Individual and Cumulative Terms of Poisson Distribution, 1962.
Owen, Handbook of Statistical Tables, 1962, pages 259-261.

Arguments

Type		Name
real(kind=wp)	::	P
real(kind=wp)	::	Alamba
real(kind=wp)	::	Ppf

Source Code

      SUBROUTINE POIPPF(P,Alamba,Ppf)
REAL(kind=wp) :: Alamba, amean, P, p0, p1, p2, pf0, Ppf, sd, x0, x1, x2, zppf
INTEGER :: i, isd, ix0, ix0p1, ix1, ix2
!
!     INPUT ARGUMENTS--P      = THE  VALUE
!                                (BETWEEN 0.0 (INCLUSIVELY)
!                                AND 1.0 (EXCLUSIVELY))
!                                AT WHICH THE PERCENT POINT
!                                FUNCTION IS TO BE EVALUATED.
!                     --ALAMBA = THE  VALUE
!                                OF THE TAIL LENGTH PARAMETER.
!                                ALAMBA SHOULD BE POSITIVE.
!     OUTPUT ARGUMENTS--PPF    = THE  PERCENT
!                                POINT FUNCTION VALUE.
!     OUTPUT--THE  PERCENT POINT  .
!             FUNCTION VALUE PPF
!             FOR THE POISSON DISTRIBUTION
!             WITH TAIL LENGTH PARAMETER = ALAMBA.
!     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
!     RESTRICTIONS--ALAMBA SHOULD BE POSITIVE.
!                 --P SHOULD BE BETWEEN 0.0 (INCLUSIVELY)
!                   AND 1.0 (EXCLUSIVELY).
!     OTHER DATAPAC   SUBROUTINES NEEDED--NORPPF, POICDF.
!     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, DEXP.
!     MODE OF INTERNAL OPERATIONS-- AND DOUBLE PRECISION.
!     COMMENT--THE  TAIL LENGTH
!              PARAMETER ALAMBA IS     NOT     RESTRICTED
!              TO ONLY INTEGER VALUES.
!              ALAMBA CAN BE SET TO ANY POSITIVE REAL
!              VALUE--INTEGER OR NON-INTEGER.
!            --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 dlamba
!
!     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 POIPPF(3f) IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****')
         WRITE (G_IO,99017) P
         Ppf = 0.0_wp
         RETURN
      ELSE
         IF ( Alamba<=0.0_wp ) THEN
            WRITE (G_IO,99002)
99002       FORMAT (' ***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO POIPPF(3f) IS NON-POSITIVE *****')
            WRITE (G_IO,99017) Alamba
            Ppf = 0.0_wp
            RETURN
         ELSE
!
!-----START POINT-----------------------------------------------------
!
            dlamba = Alamba
            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_wp
!     2) PPF = 0
!
            IF ( P/=0.0_wp ) THEN
               pf0 = DEXP(-dlamba)
               IF ( P>pf0 ) THEN
!
!     DETERMINE AN INITIAL APPROXIMATION TO THE POISSON
!     PERCENT POINT BY USE OF THE NORMAL APPROXIMATION
!     TO THE POISSON.
!     (SEE JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS,
!     page 102, FORMULA 36.1).
!
                  amean = Alamba
                  sd = SQRT(Alamba)
                  CALL NORPPF(P,zppf)
                  x2 = amean - 1.0_wp + zppf*sd
                  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 POICDF(x2,Alamba,p2)

                  IF ( p2<P ) THEN

                     ix0 = ix2
                     DO i = 1, 100000
                        ix2 = ix0 + isd
                        IF ( ix2>=ix1 ) GOTO 200
                        x2 = ix2
                        CALL POICDF(x2,Alamba,p2)
                        IF ( p2>=P ) GOTO 50
                        ix0 = ix2
                     ENDDO
                     WRITE (G_IO,99018)
                     WRITE (G_IO,99003)
                     99003 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 POICDF(x2,Alamba,p2)
                        IF ( p2<P ) GOTO 100
                        ix1 = ix2
                     ENDDO
                     WRITE (G_IO,99018)
                     WRITE (G_IO,99004)
                     99004 FORMAT (' NO LOWER BOUND FOUND AFTER 10**7 ITERATIONS')
                  ENDIF
                  GOTO 400
               ENDIF
            ENDIF
            Ppf = 0.0_wp
            RETURN
         ENDIF
 50      ix1 = ix2
         GOTO 200
      ENDIF
 100  ix0 = ix2
!
 200  IF ( ix0==ix1 ) THEN
         IF ( ix0==0 ) THEN
            ix1 = ix1 + 1
            GOTO 300
            ix0 = ix0 - 1
         ELSE
            WRITE (G_IO,99018)
            WRITE (G_IO,99005)
            99005 FORMAT (' ','LOWER AND UPPER BOUND IDENTICAL')
            GOTO 400
         ENDIF
      ENDIF
!
!     COMPUTE POISSON PROBABILITIES FOR THE
!     DERIVED LOWER AND UPPER BOUNDS.
!
 300  x0 = ix0
      x1 = ix1
      CALL POICDF(x0,Alamba,p0)
      CALL POICDF(x1,Alamba,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,99018)
                     WRITE (G_IO,99019)
                     EXIT
                  ELSE
                     x2 = ix2
                     CALL POICDF(x2,Alamba,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,99018)
                        WRITE (G_IO,99006)
                        99006 FORMAT (' BISECTION VALUE PROBABILITY (P2) LESS THAN LOWER BOUND PROBABILITY (P0)')
                        EXIT
                     ELSEIF ( p2>=p1 ) THEN
                        WRITE (G_IO,99018)
                        WRITE (G_IO,99007)
                        99007 FORMAT (' BISECTION VALUE PROBABILITY (P2) GREATER THAN UPPER BOUND PROBABILITY (P1)')
                        EXIT
                     ENDIF
                  ENDIF
               ENDIF
               WRITE (G_IO,99018)
               WRITE (G_IO,99019)
               EXIT
            ENDIF
         ENDDO
      ELSEIF ( p0==P ) THEN
         Ppf = ix0
         RETURN
      ELSEIF ( p1==P ) THEN
         Ppf = ix1
         RETURN
      ELSEIF ( p0>p1 ) THEN
         WRITE (G_IO,99018)
         WRITE (G_IO,99008)
         99008 FORMAT (' ','LOWER BOUND PROBABILITY (P0) GREATER THAN UPPER BOUND PROBABILITY (P1)')
      ELSEIF ( p0>P ) THEN
         WRITE (G_IO,99018)
         WRITE (G_IO,99009)
         99009 FORMAT (' ','LOWER BOUND PROBABILITY (P0) GREATER THAN INPUT PROBABILITY (P)')
      ELSEIF ( p1<P ) THEN
         WRITE (G_IO,99018)
         WRITE (G_IO,99010)
         99010 FORMAT (' ','UPPER BOUND PROBABILITY (P1) LESS THAN INPUT PROBABILITY (P)')
      ELSE
         WRITE (G_IO,99018)
         WRITE (G_IO,99011)
         99011 FORMAT (' ','IMPOSSIBLE BRANCH CONDITION ENCOUNTERED')
      ENDIF
!
 400  WRITE (G_IO,99012) ix0, p0
      99012 FORMAT (' ','IX0    = ',I8,10X,'P0 = ',F14.7)
      WRITE (G_IO,99013) ix1, p1
      99013 FORMAT (' ','IX1    = ',I8,10X,'P1 = ',F14.7)
      WRITE (G_IO,99014) ix2, p2
      99014 FORMAT (' ','IX2    = ',I8,10X,'P2 = ',F14.7)
      WRITE (G_IO,99015) P
      99015 FORMAT (' ','P      = ',F14.7)
      WRITE (G_IO,99016) Alamba
      99016 FORMAT (' ','ALAMBA = ',F14.7)
      RETURN
99017 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,' *****')
99018 FORMAT (' ','***** INTERNAL ERROR IN POIPPF SUBROUTINE *****')
99019 FORMAT (' ','BISECTION VALUE (X2) = LOWER BOUND (X0)')
99020 FORMAT (' ','BISECTION VALUE (X2) = UPPER BOUND (X1)')
!
END SUBROUTINE POIPPF

POIPPF Subroutine

Contents

Source Code

public subroutine POIPPF(P, Alamba, Ppf)

NAME

SYNOPSIS

DESCRIPTION

OPTIONS

EXAMPLES

AUTHOR

MAINTAINER

LICENSE

REFERENCES

Arguments

Source Code