anchor.diff:  
       dlpna = dlpna + wsum
   
  -    logw = log(fpa) + lpna
  +    logw = log(alpha_prior(alpha)) + lpna
       dlogw = dfpa / fpa + dlpna
   
     end subroutine log_weight_d
  
positive.code:     use gamma_funcs, only: digamma, trigamma, quadgamma
      use dirichlet_mod, only: alphabet_size, n_empty_bins, n_data, multi,&
           multi_z
      use dirichlet_mod, only: log_pna, log_fpa
  
      real(real64), intent(in) :: alpha
      real(real64), intent(out) :: logw, dlogw
  
      real(real64) :: fpa, dfpa, lpna, dlpna, wsum
  
      fpa = alphabet_size * trigamma(alphabet_size * alpha + one) - &
           trigamma(alpha + one)
  
      dfpa = alphabet_size**2 * quadgamma(alphabet_size * alpha + one) - &
           quadgamma(alpha + one)
  
negative.code: 
    elemental real(real64) function log_weight(alpha)
      ! un-normalized weight for alpha in the integrals; prop. to p(alpha|x)
      use dirichlet_mod, only: log_pna, alpha_prior
  
      real(real64), intent(in) :: alpha
  
      log_weight = log(alpha_prior(alpha)) + log_pna(alpha)
  
    end function log_weight
  
  