# 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@
{b, 0.1, 5, .1}]
],
{a, 0.01, 5}]


As you can see, they match very well.

• One spurious bracket removed. – Sjoerd C. de Vries Oct 23 '19 at 19:34
Clear["Global*"]



The constraints on the parameters are

assume = DistributionParameterAssumptions[dist];


The integral needs a little help:

(entropy[BetaDistribution[a, b]] = Assuming[assume && 0 < x < 1,
Integrate[
-Log[PDF[dist, x]]*PDF[dist, x] // Simplify //
FunctionExpand, {x, 0, 1}] // FullSimplify] /.
Gamma[x_ + y_] :> Gamma[x]*Gamma[y]/Beta[x, y]) // TraditionalForm
`