anchor.diff:    call initialize_dirichlet(counts, nc)
     call compute_multiplicities(counts)
   
  -  estimate = hdir(alpha)
  +  estimate = h_bayes(alpha)
   
     call dirichlet_finalize()
   
  
positive.code:     integer(int32) :: i
      real(real64)   :: a(n_multi)
  
      hdir = 0.0_real64
      do i = 1,n_multi 
         a(i) = multi(i) * (multi_z(i) + alpha) * digamma(multi_z(i) + alpha + one) 
      end do
      hdir = - sum(a)
      hdir = hdir / (n_data + alpha * alphabet_size)
      hdir = hdir + digamma(n_data + alpha * alphabet_size + one)
  
    end function hdir
  
  end module dirichlet_mod
  
  module nsb_mod
    use iso_fortran_env
    implicit none
  
    real(real64) :: log_alpha1
    real(real64) :: log_alpha2
    real(real64) :: amax
    real(real64) :: log_fwa_amax
  
  contains
  
    elemental real(real64) function log_fpa(a) 
      ! prop. to p(alpha) - the prior for alpha in NSB estimator
      use constants
      use gamma_funcs, only: trigamma
      use dirichlet_mod, only: alphabet_size
      
      real(real64), intent(in) :: a
      
  
negative.code:     end do
      xs = exp(xs)
      
      fxs = log_weight(xs)
      amax = xs(maxloc(fxs,1,fxs < largest))
      lw_max = log_weight(amax)
      
    end subroutine compute_integration_range
  
    real(real64) function m_func(x)
      use dirichlet_mod, only: h_bayes, integrand
  
      real(real64), intent(in) :: x
      real(real64) :: a, h_c,lw
  
      a = exp(x)
      call integrand(a, h_c, lw)
      m_func = exp(lw - lw_max) * h_c * a
      
    end function m_func
  
    real(real64) function m2_func(x)
      use dirichlet_mod, only: h_bayes, integrand
  
      real(real64), intent(in) :: x
      real(real64) :: a, h_c,lw
  
      a = exp(x)
      call integrand(a, h_c, lw)
      m2_func = exp(lw - lw_max) * h_c**2 * a
      
    end function m2_func
  
    real(real64) function nrm_func(x)
  
      real(real64), intent(in) :: x
      real(real64) :: a
  
      a = exp(x)
      nrm_func = exp(log_weight(a) - lw_max)  * a
      
    end function nrm_func
  
  