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?

Answer 1

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

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}