bincdf(3f) - [M_datapac:CUMULATIVE_DISTRIBUTION] compute the binomial
cumulative distribution function
SUBROUTINE BINCDF(X,P,N,Cdf)
REAL(kind=wp) :: X
REAL(kind=wp) :: P
INTEGER :: N
REAL(kind=wp) :: Cdf
BINCDF(3f) computes the cumulative distribution function value at the
double precision value X for the binomial distribution with double
precision 'Bernoulli probability' parameter = P, and integer 'number
of Bernoulli trials' parameter = N.
The binomial distribution used herein has mean = N*P and standard
deviation = sqrt(N*P*(1-P)).
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) * P**X * (1-P)**(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 = P.
X The value at which the cumulative distribution
function is to be evaluated. X should be integral-valued,
and between 0.0 (inclusively) and N (inclusively).
P The value of the 'Bernoulli probability'
parameter for the binomial distribution.
P 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.
CDF The cumulative distribution function value for the binomial distribution.
Sample program:
program demo_bincdf
use M_datapac, only : bincdf
implicit none
!call BINCDF(X,P,N,Cdf)
end program demo_bincdf
Results:
The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.
John Urban, 2022.05.31
CC0-1.0
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp) | :: | X | ||||
| real(kind=wp) | :: | P | ||||
| integer | :: | N | ||||
| real(kind=wp) | :: | Cdf |
SUBROUTINE BINCDF(X,P,N,Cdf) REAL(kind=wp) :: X REAL(kind=wp) :: P INTEGER :: N REAL(kind=wp) :: Cdf REAL(kind=wp) :: an, del, fintx INTEGER :: i, ievodd, iflag1, iflag2, imax, imin, intx, nu1, nu2 DOUBLE PRECISION :: dx , anu1 , anu2 , z , sum , term , ai , coef1 , coef2 , arg DOUBLE PRECISION :: coef DOUBLE PRECISION :: theta , sinth , costh , a , b DOUBLE PRECISION :: DSQRT , DATAN ! COMMENT--NOTE THAT EVEN THOUGH THE INPUT ! TO THIS CUMULATIVE ! DISTRIBUTION FUNCTION SUBROUTINE ! FOR THIS DISCRETE DISTRIBUTION ! SHOULD (UNDER NORMAL CIRCUMSTANCES) BE A ! DISCRETE INTEGER VALUE, ! THE INPUT VARIABLE X IS SINGLE ! PRECISION IN MODE. ! X HAS BEEN SPECIFIED AS SINGLE ! PRECISION SO AS TO CONFORM WITH THE DATAPAC ! CONVENTION THAT ALL INPUT ****DATA**** ! (AS OPPOSED TO SAMPLE SIZE, FOR EXAMPLE) ! VARIABLES TO 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 ! an = N IF ( P<0.0_wp .OR. P>1.0_wp ) THEN WRITE (G_IO,99001) 99001 FORMAT (' ', & &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE BINCDF SUBROU& &TINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****') WRITE (G_IO,99006) P Cdf = 0.0_wp RETURN ELSEIF ( N<1 ) THEN WRITE (G_IO,99002) 99002 FORMAT (' ', & &'***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE BINCDF SUBROU& &TINE IS NON-POSITIVE *****') WRITE (G_IO,99003) N 99003 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****') Cdf = 0.0_wp RETURN ELSEIF ( X<0.0_wp .OR. X>an ) THEN WRITE (G_IO,99004) N 99004 FORMAT (' ', & &'***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUMENT TO THE BINC& &DF SUBROUTINE IS OUTSIDE THE USUAL (0,N) = (0,',I0,',INTERVAL *') WRITE (G_IO,99006) X IF ( X<0.0_wp ) Cdf = 0.0_wp IF ( X>an ) Cdf = 1.0_wp RETURN ELSE intx = X + 0.0001_wp fintx = intx del = X - fintx IF ( del<0.0_wp ) del = -del IF ( del>0.001_wp ) THEN WRITE (G_IO,99005) 99005 FORMAT (' ', & &'***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUMENT TO THE BINC& &DF SUBROUTINE IS NON-INTEGRAL *****') WRITE (G_IO,99006) X ENDIF ! !-----START POINT----------------------------------------------------- ! ! TREAT IMMEDIATELY THE SPECIAL CASE OF X = N, ! IN WHICH CASE CDF = 1.0. ! ALSO TREAT IMMEDIATELY THE SPECIAL CASE OF P = 0.0_wp ! IN WHICH CASE CDF = 1.0 FOR ALL X. ! THIRDLY, TREAT THE SPECIAL CASE IN WHICH P = 1.0 ! IN WHICH CASE CDF = 0.0 FOR ALL X SMALLER THAN N ! AND CDF = 1.0 FOR ALL X EQUAL TO OR LARGER ! THAN N. ! intx = X + 0.0001_wp Cdf = 1.0_wp IF ( intx==N ) RETURN IF ( P==0.0_wp ) RETURN IF ( P==1.0_wp .AND. intx>=N ) RETURN IF ( P==1.0_wp .AND. intx<N ) Cdf = 0.0_wp IF ( P==1.0_wp .AND. intx<N ) RETURN ! ! EXPRESS THE BINOMIAL CUMULATIVE DISTRIBUTION ! FUNCTION IN TERMS OF THE EQUIVALENT F ! CUMULATIVE DISTRIBUTION FUNCTION, ! AND THEN EVALUATE THE LATTER. ! an = N dx = (P/(1.0_wp-P))*((an-X)/(X+1.0_wp)) nu1 = 2.0_wp*(X+1.0_wp) + 0.1_wp nu2 = 2.0_wp*(an-X) + 0.1_wp anu1 = nu1 anu2 = nu2 z = anu2/(anu2+anu1*dx) ! ! DETERMINE IF NU1 AND NU2 ARE EVEN OR ODD ! iflag1 = nu1 - 2*(nu1/2) iflag2 = nu2 - 2*(nu2/2) IF ( iflag1==0 ) THEN ! ! DO THE NU1 EVEN AND NU2 EVEN OR ODD CASE ! sum = 0.0D0 term = 1.0D0 imax = (nu1-2)/2 IF ( imax>0 ) THEN DO i = 1 , imax ai = i coef1 = 2.0D0*(ai-1.0D0) coef2 = 2.0D0*ai term = term*((anu2+coef1)/coef2)*(1.0D0-z) sum = sum + term ENDDO ENDIF ! sum = sum + 1.0D0 sum = (z**(anu2/2.0D0))*sum Cdf = sum RETURN ELSEIF ( iflag2==0 ) THEN ! ! DO THE NU1 ODD AND NU2 EVEN CASE ! sum = 0.0D0 term = 1.0D0 imax = (nu2-2)/2 IF ( imax>0 ) THEN DO i = 1 , imax ai = i coef1 = 2.0D0*(ai-1.0D0) coef2 = 2.0D0*ai term = term*((anu1+coef1)/coef2)*z sum = sum + term ENDDO ENDIF ! sum = sum + 1.0D0 Cdf = 1.0D0 - ((1.0D0-z)**(anu1/2.0D0))*sum RETURN ELSE ! ! DO THE NU1 ODD AND NU2 ODD CASE ! sum = 0.0D0 term = 1.0D0 arg = DSQRT((anu1/anu2)*dx) theta = DATAN(arg) sinth = arg/DSQRT(1.0D0+arg*arg) costh = 1.0D0/DSQRT(1.0D0+arg*arg) IF ( nu2/=1 ) THEN IF ( nu2/=3 ) THEN imax = nu2 - 2 DO i = 3 , imax , 2 ai = i coef1 = ai - 1.0D0 coef2 = ai term = term*(coef1/coef2)*(costh*costh) sum = sum + term ENDDO ENDIF ! sum = sum + 1.0D0 sum = sum*sinth*costh ENDIF ! a = (2.0D0/G_pi_dp)*(theta+sum) sum = 0.0D0 term = 1.0D0 IF ( nu1==1 ) b = 0.0D0 IF ( nu1/=1 ) THEN IF ( nu1/=3 ) THEN imax = nu1 - 3 DO i = 1 , imax , 2 ai = i coef1 = ai coef2 = ai + 2.0D0 term = term*((anu2+coef1)/coef2)*(sinth*sinth) sum = sum + term ENDDO ENDIF ! sum = sum + 1.0D0 sum = sum*sinth*(costh**N) coef = 1.0D0 ievodd = nu2 - 2*(nu2/2) imin = 3 IF ( ievodd==0 ) imin = 2 IF ( imin<=nu2 ) THEN DO i = imin , nu2 , 2 ai = i coef = ((ai-1.0D0)/ai)*coef ENDDO ENDIF ! coef = coef*anu2 IF ( ievodd /= 0 ) coef = coef*(2.0D0/G_pi_dp) ! b = coef*sum ENDIF ! Cdf = 1.0D0 - (a-b) ENDIF ENDIF 99006 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,' *****') ! END SUBROUTINE BINCDF