Expectation[-q*Log[q], q \[Distributed] BetaDistribution[a, b]]
(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.