How to calculate the mode of a probability distribution

Question

I was wondering if there is any command to get the mode of the probability distributions. Calculating manually, the mode is (alpha-1)/(alpha+beta-2)

While the mean, median or variance of the beta dist can be found. For example, the mean can be found using this following command:

Mean[BetaDistribution[α, β]]

As suggested by @kguler , this command ArgMax[PDF[BetaDistribution[2, 3], x], x] works perfectly fine. However, when I changed the value of the parameters to something like thousands, the running process took a very long time.

Is there any other command that will produce the results much faster?

Cheers.

Answer 1

Experience working with distributions suggests analyzing the logarithm of the density function, rather than the density itself. Because the log is a monotonic increasing transformation, the mode of the log density occurs at the same value as the mode of the density. (This approach has general application, not just for beta distributions.)

Let's develop this solution in steps, asking Mathematica to simplify as much as possible at each stage. First, the distribution itself:

f[x_, {a_, b_}] := PDF[BetaDistribution[a, b]][x]

(Bear in mind that this is defined only for positive values of $a$ and $b$.)

Now, its logarithm. Take care to avoid places where the distribution will be zero or undefined:

g[x_, {a_, b_}] := Evaluate[Simplify[Log[f[x, {a, b}]], Assumptions -> 0 < x < 1]]

We plan to inspect the zeros of its derivative for the mode:

h[x_, {a_, b_}] := Evaluate[Simplify[D[g[x, {a, b}], {x}]]]

The critical points will include all zeros of the derivative:

criticalPoints = Solve[h[x, {a, b}] == 0, {x}]

$\left\{\left\{x\to \frac{-1+a}{-2+a+b}\right\}\right\}$

---

If we're working blindly with an unfamiliar distribution, we had better check the endpoints to make sure that one of them isn't actually a mode, too. This can be problematic, because Mathematica often will have trouble finding the limiting values for symbolic parameters $a$ and $b$. One alternative is to explore the situation graphically and manipulate $a$ and $b$ over reasonable ranges:

Manipulate[
 TableForm[{"Value at 0: " <> ToString[Limit[f[x, {a, b}], x -> 0, Direction -> -1]], 
   "Value at 1: " <> ToString[Limit[f[x, {a, b}], x -> 1, Direction -> 1]], 
   Plot[f[x, {a, b}], {x, 0, 1}, ImageSize -> 400]}], 
     {{a, 1/2}, 0, 2}, {{b, 3/2}, 0, 2}]

This will quickly show that

• When $a \lt 1$, there is a mode at $0$.

• When $b \lt 1$, there is a mode at $1$.

• Otherwise, when both $a \ge 1$ and $b \ge 1$, the mode is at $\frac{a-1}{[a-1] + [b-1]}$, as given by the value of criticalPoints.

When both $a \lt 1$ and $b \lt 1$, the distribution is bimodal ("U shaped").

Answer 2

You can just calculate the derivative of the corresponding PDF and then set it to zero; for generic α, β this will not give a very nice result without additional manipulations (however the general solution is there).

Reduce[Simplify[D[PDF[BetaDistribution[α, β], x], x], 
   Assumptions -> {0 < x < 1}] == 0, x]

For numeric values of the distribution parameters you get a much cleaner answer in most cases :

Reduce[Simplify[D[PDF[BetaDistribution[1/2, 1/3], x], x], 
 Assumptions -> {0 < x < 1}] == 0 , x]
(* x == 3/7 *) 

For numeric values however I'd suggest using the answer provided by @kguler (for α = 1, β = 2 my method fails because the PDF is a linear function).