program bbml_prb
!
!*******************************************************************************
!
!! BBML_PRB controls the test of BBML.
!
  implicit none
!
  real a
  real a_est
  real b
  real b_est
  integer c
  character ( len = 8 ) date
  real mu
  real mu_est
  integer sample_log
  integer sample_num
  real theta
  real theta_est
  character ( len = 10 ) time
!
  a = 2.0E+00
  b = 3.0E+00
  c = 4

  call beta_binomial_check ( a, b, c )

  mu = a / ( a + b )
  theta = 1.0E+00 / ( a + b )

  call date_and_time ( date, time )

  write ( *, * ) ' '
  write ( *, * ) 'BBML_PRB'
  write ( *, * ) '  Demonstration of "BBML",'
  write ( *, * ) '  Applied Statistics Algorithm 189'
  write ( *, * ) '  Estimate the parameters A and B of a '
  write ( *, * ) '  beta binomial distribution.'
  write ( *, * ) ' '
  write ( *, * ) '  Today''s date: ', date
  write ( *, * ) '  Today''s time: ', time
  write ( *, * ) ' '
  write ( *, * ) 'Exact values:'
  write ( *, * ) ' '
  write ( *, * ) '          A             B                MU         THETA'
  write ( *, * ) ' '
  write ( *, '(6x,4g14.6)' ) a, b, mu, theta
  write ( *, * ) ' '
  write ( *, * ) 'Estimated values:'
  write ( *, * ) ' '
  write ( *, * ) 'Samples   A_est         B_est        MU_est     THETA_est'
  write ( *, * ) ' '

  do sample_log = 2, 5

    sample_num = 10**sample_log

    call analyze ( sample_num, a, b, c, mu_est, theta_est )
!
!  Convert the BBML "THETA" and "MU" parameters to "A" and "B".
!
    a_est = mu_est / theta_est
    b_est = ( 1.0E+00 - mu_est ) / theta_est

    write ( *, '(i6,4g14.6)' ) sample_num, a_est, b_est, mu_est, theta_est

    end do

  write ( *, * ) ' '
  write ( *, * ) 'BBML_PRB'
  write ( *, * ) '  Normal end of BBML tests.'

  stop
end
subroutine analyze ( sample_num, a, b, c, mu_est, theta_est )
!
!*******************************************************************************
!
!! ANALYZE generates data and analyzes it with BBML.
!
!
!  Discussion:
!
!    The routine to generate the samples uses parameter A, B and C,
!    while BBML prefers a related form MU, THETA.  The calling routine
!    has to figure out how these are related.
!
!  Modified:
!
!    09 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer SAMPLE_NUM, the number of samples to generate.
!
!    Input, real A, real B, integer C, the parameters to use in
!    the beta binomial distribution.
!
!    Output, real MU_EST, THETA_EST, the estimates of MU and THETA
!    produced by BBML.
!
  implicit none
!
  integer, parameter :: mrl = 4
!
  integer sample_num
!
  real a
  real b
  integer c
  real ccrit
  integer ifault
  integer in(sample_num)
  integer iseed
  integer iter
  real mu_est
  real mu_se
  integer rl(mrl,3)
  real rnl
  integer sample_i
  real theta_est
  real theta_se
  integer x(sample_num)
!
!  Generate the sample data using the exact parameters A, B.
!
  call get_seed ( iseed )

  do sample_i = 1, sample_num
    call beta_binomial_sample ( a, b, c, iseed, x(sample_i) )
  end do
!
!  Analyze the sample data, trying to estimate the parameters
!  in the form "MU", "THETA".  Note that BBML expects to be told
!  the value of C (although C could vary from one sample to the next).
!
  in(1:sample_num) = c

  iter = 10
  ccrit = 0.001E+00

  call bbml ( sample_num, x, in, rl, mrl, iter, ccrit, mu_est, &
    theta_est, mu_se, theta_se, rnl, ifault )

  return
end
function beta ( a, b )
!
!*******************************************************************************
!
!! BETA returns the value of the Beta function.
!
!
!  Formula:
!
!    BETA(A,B) = ( GAMMA ( A ) * GAMMA ( B ) ) / GAMMA ( A + B )
!              = Integral ( 0 <= T <= 1 ) T**(A-1) (1-T)**(B-1) dT.
!
!  Modified:
!
!    10 July 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real A, B, the parameters of the function.
!    0.0 < A,
!    0.0 < B.
!
!    Output, real BETA, the value of the function.
!
  implicit none
!
  real a
  real b
  real beta
  real gamma_log
!
!  Check.
!
  if ( a <= 0.0E+00 .or. b <= 0.0E+00 ) then
    write ( *, * ) ' '
    write ( *, * ) 'BETA - Fatal error!'
    write ( *, * ) '  Both A and B must be greater than 0.'
    stop
  end if

  beta = exp ( gamma_log ( a ) + gamma_log ( b ) - gamma_log ( a + b ) )

  return
end
subroutine beta_binomial_cdf_inv ( cdf, a, b, c, x )
!
!*******************************************************************************
!
!! BETA_BINOMIAL_CDF_INV inverts the Beta Binomial CDF.
!
!
!  Discussion:
!
!    A simple-minded discrete approach is used.
!
!  Modified:
!
!    07 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real CDF, the value of the CDF.
!
!    Input, real A, B, parameters of the PDF.
!    0.0 < A,
!    0.0 < B.
!
!    Input, integer C, a parameter of the PDF.
!    0 <= C.
!
!    Output, integer X, the smallest X whose cumulative density function
!    is greater than or equal to CDF.
!
  implicit none
!
  real a
  real b
  real beta
  integer c
  real cdf
  real cum
  real pdf
  integer x
  integer y
!
  if ( cdf <= 0.0E+00 ) then

    x = 0

  else

    cum = 0.0E+00

    do y = 0, c

      pdf = beta ( a + real ( y ), b + real ( c - y ) ) / ( real ( c + 1 ) &
        * beta ( real ( y + 1), real ( c - y + 1 ) ) * beta ( a, b ) )

      cum = cum + pdf

      if ( cdf <= cum ) then
        x = y
        return
      end if

    end do

    x = c

  end if

  return
end
subroutine beta_binomial_check ( a, b, c )
!
!*******************************************************************************
!
!! BETA_BINOMIAL_CHECK checks the parameters of the Beta Binomial PDF.
!
!
!  Modified:
!
!    07 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real A, B, parameters of the PDF.
!    0.0 < A,
!    0.0 < B.
!
!    Input, integer C, a parameter of the PDF.
!    0 <= C.
!
  implicit none
!
  real a
  real b
  integer c
!
  if ( a <= 0.0E+00 ) then
    write ( *, * ) ' '
    write ( *, * ) 'BETA_BINOMIAL_CHECK - Fatal error!'
    write ( *, * ) '  A <= 0.'
    stop
  end if

  if ( b <= 0.0E+00 ) then
    write ( *, * ) ' '
    write ( *, * ) 'BETA_BINOMIAL_CHECK - Fatal error!'
    write ( *, * ) '  B <= 0.'
    stop
  end if

  if ( c < 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'BETA_BINOMIAL_CHECK - Fatal error!'
    write ( *, * ) '  C < 0.'
    stop
  end if

  return
end
subroutine beta_binomial_sample ( a, b, c, iseed, x )
!
!*******************************************************************************
!
!! BETA_BINOMIAL_SAMPLE samples the Beta Binomial CDF.
!
!
!  Modified:
!
!    07 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real A, B, parameters of the PDF.
!    0.0 < A,
!    0.0 < B.
!
!    Input, integer C, a parameter of the PDF.
!    0 <= C.
!
!    Input/output, integer ISEED, a seed for the random number generator.
!
!    Output, integer X, a sample of the PDF.
!
  implicit none
!
  real a
  real b
  integer c
  real cdf
  integer iseed
  integer x
!
  call r_random ( cdf, 0.0E+00, 1.0E+00, iseed )

  call beta_binomial_cdf_inv ( cdf, a, b, c, x )

  return
end
function gamma_log ( x )
!
!*******************************************************************************
!
!! GAMMA_LOG calculates the natural logarithm of GAMMA ( X ) for positive X.
!
!
!  Discussion:
!
!    Computation is based on an algorithm outlined in references 1 and 2.
!    The program uses rational functions that theoretically approximate
!    LOG(GAMMA(X)) to at least 18 significant decimal digits.  The
!    approximation for X > 12 is from reference 3, while approximations
!    for X < 12.0 are similar to those in reference 1, but are unpublished.
!    The accuracy achieved depends on the arithmetic system, the compiler,
!    intrinsic functions, and proper selection of the machine-dependent
!    constants.
!
!  Modified:
!
!    16 June 1999
!
!  Authors:
!
!    W. J. Cody and L. Stoltz
!    Argonne National Laboratory
!
!  References:
!
!    # 1)
!    W. J. Cody and K. E. Hillstrom,
!    Chebyshev Approximations for the Natural Logarithm of the Gamma Function,
!    Math. Comp.
!    Volume 21, 1967, pages 198-203.
!
!    # 2)
!    K. E. Hillstrom,
!    ANL/AMD Program ANLC366S, DGAMMA/DLGAMA,
!    May 1969.
!
!    # 3)
!    Hart, Et. Al.,
!    Computer Approximations,
!    Wiley and sons, New York, 1968.
!
!  Parameters:
!
!    Input, real X, the argument of the Gamma function.  X must be positive.
!
!    Output, real GAMMA_LOG, the logarithm of the Gamma function of X.
!    If X <= 0.0, or if overflow would occur, the program returns the
!    value XINF, the largest representable floating point number.
!
!*******************************************************************************
!
!  Explanation of machine-dependent constants
!
!  BETA   - radix for the floating-point representation.
!
!  MAXEXP - the smallest positive power of BETA that overflows.
!
!  XBIG   - largest argument for which LN(GAMMA(X)) is representable
!           in the machine, i.e., the solution to the equation
!             LN(GAMMA(XBIG)) = BETA**MAXEXP.
!
!  XINF   - largest machine representable floating-point number;
!           approximately BETA**MAXEXP.
!
!  EPS    - The smallest positive floating-point number such that
!           1.0+EPS .GT. 1.0
!
!  FRTBIG - Rough estimate of the fourth root of XBIG
!
!
!  Approximate values for some important machines are:
!
!                            BETA      MAXEXP         XBIG
!
!  CRAY-1        (S.P.)        2        8191       9.62E+2461
!  Cyber 180/855
!    under NOS   (S.P.)        2        1070       1.72E+319
!  IEEE (IBM/XT,
!    SUN, etc.)  (S.P.)        2         128       4.08E+36
!  IEEE (IBM/XT,
!    SUN, etc.)  (D.P.)        2        1024       2.55D+305
!  IBM 3033      (D.P.)       16          63       4.29D+73
!  VAX D-Format  (D.P.)        2         127       2.05D+36
!  VAX G-Format  (D.P.)        2        1023       1.28D+305
!
!
!                            XINF        EPS        FRTBIG
!
!  CRAY-1        (S.P.)   5.45E+2465   7.11E-15    3.13E+615
!  Cyber 180/855
!    under NOS   (S.P.)   1.26E+322    3.55E-15    6.44E+79
!  IEEE (IBM/XT,
!    SUN, etc.)  (S.P.)   3.40E+38     1.19E-7     1.42E+9
!  IEEE (IBM/XT,
!    SUN, etc.)  (D.P.)   1.79D+308    2.22D-16    2.25D+76
!  IBM 3033      (D.P.)   7.23D+75     2.22D-16    2.56D+18
!  VAX D-Format  (D.P.)   1.70D+38     1.39D-17    1.20D+9
!  VAX G-Format  (D.P.)   8.98D+307    1.11D-16    1.89D+76
!
  implicit none
!
  real, parameter :: d1 = - 5.772156649015328605195174e-1
  real, parameter :: d2 =   4.227843350984671393993777e-1
  real, parameter :: d4 =   1.791759469228055000094023e0
  real, parameter :: EPS = 1.19e-7
  real, parameter :: FRTBIG = 1.42e9
  real, parameter :: PNT68 = 0.6796875
  real, parameter :: SQRTPI = 0.9189385332046727417803297
  real, parameter :: XBIG = 4.08e36
  real, parameter :: XINF = 3.401e38
!
  real, parameter, dimension ( 7 ) :: c = (/ &
    -1.910444077728e-03, &
     8.4171387781295e-04, &
    -5.952379913043012e-04, &
     7.93650793500350248e-04, &
    -2.777777777777681622553e-03, &
     8.333333333333333331554247e-02, &
     5.7083835261e-03 /)
  real corr
  integer i
  real gamma_log
  real, parameter, dimension ( 8 ) :: p1 = (/ &
    4.945235359296727046734888e0, &
    2.018112620856775083915565e2, &
    2.290838373831346393026739e3, &
    1.131967205903380828685045e4, &
    2.855724635671635335736389e4, &
    3.848496228443793359990269e4, &
    2.637748787624195437963534e4, &
    7.225813979700288197698961e3 /)
  real, parameter, dimension ( 8 ) :: p2 = (/ &
    4.974607845568932035012064e0, &
    5.424138599891070494101986e2, &
    1.550693864978364947665077e4, &
    1.847932904445632425417223e5, &
    1.088204769468828767498470e6, &
    3.338152967987029735917223e6, &
    5.106661678927352456275255e6, &
    3.074109054850539556250927e6 /)
  real, parameter, dimension ( 8 ) :: p4 = (/ &
    1.474502166059939948905062e4, &
    2.426813369486704502836312e6, &
    1.214755574045093227939592e8, &
    2.663432449630976949898078e9, &
    2.940378956634553899906876e10, &
    1.702665737765398868392998e11, &
    4.926125793377430887588120e11, &
    5.606251856223951465078242e11 /)
  real, parameter, dimension ( 8 ) :: q1 = (/ &
    6.748212550303777196073036e1, &
    1.113332393857199323513008e3, &
    7.738757056935398733233834e3, &
    2.763987074403340708898585e4, &
    5.499310206226157329794414e4, &
    6.161122180066002127833352e4, &
    3.635127591501940507276287e4, &
    8.785536302431013170870835e3 /)
  real, parameter, dimension ( 8 ) :: q2 = (/ &
    1.830328399370592604055942e2, &
    7.765049321445005871323047e3, &
    1.331903827966074194402448e5, &
    1.136705821321969608938755e6, &
    5.267964117437946917577538e6, &
    1.346701454311101692290052e7, &
    1.782736530353274213975932e7, &
    9.533095591844353613395747e6 /)
  real, parameter, dimension ( 8 ) :: q4 = (/ &
    2.690530175870899333379843e3, &
    6.393885654300092398984238e5, &
    4.135599930241388052042842e7, &
    1.120872109616147941376570e9, &
    1.488613728678813811542398e10, &
    1.016803586272438228077304e11, &
    3.417476345507377132798597e11, &
    4.463158187419713286462081e11 /)
  real res
  real x
  real xden
  real xm1
  real xm2
  real xm4
  real xnum
  real xsq
!
!  Return immediately if the argument is out of range.
!
  if ( x <= 0.0 .or. x > XBIG ) then
    gamma_log = XINF
    return
  end if

  if ( x <= EPS ) then

    res = - log ( x )

  else if ( x <= 1.5 ) then

    if ( x < PNT68 ) then
      corr = - log ( x )
      xm1 = x
    else
      corr = 0.0
      xm1 = ( x - 0.5 ) - 0.5
    end if

    if ( x <= 0.5 .or. x >= PNT68 ) then

      xden = 1.0
      xnum = 0.0

      do i = 1, 8
        xnum = xnum * xm1 + p1(i)
        xden = xden * xm1 + q1(i)
      end do

      res = corr + ( xm1 * ( d1 + xm1 * ( xnum / xden ) ) )

    else

      xm2 = ( x - 0.5 ) - 0.5
      xden = 1.0
      xnum = 0.0
      do i = 1, 8
        xnum = xnum * xm2 + p2(i)
        xden = xden * xm2 + q2(i)
      end do

      res = corr + xm2 * ( d2 + xm2 * ( xnum / xden ) )

    end if

  else if ( x <= 4.0 ) then

    xm2 = x - 2.0
    xden = 1.0
    xnum = 0.0
    do i = 1, 8
      xnum = xnum * xm2 + p2(i)
      xden = xden * xm2 + q2(i)
    end do

    res = xm2 * ( d2 + xm2 * ( xnum / xden ) )

  else if ( x <= 12.0 ) then

    xm4 = x - 4.0
    xden = - 1.0
    xnum = 0.0
    do i = 1, 8
      xnum = xnum * xm4 + p4(i)
      xden = xden * xm4 + q4(i)
    end do

    res = d4 + xm4 * ( xnum / xden )

  else

    res = 0.0

    if ( x <= FRTBIG ) then

      res = c(7)
      xsq = x * x

      do i = 1, 6
        res = res / xsq + c(i)
      end do

    end if

    res = res / x
    corr = log ( x )
    res = res + SQRTPI - 0.5 * corr
    res = res + x * ( corr - 1.0 )

  end if

  gamma_log = res

  return
end
subroutine get_seed ( iseed )
!
!*******************************************************************************
!
!! GET_SEED returns a random seed for the random number generator.
!
!
!  Discussion:
!
!    The seed depends on the current time, and ought to be (slightly)
!    different every millisecond.  Once the seed is obtained, a random
!    number generator should be called a few times to further process
!    the seed.
!
!  Modified:
!
!    09 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer ISEED, a pseudorandom seed value.
!
  implicit none
!
  integer, parameter :: I_MAX = 2147483647
!
  character ( len = 8 ) date
  integer iseed
  double precision temp
  character ( len = 10 ) time
  integer values(8)
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  temp = 0.0

  temp = temp + dble ( values(2) - 1 ) / 11.0
  temp = temp + dble ( values(3) - 1 ) / 30.0
  temp = temp + dble ( values(5) ) / 23.0
  temp = temp + dble ( values(6) ) / 59.0
  temp = temp + dble ( values(7) ) / 59.0
  temp = temp + dble ( values(8) ) / 999.0
  temp = temp / 6.0

  if ( temp <= 0.0 ) then
    temp = 1.0 / 3.0
  else if ( temp >= 1.0 ) then
    temp = 2.0 / 3.0
  end if

  iseed = int ( dble ( I_MAX ) * temp )
!
!  Never use a seed of 0 or I_MAX.
!
  if ( iseed == 0 ) then
    iseed = 1
  end if

  if ( iseed == I_MAX ) then
    iseed = I_MAX - 1
  end if

  return
end
subroutine r_random ( r, rlo, rhi, iseed )
!
!*******************************************************************************
!
!! R_RANDOM returns a random real in a given range.
!
!
!  Modified:
!
!    26 October 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real R, the randomly chosen value.
!
!    Input, real RLO, RHI, the minimum and maximum values.
!
!    Input/output, integer ISEED, a seed for the random number generator.
!
  implicit none
!
  integer iseed
  real r
  real rhi
  real rlo
  real t
  real uniform_01_sample
!
!  Pick a random number in (0,1).
!
  t = uniform_01_sample ( iseed )
!
!  Set R.
!
  r = ( 1.0 - t ) * rlo + t * rhi

  return
end
function uniform_01_sample ( iseed )
!
!*******************************************************************************
!
!! UNIFORM_01_SAMPLE is a portable random number generator.
!
!
!  Formula:
!
!    ISEED = ISEED * (7**5) mod (2**31 - 1)
!    UNIFORM_01_SAMPLE = ISEED * / ( 2**31 - 1 )
!
!  Parameters:
!
!    Input/output, integer ISEED, the integer "seed" used to generate
!    the output random number, and updated in preparation for the
!    next one.  ISEED should not be zero.
!
!    Output, real UNIFORM_01_SAMPLE, a random value between 0 and 1.
!
!  Local Parameters:
!
!    IA = 7**5
!    IB = 2**15
!    IB16 = 2**16
!    IP = 2**31-1
!
  implicit none
!
  integer, parameter :: ia = 16807
  integer, parameter :: ib15 = 32768
  integer, parameter :: ib16 = 65536
  integer, parameter :: ip = 2147483647
  integer iprhi
  integer iseed
  integer ixhi
  integer k
  integer leftlo
  integer loxa
  real uniform_01_sample
!
!  Don't let ISEED be 0.
!
  if ( iseed == 0 ) then
    iseed = ip
    end if
!
!  Get the 15 high order bits of ISEED.
!
  ixhi = iseed / ib16
!
!  Get the 16 low bits of ISEED and form the low product.
!
  loxa = ( iseed - ixhi * ib16 ) * ia
!
!  Get the 15 high order bits of the low product.
!
  leftlo = loxa / ib16
!
!  Form the 31 highest bits of the full product.
!
  iprhi = ixhi * ia + leftlo
!
!  Get overflow past the 31st bit of full product.
!
  k = iprhi / ib15
!
!  Assemble all the parts and presubtract IP.  The parentheses are
!  essential.
!
  iseed = ( ( ( loxa - leftlo * ib16 ) - ip ) + ( iprhi - k * ib15 ) * ib16 ) &
    + k
!
!  Add IP back in if necessary.
!
  if ( iseed < 0 ) then
    iseed = iseed + ip
    end if
!
!  Multiply by 1 / (2**31-1).
!
  uniform_01_sample = real ( iseed ) * 4.656612875e-10

  return
end