random_binomial1 Function

FUNCTION random_binomial1(n, p, first) RESULT(ival)
  ! 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.
  INTEGER, INTENT(IN) :: n
  REAL, INTENT(IN)    :: p
  LOGICAL, INTENT(IN) :: first
  INTEGER             :: ival
  !     Local variables
  INTEGER         :: ru, rd
  INTEGER, SAVE   :: r0
  REAL            :: u, pd, pu
  REAL, SAVE      :: odds_ratio, p_r
  REAL, PARAMETER :: zero = 0.0, one = 1.0
  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
     ival = 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 * (n-rd))
        u = u - pd
        IF (u < zero) THEN
           ival = rd
           RETURN
        END IF
     END IF
     ru = ru + 1
     IF (ru <= n) THEN
        pu = pu * (n-ru+1) * odds_ratio / ru
        u = u - pu
        IF (u < zero) THEN
           ival = ru
           RETURN
        END IF
     END IF
  END DO
  !     This point should not be reached, but just in case:
  ival = r0
  RETURN
END FUNCTION random_binomial1