anchor.diff:  
     end function log_weight
   
  -  elemental real(real64) function h_bayes(alpha)
  +  elemental real(real64) function h_dir(alpha)
       ! posterior average of the entropy given data and a specific alpha value
       ! computed from histogram multiplicities
       use gamma_funcs, only: digamma
  
positive.code:     real(real64), intent(in) :: alpha
      integer(int32) :: i_
  
      h_bayes = - sum(multi * (multi_z + alpha) * digamma(multi_z + alpha + one))
      h_bayes = h_bayes / (n_data + alpha * alphabet_size)
      h_bayes = h_bayes + digamma(n_data + alpha * alphabet_size + one)
  
    end function h_bayes
  
    elemental real(real64) function integrand(alpha, amax, order)
      ! posterior average of the entropy given the data and alpha
  
negative.code:   call initialize_dirichlet(counts, nc)
    call compute_multiplicities(counts)
  
    estimate = h_dir(alpha)
  
    call dirichlet_finalize()
  
  