subroutine bbl ( mu, theta, rl, mrl, lm, rnl )

!*****************************************************************************80
!
!! BBL calculates the beta binomial log likelihood.
!
!  Modified:
!
!    30 March 1999
!
!  Author:
!
!    D Smith
!
!  Reference:
!
!    D Smith,
!    Algorithm AS 189,
!    Maximum Likelihood Estimation of the Parameters of the Beta
!    Binomial Distribution,
!    Applied Statistics, 
!    Volume 32, Number 2, 1983, pages 196-204.
!
!  Parameters:
!
!    Input, real MU, the estimated value of MU.
!
!    Input, real THETA, the estimated value of THETA.
!
!    Input, integer RL(MRL,3), array of coefficients of 
!    (MU + R * THETA), ( 1 - MU + R * THETA) and ( 1 + R * THETA ) terms.
!
!    Input, integer MRL, the first dimension of the RL array, 
!    which must be at least the maximum of the values in IN(*).
!
!    Input, integer LM(3), contain the values Max ( IX(J) - 1 ),
!    Max ( IN(J) - IX(J) - 1 ), and Max ( IN(J) - 1 ).
!
!    Output, real RNL, the log likelihood. 
!
  implicit none

  integer mrl

  real a
  integer i
  integer lm(3)
  real mu
  integer mlm
  integer rl(mrl,3)
  real rnl
  real theta

  rnl = 0.0E+00
  mlm = lm(3)

  do i = 1, mlm

    a = real ( i - 1 ) * theta

    if ( i <= lm(1) ) then
      rnl = rnl + real ( rl(i,1) ) * log ( mu + a )
    end if

    if ( i <= lm(2) ) then
      rnl = rnl + real ( rl(i,2) ) * log ( 1.0E+00 - mu + a )
    end if

    rnl = rnl - real ( rl(i,3) ) * log ( 1.0E+00 + a )

  end do

  return
end
subroutine bbme ( n, ix, in, inf, mu, theta )

!*****************************************************************************80
!
!! BBME estimates MU and THETA of the beta binomial distribution.
!
!  Discussion:
!
!    The method of moments is used.
!
!  Modified:
!
!    02 February 2003
!
!  Author:
!
!    D Smith
!
!  Reference:
!
!    D Smith,
!    Algorithm AS 189,
!    Maximum Likelihood Estimation of the Parameters of the Beta
!    Binomial Distribution,
!    Applied Statistics, 
!    Volume 32, Number 2, 1983, pages 196-204.
!
!  Parameters:
!
!    Input, integer N, the number of observations or trials.
!
!    Input, integer IX(N), contains the number of successes for 
!    each trial.
!
!    Input, integer IN(N), contains the number tested on each trial.
!
!    Input, real INF, the maximum acceptable value for THETA.
!
!    Output, real MU, the estimate for MU.
!
!    Output, real THETA, the estimate for THETA.
!
  implicit none

  integer n

  real d1
  real d2
  integer i
  integer in(n)
  real inf
  integer ix(n)
  logical j
  real mu
  real p(n)
  real r
  real s
  real theta
  real tp
  real w(n)
  real w_sum

  j = .false.
  w(1:n) = real ( in(1:n) )
  p(1:n) = real ( ix(1:n) ) / real ( in(1:n) )

  do

    w_sum = sum ( w(1:n) )
    tp = dot_product ( w(1:n), p(1:n) )

    tp = tp / w_sum

    s = 0.0E+00
    d1 = 0.0E+00
    d2 = 0.0E+00
    do i = 1, n
      r = p(i) - tp
      s = s + w(i) * r**2
      r = w(i) * ( 1.0E+00 - w(i) / w_sum )
      d1 = d1 + r / real ( in(i) )
      d2 = d2 + r
    end do

    s = real ( n - 1 ) * s / real ( n )
    r = tp * ( 1.0E+00 - tp )

    if ( r == 0.0E+00 ) then
      exit
    end if

    r = ( s - r * d1 ) / ( r * ( d2 - d1 ) )
    r = max ( r, 0.0E+00 )

    if ( j ) then
      exit
    end if

    w(1:n) = w(1:n) / ( 1.0E+00 + r * ( w(1:n) - 1.0E+00 ) )

    j = .true.

  end do
!
!  Set the estimates.
!
  mu = tp

  if ( r < 1.0E+00 ) then

    theta = r / ( 1.0E+00 - r )

    theta = min ( theta, inf )

  else
 
    theta = inf

  end if

  return
end
subroutine bbml ( n, ix, in, rl, mrl, iter, ccrit, mu, theta, mu_se, theta_se, &
  rnl, ifault )

!*****************************************************************************80
!
!! BBML estimates the parameters of a beta binomial distribution
!
!  Definition:
!
!    The beta binomial probability density function for X successes
!    out of N trials is
!
!      PDF(X) ( N, MU, THETA ) =
!        C(N,X) * Product ( 0 <= R <= X - 1 ) ( MU + R * THETA )
!               * Product ( 0 <= R <= N - X - 1 ) ( 1 - MU + R * THETA )
!               / Product ( 0 <= R <= N - 1 )  ( 1 + R * THETA )
!
!    where 
!
!      C(N,X) is the combinatorial coefficient;
!      MU is the expectation of the underlying Beta distribution;
!      THETA is a shape parameter.  
!
!    A THETA value of 0 results in a PDF equivalent to the binomial 
!    distribution:
!
!      PDF(X) ( N, MU, 0 ) = C(N,X) * MU**X * ( 1 - MU )**(N-X)
!
!    This PDF can be reformulated as:
!
!      PDF2(X)(A,B,C) = Beta(A+X,B+C-X) 
!        / ( (C+1) * Beta(X+1,C-X+1) * Beta(A,B) )  for 0 <= X <= C.
!
!    Given A, B, C for PDF2, the equivalent PDF has:
!
!      N     = C
!      MU    = A / ( A + B )
!      THETA = 1 / ( A + B )
!
!    Given N, MU, THETA for PDF, the equivalent PDF2 has:
!
!      A = MU / THETA
!      B = ( 1 - MU ) / THETA
!      C = N
!
!  Modified:
!
!    30 March 1999
!
!  Author:
!
!    D Smith
!
!  Reference:
!
!    D Smith,
!    Algorithm AS 189,
!    Maximum Likelihood Estimation of the Parameters of the Beta
!    Binomial Distribution,
!    Applied Statistics, 
!    Volume 32, Number 2, 1983, pages 196-204.
!
!  Parameters:
!
!    Input, integer N, the number of observations or trials.
!
!    Input, integer IX(N), contains the number of successes for 
!    each trial.
!
!    Input, integer IN(N), contains the number tested on each trial.
!
!    Workspace, integer RL(MRL,3), array of coefficients of 
!    (MU + R * THETA), ( 1 - MU + R * THETA) and ( 1 + R * THETA ) terms.
!
!    Input, integer MRL, the first dimension of the RL array, 
!    which must be at least the maximum of the values in IN(*).
!
!    Input/output, integer ITER;
!    On input, the maximum number of iterations allowed.
!    On output, the number of iterations taken.
!    
!    Input, real CCRIT, the convergence criterion.  The iteration is 
!    judged to have converged when abs ( delta MU ) and 
!    abs ( delta THETA) are less than or equal to CCRIT.
!
!    Output, real MU, the maximum likelihood estimate of MU, the mean
!    of the beta binomial distribution.
!
!    Output, real THETA, the maximum likelihood estimate of THETA, the
!    shape parameter of the beta binomial distribution.
!
!    Output, real MU_SE, the standard error of the estimate of MU; returned
!    as -1.0 if it cannot be calculated.
!
!    Output, real THETA_SE, the standard error of the estimate of THETA;
!    returned as -1.0 if it cannot be calculated.
!
!    Output, real RNL, the log likelihood for the maximum likelihood
!    estimates.
!
!    Output, integer IFAULT, error flag.
!    0, no error.
!    1, N <= 1;
!    2, IX(I) = 0 for all I;
!    3, IX(I) = IN(I) for all I;
!    4, max ( IN(I) ) > MRL;
!    5, either IX(I) < 0 or IN(I) < IX(I) for some I;
!    6, MU went outside the range of [0,1], or THETA went outside the
!       range [0,INF], where INF represents Infinity;
!    7, if the maximum number of iterations was exceeded;
!    8, if the damped Newton-Raphson iteration failed.
!
  implicit none

  integer n

  real a
  real b
  real c
  real ccrit
  real d
  real del
  real dum
  real e
  real eps
  real f
  real fd(2)
  integer ifault
  integer in(n)
  real, parameter :: inf = 1.0E+06
  integer iter
  integer iter_max
  integer ix(n)
  integer lm(3)
  logical mc
  real mu
  real mu_se
  integer mrl
  integer nnd
  integer rd1(2,2)
  integer rd2(2,3)
  integer rd3(2,4)
  integer rl(mrl,3)
  real rnl
  real sd(3)
  real td(4)
  real theta
  real theta_se
  real ub(2)
!
  rd1(1,1) =   1
  rd1(2,1) = - 1
  rd1(1,2) =   1
  rd1(2,2) =   1

  rd2(1,1) = - 1
  rd2(2,1) = - 1
  rd2(1,2) = - 1
  rd2(2,2) =   1
  rd2(1,3) = - 1
  rd2(2,3) = - 1

  rd3(1,1) =   2
  rd3(2,1) = - 2
  rd3(1,2) =   2
  rd3(2,2) =   2
  rd3(1,3) =   2
  rd3(2,3) = - 2
  rd3(1,4) =   2
  rd3(2,4) =   2

  iter_max = iter
  iter = 0
  mc = .true.
  ub(1) = 0.01E+00
  ub(2) = 0.01E+00
!
!  Set the arrays RL and LM.
!
  call set ( n, ix, in, rl, mrl, lm, ifault )

  if ( ifault /= 0 ) then
    return
  end if

  mu_se = - 1.0E+00
  theta_se = - 1.0E+00
  nnd = 0
!
!  Calculate initial estimates by the method of moments.
!
  call bbme ( n, ix, in, inf, mu, theta )

  if ( theta == inf ) then
    go to 50
  end if
!
!  Newton-Raphson iteration on first derivatives.
!
5 continue

  if ( iter_max < iter ) then
    ifault = 7
    go to 60
  end if
!
!  Calculate the first derivatives of the log likelihood.
!
  call gder ( mu, theta, rl, mrl, lm, 2, rd1, fd )
!
!  Calculate the second derivatives of the log likelihood.
!
  call gder ( mu, theta, rl, mrl, lm, 3, rd2, sd )
!
!  Calculate the third derivatives of the log likelihood.
!
  call gder ( mu, theta, rl, mrl, lm, 4, rd3, td )
!
!  Calculate the increments.
!
  dum = sd(1) * sd(3) - sd(2)**2

  if ( sd(1) < 0.0E+00 .and. dum > 0.0E+00 ) then
    go to 15
  end if
!
!  Non negative definite matrix.
!
  nnd = nnd + 1
!
!  SD(1) is always negative so a gradient step is made on MU.
!
  a = mu - fd(1) / sd(1)
  b = theta

  if ( fd(2) /= 0.0E+00 ) then
    b = b + sign ( ub(2), fd(2) )
  end if

  if ( a <= 0.0E+00 ) then
    a = 0.0001E+00
  else if ( a >= 1.0E+00 ) then
    a = 0.9999E+00
  end if

  b = max ( b, 0.0E+00 )
  b = min ( b, inf )

  call bbl ( mu, theta, rl, mrl, lm, c )
  call bbl ( a, b, rl, mrl, lm, d )

  if ( nnd > 10 .or. c >= d ) then
    go to 40
  end if

  iter = iter + 1
  mu = a
  theta = b
  go to 5

15 continue

  del = ( fd(2) * sd(2) - fd(1) * sd(3) ) / dum
  eps = ( fd(1) * sd(2) - fd(2) * sd(1) ) / dum
!
!  Check to see if the Lipschitz condition is satisfied.
!
  a = sd(2) * td(2) - td(1) * sd(3)
  b = sd(2) * td(3) - td(2) * sd(3)
  c = td(1) * sd(2) - td(2) * sd(1)
  d = sd(2) * td(2) - sd(1) * td(3)
  e = sd(2) * td(4) - td(3) * sd(3)
  f = td(3) * sd(2) - td(4) * sd(1)

  a = del * a + eps * b
  c = del * c + eps * d
  e = del * b + eps * e
  f = del * d + eps * f

  dum = ( a**2 + c**2 + e**2 + f**2 ) / dum**2

  if ( dum >= 1.0E+00 ) then
    go to 20
  end if

  if ( abs ( del ) <= ccrit .and. abs ( eps ) <= ccrit ) then
    mc = .false.
  end if

  go to 45
!
!  Failure of the Lipschitz condition. 
!  A step in the direction of the gradient is made.
!
20 continue

  a = fd(1)**2
  b = fd(2)**2
  c = a * sd(1) + 2.0E+00 * sd(2) * fd(1) * fd(2) + b * sd(3)

  if ( c /= 0.0E+00 ) then

    c = - ( a + b ) / c
    del = c * fd(1)
    eps = c * fd(2)

    if ( abs ( del ) > ub(1) ) then
      del = sign ( ub(1), del )
    end if

    ub(1) = 2.0E+00 * abs ( del )

    if ( abs ( eps ) > ub ( 2 ) ) then
      eps = sign ( ub(2), eps )
    end if

    ub(2) = 2.0E+00 * abs ( eps )

  else

    if ( fd(1) /= 0.0E+00 ) then 
      del = sign ( ub(1), fd(1) )
    else
      del = 0.0E+00
    end if

    if ( fd(2) /= 0.0E+00 ) then
      eps = sign ( ub(2), fd(2) )
    else
      eps = 0.0E+00
    end if

  end if

  call bbl ( mu, theta, rl, mrl, lm, c )
!
!  Begin loop.
!
35 continue

  a = mu + del
  b = theta + eps

  if ( a <= 0.0E+00 ) then
    a = 0.0001E+00
  end if

  if ( a >= 1.0E+00 ) then
    a = 0.9999E+00
  end if

  del = a - mu

  b = max ( b, 0.0E+00 )
  b = min ( b, inf )
  
  eps = b - theta
  call bbl ( a, b, rl, mrl, lm, d )
!
!  Check to see if gradient step has increased log likelihood.
!
  if ( d > c ) then
    go to 45
  end if

  del = del / 2.0E+00
  eps = eps / 2.0E+00

  if ( abs ( del ) > ccrit .or. abs ( eps ) > ccrit ) then
    go to 35
  end if

40    continue

  ifault = 8
  go to 60

45    continue

  iter = iter + 1
  a = mu + del
  b = theta + eps

  if ( a > 0.0E+00 .and. a < 1.0E+00 .and. b >= 0.0E+00 .and. b <= inf ) then
    go to 55
  end if

  if ( a <= 0.0E+00 ) then
    mu = 0.0E+00
  end if

  if ( a >= 1.0E+00 ) then
    mu = 1.0E+00
  end if

  if ( b < 0.0E+00 ) then
    theta = 0.0E+00
  end if

  if ( b > inf ) then
    theta = inf
  end if

50 continue

  ifault = 6
  go to 60

55 continue

  mu = a
  theta = b

  if ( mc ) then
    go to 5
  end if
!
!  Calculate the standard errors.
!
  if ( sd(1) < 0.0E+00 ) then
    mu_se = sqrt ( - 1.0E+00 / sd(1) )
  end if

  if ( sd(3) < 0.0E+00 ) then
    theta_se = sqrt ( - 1.0E+00 / sd(3) )
  end if
!
!  Calculate the log likelihood.
!
60 continue

  call bbl ( mu, theta, rl, mrl, lm, rnl )

  return
end
subroutine gder ( mu, theta, rl, mrl, lm, ider, rd, pd )

!*****************************************************************************80
!
!! GDER computes derivatives of the log likelihood.
!
!  Modified:
!
!    30 March 1999
!
!  Author:
!
!    D Smith
!
!  Reference:
!
!    D Smith,
!    Algorithm AS 189,
!    Maximum Likelihood Estimation of the Parameters of the Beta
!    Binomial Distribution,
!    Applied Statistics, 
!    Volume 32, Number 2, 1983, pages 196-204.
!
!  Parameters:
!
!    Input, real MU, the estimated value of MU.
!
!    Input, real THETA, the estimated value of THETA.
!
!    Input, integer RL(MRL,3), array of coefficients of 
!    (MU + R * THETA), ( 1 - MU + R * THETA) and ( 1 + R * THETA ) terms.
!
!    Input, integer MRL, the first dimension of the RL array, 
!    which must be at least the maximum of the values in IN(*).
!
!    Input, integer LM(3), contain the values Max ( IX(J) - 1 ),
!    Max ( IN(J) - IX(J) - 1 ), and Max ( IN(J) - 1 ).
!
!    Input, integer IDER, 1 plus the order of the derivative
!    desired.  IDER can be 2, 3 or 4.
!
!    Input, integer RD(2,IDER), an array of coefficients.
!
!    Output, real PD(IDER), the derivatives of the log likelihood. 
!
  implicit none

  integer ider
  integer mrl

  real a
  real b
  real c
  real d
  integer i
  integer j
  integer k
  integer kk
  integer lm(3)
  real mu
  integer mlm
  real pd(ider)
  integer rd(2,ider)
  integer rl(mrl,3)
  real theta

  mlm = lm(3)
  kk = ider - 1

  pd(1:ider) = 0.0E+00

  do i = 1, mlm

    c = real ( i - 1 )
    a = c * theta

    do j = 1, 3

      if ( i <= lm(j) ) then

        if ( j == 1 ) then
          d = mu + a
        else if ( j == 2 ) then
          d = 1.0E+00 - mu + a
        else if ( j == 3 ) then
          d = 1.0E+00 + a
        end if

        b = real ( rl(i,j) ) / d**kk

        if ( j /= 3 ) then

          do k = 1, ider
            pd(k) = pd(k) + real ( rd(j,k) ) * b
            b = b * c
          end do

        else

          d = - real ( rd(1,1) ) * b * c**kk
          pd(ider) = pd(ider) + d

        end if

      end if

    end do

  end do

  return
end
subroutine set ( n, ix, in, rl, mrl, lm, ifault )

!*****************************************************************************80
!
!! SET sets up the arrays RL and LM.
!
!  Modified:
!
!    30 March 1999
!
!  Author:
!
!    D Smith
!
!  Reference:
!
!    D Smith,
!    Algorithm AS 189,
!    Maximum Likelihood Estimation of the Parameters of the Beta
!    Binomial Distribution,
!    Applied Statistics, 
!    Volume 32, Number 2, 1983, pages 196-204.
!
!  Parameters:
!
!    Input, integer N, the number of observations or trials.
!
!    Input, integer IX(N), contains the number of successes for 
!    each trial.
!
!    Input, integer IN(N), contains the number tested on each trial.
!
!    Output, integer RL(MRL,3), array of coefficients of 
!    (MU + R * THETA), ( 1 - MU + R * THETA) and ( 1 + R * THETA ) terms.
!
!    Input, integer MRL, the first dimension of the RL array, 
!    which must be at least the maximum of the values in IN(*).
!
!    Output, integer LM(3), contain the values Max ( IX(J) - 1 ),
!    Max ( IN(J) - IX(J) - 1 ), and Max ( IN(J) - 1 ).
!
!    Output, integer IFAULT, error flag.
!    0, no error.
!    1, N <= 1;
!    2, IX(I) = 0 for all I;
!    3, IX(I) = IN(I) for all I;
!    4, max ( IN(I) ) > MRL;
!    5, either IX(I) < 0 or IN(I) < IX(I) for some I.
!
  implicit none

  integer mrl
  integer n

  integer i
  integer ifault
  integer in(n)
  integer ix(n)
  integer j
  integer jj
  integer k
  integer lm(3)
  integer mar
  integer rl(mrl,3)
!
!  Check the input data.
!
  if ( n <= 1 ) then
    ifault = 1
    return
  end if

  ifault = 2
  do i = 1, n
    if ( ix(i) > 0 ) then
      ifault = 0
    end if
  end do

  if ( ifault == 2 ) then
    return
  end if

  ifault = 3
  do i = 1, n
    if ( ix(i) < in(i) ) then
      ifault = 0
    end if
  end do

  if ( ifault == 3 ) then
    return
  end if
!
!  Form the matrix of counts.
!
  ifault = 4

  lm(1:3) = 0
  rl(1:mrl,1:3) = 0

  do i = 1, n

    jj = ix(i)
    mar = 1
 
10  continue

    if ( jj < 0 ) then
      ifault = 5
      return
    end if

    if ( jj > 0 ) then

      if ( jj > mrl ) then
        return
      end if

      if ( jj > lm(mar) ) then
        lm(mar) = jj
      end if

      rl(jj,mar) = rl(jj,mar) + 1

    end if

    if ( mar == 1 ) then
      jj = in(i) - ix(i)
      mar = 2
      go to 10
    else if ( mar == 2 ) then
      jj = in(i)
      mar = 3
      go to 10
    end if

  end do

  ifault = 0
!
!  Evaluate number of calls to different terms of likelihood function.
!
  do i = 1, 3

    jj = lm(i) - 1
    k = jj
    do j = 1, jj
      rl(k,i) = rl(k,i) + rl(k+1,i)
      k = k - 1
    end do

  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end