Find the maximum of a SmoothKernelDistribution

Question

I am trying to find the maximum of a PDF calculated using the SmoothKernelDistribution function. I do so by solving the first order condition but Mathematica seems unable to solve it. As a minimum working example, I've using so far the following:

data = RandomVariate[BetaDistribution[4, 5], 10000];
dist = SmoothKernelDistribution[data];
Solve[D[PDF[dist,x]]==0,x]

The result I obtain is:

Solve::inex: Solve was unable to solve the system with inexact 
coefficients or the system obtained by direct rationalization of 
inexact numbers present in the system. Since many of the methods used 
by Solve require exact input, providing Solve with an exact version of 
the system may help. >>

Moreover, FindRoot proposed as solution {x -> 1.79065}. How should I proceed?

Solve and D are or symbolic computation. Here you are dealing with numerical methods. Use FindMaximum, NMaximize, NArgMax, etc. Take a look at this before you start, in case it will become necessary. Also try FindMaximum[PDF[dist, x], {x, 0.5}]. — Szabolcs, Oct 10 '16 at 13:05
I don't understand the question about FindRoot. It seems to have given you the result. Thus you have the location of the maximum. — Szabolcs, Oct 10 '16 at 13:06
Thanks,@Szabolcs! Regarding the FindRoot it gave an answer but completely off the domain of the function. About the numerical nature of the problem, I thought that the underlying structure of the SmoothKernelDistribution was a spline function. Anyway, I'd like to have an analytical result since this is part of a wider problem. — Keizer, Oct 10 '16 at 13:11
You are right. FindRoot and FindMaximum have options to keep them within the domain, see the documentation. E.g. FindMaximum[..., {x, 0, 1}]. With NMaximize you can either use a constraint, NMaximize[{PDF[dist, x], 0 < x < 1}, x], or a domain: NMaximize[PDF[dist, Indexed[x, 1]], x \[Element] Interval[{0, 1}]]. In the second case x is considered a "1D vector", hence the need to index it. — Szabolcs, Oct 10 '16 at 13:16
About the numerical nature: Yes, the distinction between exact and numerical is sometimes blurred in Mathematica. It is possible to get the exact derivative of an InterpolatingFunction. But Solve cannot work with it. It is possible in theory to convert it to an explicit piecewise polynomial and attempt symbolic manipulation, but it's not straightforward and I don't recommend it. — Szabolcs, Oct 10 '16 at 13:18

score 1 · Accepted Answer · answered Oct 12 '16 at 06:35

Here are a few ways to approximately locate the maximum of the PDF. Let's generate the data and the associated distribution first:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]; (* for reproducibility *)
            data = RandomVariate[BetaDistribution[4, 5], 10000]];

dist = SmoothKernelDistribution[data];
pdf = PDF[dist];

A plot of the PDF shows that it is unimodal, which makes things a lot easier.

In fact, the plotting functionality, through its powerful MeshFunctions option, allows one to approximately locate the maximum:

plt = Plot[pdf[x], {x, 0, 1}, Axes -> None, Frame -> True, Mesh -> {{0}}, 
           MeshFunctions -> {pdf'[#] &}, MeshStyle -> PointSize[Medium]]

plot of the PDF with marked maximum

To find that point, we need to peer into the internal structure of the plot and extract the required coordinates. Here's one way to do this:

pmax = First[Cases[Normal[plt], Point[pt_] :> pt, ∞]]
   {0.420995139971432, 2.380119517396955}

Another method uses an undocumented property of a DataDistribution[] object. In particular, we can extract the "active" domain of such an object like so:

dom = dist["Domain"]
   Interval[{-0.07888027350955408, 1.031712484891653}]

which can then be fed to FindMaximum[]:

FindMaximum[pdf[x], Prepend[First[dom], x] // Evaluate]
   {2.3801667157377957, {x -> 0.4209951363383487}}

Still another method hinges on the fact that the PDF is internally represented as an InterpolatingFunction[], which allows us to extract the interpolation points used by the function. PeakDetect[] can then be used to locate the maximum:

iF = First[Cases[pdf, _InterpolatingFunction, ∞, Heads -> True]];
x0 = First[Pick[First[iF["Coordinates"]], PeakDetect[iF["ValuesOnGrid"]], 1]];
{x0, pdf[x0]}
   {0.42099513633834723, 2.380166715737796}

Find the maximum of a SmoothKernelDistribution

1 Answers1