Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?


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.

share|improve this question
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.

Instead this should have been:

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:

Mathematica graphics

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}];


  Plot[entropyBeta[a, b], {b, 0.01, 5}],
    NExpectation[-Log[PDF[BetaDistribution[a, b], q]], q \[Distributed] BetaDistribution[a, b]]], 
   {b, 0.1, 5, .1}]
 {a, 0.01, 5}]

Mathematica graphics

As you can see, they match very well.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.