anchor.diff:      - alphabet_size * digamma(alpha) &
            - alphabet_size * digamma(n_data + alpha * alphabet_size)
   
  -    wsum = sum(hzs * (digamma(hns + alpha)))
  +    wsum = sum(hz * (digamma(hn + alpha)))
   
       dlpna = dlpna + wsum
   
  
positive.code:     real(real64) :: c, nu, ni, xi, jsum
  
      nu = n_data + alpha * alphabet_size
      phi = digamma(hns + alpha + one) - &
           digamma(nu + two)
      c = trigamma(nu + two)
  
      h_var = 0.0
      do i_ = 0, size(hzs)-1
         ni = hns(i_) + alpha
         xi = phi(i_)
         jsum = sum(hzs * ni * (hns + alpha) * &
              (xi * phi - c))
         h_var = h_var + hzs(i_) * jsum
         h_var = h_var - hzs(i_) * ni**2 * (xi**2 - c)
         xi = xi + 1 / (ni + one)
         h_var = h_var + hzs(i_) * (ni + one) * ni * &
              (xi**2 + trigamma(ni + two) - c)
      end do
  
  
negative.code:     log_pna_u = log_gamma(alpha * alphabet_size) &
           - alphabet_size * log_gamma(alpha) &
           - log_gamma(n_data + alpha * alphabet_size) &
           + sum(hz * (log_gamma(hn + alpha)))
  
    end function log_pna_u
  
  