anchor.diff:           - alphabet_size * log_gamma(alpha) &
            - log_gamma(n_data + alpha * alphabet_size)
   
  -    wsum = sum(multi * (log_gamma(multi_z + alpha) - log_gamma(multi_z + one)))
  +    wsum = n_empty_bins * (log_gamma(alpha) - log_gamma(one))
  +    wsum = wsum + &
  +         sum(multi * (log_gamma(multi_z + alpha) - log_gamma(multi_z + one)))
   
       log_pna = log_pna + wsum
       
  
positive.code:     allocate(multi_z(n_multi),stat=err)
      allocate(multi(n_multi),stat=err)
      k_ = 0
      do i_ = 1, nmax + 1
         if (multi0(i_) > 0) then 
            k_ = k_ + 1
            multi_z(k_) = i_ - 1
            multi(k_) = multi0(i_)
         end if
      end do
  
negative.code:          - alphabet_size * log_gamma(alpha) &
           - log_gamma(n_data + alpha * alphabet_size)
  
      asum = n_empty_bins * alpha * digamma(alpha + one)
      asum = asum + sum(multi * (multi_z + alpha) &
           * digamma(multi_z + alpha + one))
  
      bsum = n_empty_bins * (log_gamma(alpha) - log_gamma(one))
      bsum = bsum + sum(multi * (log_gamma(multi_z + alpha) &
           - log_gamma(multi_z + one)))
  
      lpna = lpna + bsum
  