I implemented a Kernel Density Estimator. I have a multivariate dataset that I use with it and I would like to find the point of highest likelihood. A way I thought about is sampling n points using the KDE as my PDF. Then for every point I could use a gradient descent technique and hope to find the maximum.
Is there a better / common way how to do this?
EDIT: This addresses the original question and is retained because it is relevant to that and also to the title of the question, which refers to the peak density. It isn't pertinent to the revised question text, which refers explicitly to the multivariate case.
As the peak of the (estimated) density corresponds to the (main) mode, an alternative is to estimate that mode directly.
A "half-sample mode" can be estimated using recursive selection of the half-sample with the shortest length. Although the idea has longer roots, note particularly the detailed and thorough discussion of Bickel and Frühwirth (2006).
First some more general comments: The mode is often disparaged or neglected by comparison with its siblings the mean and median, but it can be of distinct interest or even use, especially whenever distributions are unimodal but asymmetric. (Modes also have a long history, as readers of Thucydides will recall.)
If a variable is categorical or counted, the mode can usually be read off a frequency table, subject to the occurrence of ties. The same approach can be applied to any variable, subject to the resolution of measurement.
However, the main question is how to get at an estimate of the mode whenever a variable is measured with a resolution such that counting is not a reliable method, if especially all or almost all measurements are distinct. Many people will have been brought up to look at a histogram and read off an approximate value, and may have the impression that not much more can or should be done. Looking at a graph is naturally always a good idea to put any estimate of mode in context. A more modern way of doing it is to get a kernel estimate of the density and modes have been estimated in that way. Either of these approaches suffers from some arbitrariness, for example over bin origin and width or kernel type and width. This shouldn't usually matter, but sometimes a direct method would be useful.
Less obvious than looking for a peak in density, but still worth a try, is to look for a shoulder on a quantile plot.
Kernel estimation is an excellent method, especially when bimodality or multimodality is a possibility. The suggestion, however, is that it may be regarded as an independent method of assessing modality.
An idea of estimating the mode as the midpoint of the shortest interval that contains a fixed number of observations goes back at least to Dalenius (1965). See also Robertson and Cryer (1974), Bickel (2002) and Bickel and Frühwirth (2006) on other estimators of the mode.
The order statistics of a sample of $n$ values of $x$ are defined by
$x_{(1)} \le x_{(2)} \le \cdots \le x_{(n-1)} \le x_{(n)}$.
The half-sample mode can be defined using two rules.
Rule 1 If $n = 1$, the half-sample mode is $x_{(1)}$. If $n = 2$, the half-sample mode is $(x_{(1)} + x_{(2)}) / 2$. If $n = 3$, the half-sample mode is $(x_{(1)} + x_{(2)}) / 2$ if $x_{(1)}$ and $x_{(2)}$ are closer than $x_{(2)}$ and $x_{(3)}$, $(x_{(2)} + x_{(3)}) / 2$ if the opposite is true, and $x_{(2)}$ otherwise.
Rule 2 If $n \ge 4$, we apply recursive selection until left with $3$ or fewer values. First let $h_1 = \lfloor n / 2 \rfloor$. The shortest half of the data from rank $k$ to rank $k + h_1$ is identified to minimise
$x_{(k + h_1)} - x_{(k)}$
over $k = 1, \cdots, n - h_1$. Then the shortest half of those $h_1 + 1$ values is identified using $h_2 = \lfloor h_1 / 2 \rfloor$, and so on. To finish, use Rule 1.
Bickel, D.R. 2002. Robust estimators of the mode and skewness of continuous data. Computational Statistics & Data Analysis 39: 153-163.
Bickel, D.R. and R. Frühwirth. 2006. On a fast, robust estimator of the mode: comparisons to other estimators with applications. Computational Statistics & Data Analysis 50: 3500-3530.
Dalenius, T. 1965. The mode - A neglected statistical parameter.
Journal, Royal Statistical Society A 128: 110-117.
Robertson, T. and J.D. Cryer. 1974. An iterative procedure for estimating the mode. Journal, American Statistical Association 69: 1012-1016.
Note: a Stata implementation is named hsmode (SSC) for which see here. The notes here are based on the help for that program, which includes further discussion and references. Stata users can install using ssc install hsmode. For other implementations see this webpage
A Kernel Density Estimator is defined as $$f(x;h) = \frac{1}{nh}\sum_{i=1}^{n}{K(\frac{x-x_i}{h})}$$ where $K$ is defined as $$\forall x:K(x) > 0, K(-x)=K(x)$$ $$\int_{-\infty}^{\infty}K(x)dx = 1$$ $$\mathbf{E}[K] = 0$$ Now the problem you posed is to find the maximum of the function $f$. The generality of the $K$ function makes the maximum hard to find analytically. Like you mentioned above I think, gradient descent is most likely the way to go. Luckily the expression for the gradient is rather simplistic. $$\frac{\partial f(x;h)}{\partial x} = \frac{1}{nh^2}\sum_{i=1}^{n}{\frac{dK}{dx}}$$ But, the optimization problem is not convex. And therefore gradient descent will not converge to the same point from different starting points. You have a couple options here. One, is what you stated. Sample a set of points and run gradient descent until convergence. And then select the maximum. But this seems computationally expensive, especially in a multidimensional case.
Another approach that is not as computationally expensive and is an approximation to the algorithm you stated above, is to select a point that you believe is close enough to the maximum and then compute gradient descent from there.
From the definition of K above, we can say that the maximum of K (if $K$ is symmetrical around 0) will be at 0. Therefore in function $f$, $x=x_i$. So intuitively a "good enough" starting point should be a point that minimizes the distances from all the data-set points. And this simply will be the mean or expected value of all your samples.
So to reiterate my approach for a computationally inexpensive approximation to the maximum. You first find the expected value of your data-set. Set that as your starting point in your gradient descent algorithm. And you run the optimization until converges.
