# How do I compute the entropy of the beta distribution?

I tried

Expectation[-q*Log[q], q \[Distributed] BetaDistribution[a, b]]


and got

(a (HarmonicNumber[a] - HarmonicNumber[a + b]))/(a + b)


(I get the same result using integration.) This seems pretty different from the expression in Wikipedia, even after substituting the identity (mentioned at MathWorld that H[n] = gamma + PolyGamma[0,n+1], where gamma is Euler's constant and PolyGamma[0, n+1] = Digamma[n+1].

Any ideas which identity I'm missing?

Kevin

P.S. I trust the Wikipedia answer, I have seen this result before, but I want to re-derive it with Mathematica as a sanity check, before I go on to derive other expressions of interest.

-
You could plot the Wikipedia result and yours for a number of parameter values to see whether they might be the same after all. –  Sjoerd C. de Vries Nov 6 '13 at 7:02

Your definition of entropy is incorrect. It's $E(-\ln(P(x)))$, with $E$ the expectation operator and $P$ the probability mass function of the random variable $x$. I believe you may have been mixing up a few things.

The formal definition of the expectation is $E(x)=\int{x P(x)dx}$. I assume that you had this in mind and you further confused your random variable $q$ with its PDF to write

Expectation[-q*Log[q], q \[Distributed] BetaDistribution[a, b]]


for the entropy.

Expectation[-Log[PDF[BetaDistribution[a, b], q]], q \[Distributed] BetaDistribution[a, b]]


Unfortunately, I couldn't prod Mathematica to show the equality between Wikipedia's result and this formal definition symbolically. However, numerically it works out quite well.

This is Wikipedia's version:

In Mathematica terms:

Ψ[n_] := HarmonicNumber[n - 1] - EulerGamma

entropyBeta[a_, b_] =
FullSimplify[Log[Beta[a, b]] - (a - 1) Ψ[a] - (b - 1) Ψ[b] + (a + b - 2) Ψ[a + b],
Assumptions -> {a > 0, b > 0}];


Comparing:

Manipulate[
Show[
Plot[entropyBeta[a, b], {b, 0.01, 5}],
DiscretePlot[
Quiet@