# How to calculate the mode of a probability distribution

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.

-
Welcome to Mathematica.SE! It might be easier for everyone if you pick a more distinctive user name. As to your question, Commonest gives the mode of a list, so you could do this on a sample, but I assume what you really want is to get the maximum of the PDF by taking the derivative and solving for zero? – Verbeia Oct 11 '12 at 12:13
hi, thanks for the suggestion! I will change my username soon :) I have tried use commonest but here is what I got: Commonest::arg1: The first argument is expected to be a list. – yyasinta Oct 11 '12 at 12:16
For numeric parameters, you can use mode=ArgMax[PDF[your distribution, x],x]: for example,ArgMax[PDF[BetaDistribution[2, 3], x], x] gives 1/3. – kglr Oct 11 '12 at 12:21
my only problem with this command is that it takes ages in my Mathematica to run it with very large numeric parameters value. :( – yyasinta Oct 11 '12 at 12:57

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").

-
@whubber, thank you for such a detail and useful explanation. However, where should I input my numeric parameters into? I have tried to input the value into last equation and when I define the dist, but it doesnt work. – yyasinta Oct 12 '12 at 0:26
They are a and b. As an example of one way to obtain the answer for, say, $a=10000$ and $b=15000$, run the steps up to obtaining criticalPoints and then use criticalPoints /. {a -> 10000, b -> 15000}. The whole point is that Mathematica has given you this simple formula for the mode, so once you have it you might as well use it. – whuber Oct 12 '12 at 13:56

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).

-
will this work if I substitute the parameters with numerical values? Cause I tried to use your suggestion with alpha=1 and beta=2 but then it gives me this message:PolynomialGCD::lrgexp: Exponent is out of bounds for function PolynomialGCD. – yyasinta Oct 11 '12 at 12:23