anchor.diff:  
     end function log_pna
   
  +  pure real(real64) function log_pna_u(alpha)
  +    ! log of "unnormalized" pna. keep only alpha-dependent terms
  +    use constants
  +
  +    real(real64), intent(in) :: alpha
  +
  +    log_pna_u = log_gamma(alpha * alphabet_size) &
  +         - alphabet_size * log_gamma(alpha) &
  +         - log_gamma(n_data + alpha * alphabet_size) &
  +         + sum(multi * (log_gamma(multi_z + alpha)))
  +
  +  end function log_pna_u
  +
  +
     elemental real(real64) function alpha_prior(alpha)
       ! prop. to p(alpha) - the prior for alpha in NSB estimator
       use constants
  
positive.code: 
      ! log weight
      prior = alpha_prior(alpha)
      logw = log(prior) + log_pna(alpha)
  
      ! log weight derivative
      dprior = alphabet_size**2 * quadgamma(alphabet_size * alpha + one) - &
  
negative.code:     ! un-normalized weight for alpha in the integrals; prop. to p(alpha|x)
      real(real64), intent(in) :: alpha
  
      log_weight = log(alpha_prior(alpha)) + log_pna_u(alpha)
  
    end function log_weight
  
  