'Compute the entropy of the multinomial distribution.
The entropy is computed using this expression:
.. math::
f(x) = - \log n! - n\sum_{i=1}^k p_i \log p_i +
\sum_{i=1}^k \sum_{x=0}^n \binom n x p_i^x(1-p_i)^{n-x} \log x!
Parameters
%(_doc_default_callparams)s
Returns
h : scalar
Entropy of the Multinomial distribution
Notes
%(_doc_callparams_note)s'

def entropy(self, n, p):

    (n, p, npcond) = self._process_parameters(n, p)
    x = np.r_[1:(np.max(n) + 1)]
    term1 = (n * np.sum(entr(p), axis=(-1)))
    term1 -= gammaln((n + 1))
    n = n[..., np.newaxis]
    new_axes_needed = ((max(p.ndim, n.ndim) - x.ndim) + 1)
    x.shape += ((1,) * new_axes_needed)
    term2 = np.sum((binom.pmf(x, n, p) * gammaln((x + 1))), axis=((-1), ((-1) - new_axes_needed)))
    return self._checkresult((term1 + term2), npcond, np.nan)
