PROGRAM test_kemp_binomial IMPLICIT NONE INTEGER :: n, nrep, i, ix REAL :: p, mean = 0.0, sumsq = 0.0, stdev = 0.0, dev LOGICAL :: first INTERFACE FUNCTION random_binomial(n, p, first) RESULT(ran_binom) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: p LOGICAL, INTENT(IN) :: first INTEGER :: ran_binom END FUNCTION random_binomial END INTERFACE WRITE(*, *) 'Enter n & p: ' READ(*, *) n, p WRITE(*, *) 'How many random deviates?: ' READ(*, *) nrep first = .TRUE. DO i = 1, nrep ix = random_binomial(n, p, first) dev = ix - mean mean = mean + dev/i sumsq = sumsq + dev*(ix - mean) first = .FALSE. END DO IF (nrep > 1) stdev = SQRT(sumsq / (nrep-1)) WRITE(*, '(1x, "Mean = ", f8.2, 5x, "Std. devn. = ", f8.2)') mean, stdev STOP END PROGRAM test_kemp_binomial FUNCTION random_binomial(n, p, first) RESULT(ran_binom) ! FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method. ! This algorithm is suitable when many random variates are required ! with the SAME parameter values for n & p. ! P = BERNOULLI SUCCESS PROBABILITY ! (0 <= REAL <= 1) ! N = NUMBER OF BERNOULLI TRIALS ! (1 <= INTEGER) ! FIRST = .TRUE. for the first call using the current parameter values ! = .FALSE. if the values of (n,p) are unchanged from last call ! Reference: Kemp, C.D. (1986). `A modal method for generating binomial ! variables', Commun. Statist. - Theor. Meth. 15(3), 805-813. IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: p LOGICAL, INTENT(IN) :: first INTEGER :: ran_binom INTERFACE FUNCTION bin_prob(n, p, r) RESULT(prob) IMPLICIT NONE INTEGER, INTENT(IN) :: n, r REAL, INTENT(IN) :: p REAL :: prob END FUNCTION bin_prob END INTERFACE ! Local variables INTEGER :: ru, rd INTEGER, SAVE :: r0 REAL :: u, pd, pu, zero = 0.0, one = 1.0 REAL, SAVE :: odds_ratio, p_r IF (first) THEN r0 = (n+1)*p p_r = bin_prob(n, p, r0) odds_ratio = p / (one - p) END IF CALL RANDOM_NUMBER(u) u = u - p_r IF (u < zero) THEN ran_binom = r0 RETURN END IF pu = p_r ru = r0 pd = p_r rd = r0 DO rd = rd - 1 IF (rd >= 0) THEN pd = pd * (rd+1) / (odds_ratio * REAL(n-rd)) u = u - pd IF (u < zero) THEN ran_binom = rd RETURN END IF END IF ru = ru + 1 IF (ru <= n) THEN pu = pu * (n-ru+1) * odds_ratio / REAL(ru) u = u - pu IF (u < zero) THEN ran_binom = ru RETURN END IF END IF END DO ! This point should not be reached, but just in case: ran_binom = r0 RETURN END FUNCTION random_binomial FUNCTION bin_prob(n, p, r) RESULT(prob) ! Calculate a binomial probability IMPLICIT NONE INTEGER, INTENT(IN) :: n, r REAL, INTENT(IN) :: p REAL :: prob ! Local variable REAL :: one = 1.0 prob = EXP( lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1)) & + r*LOG(p) + (n-r)*LOG(one - p) ) RETURN CONTAINS FUNCTION lngamma(z) RESULT(lanczos) ! Uses Lanczos-type approximation to ln(gamma) for z > 0. ! Reference: ! Lanczos, C. 'A precision approximation of the gamma ! function', J. SIAM Numer. Anal., B, 1, 86-96, 1964. ! Accuracy: About 14 significant digits except for small regions ! in the vicinity of 1 and 2. ! Programmer: Alan Miller ! 1 Creswick Street, Brighton, Vic. 3187, Australia ! Latest revision - 14 October 1996 IMPLICIT NONE INTEGER, PARAMETER :: doub_prec = SELECTED_REAL_KIND(15, 60) REAL(doub_prec), INTENT(IN) :: z REAL(doub_prec) :: lanczos ! Local variables REAL(doub_prec) :: a(9) = (/ 0.9999999999995183D0, 676.5203681218835D0, & -1259.139216722289D0, 771.3234287757674D0, & -176.6150291498386D0, 12.50734324009056D0, & -0.1385710331296526D0, 0.9934937113930748D-05, & 0.1659470187408462D-06 /), zero = 0.D0, & one = 1.d0, lnsqrt2pi = 0.9189385332046727D0, & half = 0.5d0, sixpt5 = 6.5d0, seven = 7.d0, tmp INTEGER :: j IF (z <= zero) THEN WRITE(*, *)'Error: zero or -ve argument for lngamma' RETURN END IF lanczos = zero tmp = z + seven DO j = 9, 2, -1 lanczos = lanczos + a(j)/tmp tmp = tmp - one END DO lanczos = lanczos + a(1) lanczos = LOG(lanczos) + lnsqrt2pi - (z + sixpt5) + (z - half)*LOG(z + sixpt5) RETURN END FUNCTION lngamma END FUNCTION bin_prob