Stack Exchange network consists of 183 Q&A communities including Stack Overflow, the largest, most trusted online community for developers to learn, share their knowledge, and build their careers.
What are the best methods for fitting the 'mode' of data sampled from a continuous distribution?
Since the mode is technically undefined (right?) for a continuous distribution, I'm really asking 'how do you find the most common value'?
If you assume the parent distribution is gaussian, you could bin the data and find say the mode is the bin location with the greatest counts. However, how do you determine the bin size? Are there robust implementations available? (i.e., robust to outliers). I use python/scipy/numpy, but I can probably translate R without too much difficulty.
$\begingroup$I'm not sure if the mode is technically defined this way, but the global mode of a continuous distribution is usually taken to mean the point with the highest density.$\endgroup$
$\begingroup$Maybe fit a kernel density estimate for your data and estimate the mode as the peak of that? This seems like a reasonable approach but I'm not familiar with the literature on this problem.$\endgroup$
$\begingroup$If you don't assume the parent distribution is gaussian, is it still possible to bin the data and take the mode to be the bin location with the largest count? Why or why not? On a more serious note, why not find the deciles $x_0=x_{\min},x_1,x_2,\ldots,x_9,x_{10}=x_{\max}$ so that $10\%$ of the samples are in the interval $x_{i+1}-x_i$, and so it is likely that the mode is in the shortest interdecile interval $\min_{1 \leq j \leq 10} x_{j+1}-x_j$? Then take the bin size to be, say, one-fourth of this shortest interdecile interval.$\endgroup$
$\begingroup$What assumptions can you make about the parent distribution, keflavich? If they are parametric, it's best to estimate the parameters and then estimate the mode from those parameters. (E.g., the sample mean estimates the mode of a normal distribution.) If not, binning can be a poor method. Instead, a sequence of kernel estimators with varying halfwidth can be used to provide a sequence of estimators; typically, if the underlying distribution is unimodal, the modes of the kernel smooths will appear to converge towards a unique mode as the halfwidths get large and that can be your estimate.$\endgroup$
In R, applying the method that isn't based on parametric modelling of the underlying distribution and uses the default kernel estimator of density to 10000 gamma distributed variables:
x <- rgamma(10000, 2, 5)
z <- density(x)
plot(z) # always good to check visually
z$x[z$y==max(z$y)]
returns 0.199 which is the value of x estimated to have the highest density (the density estimates are stored as "z$y").
$\begingroup$The only thing I would do differently from that is use a different bandwidth. The default bandwidth for density() is not particularly good. density(x,bw="SJ") is better. Even better would be to use a bandwidth designed for mode estimation. See sciencedirect.com/science/article/pii/0167715295000240 for some discussion.$\endgroup$
$\begingroup$This is the right answer to the question I had asked ~10 years ago - basically, find the peak point of a kernel density estimate.$\endgroup$
Suppose you make a histogram, of bin size b, and the largest bin has k entries, from your total sample of size n. Then the average PDF within that bin can be estimated as b*k/n.
The problem is that another bin, which has fewer total members, could have a high spot density. You can only know about this if you have a reasonable assumption about the rate of change of the PDF. If you do, then you can estimate the probability that the second largest bin actually contains the mode.
The underlying problem is this. A sample provides good knowledge of the CDF, by the Kolmogorov-Smirnov theorem, and so a good estimate of median and other quantiles. But knowing an approximation to a function in L1 does not provide approximate knowledge of its derivative. So no sample provides good knowledge of the PDF, without additional assumptions.
$\begingroup$Surely in practice it provides some approximate knowledge of the derivative. In an interval, the change in the CDF does directly give you the mean of the PDF in that interval, and if it is merely continuous, then that mean value must even be obtained. You might need some Lipschitz or similar condition to provide good bounds, but the samples closer to other samples are certainly more deserving to be called modes. I think in practice the best you can do is adaptive Kernel Density Estimation, with smaller bandwidth at areas more highly sampled.$\endgroup$
$\begingroup$KDE involves assumptions about smoothness. These assumptions imply information about the derivative - but the data do not. Here's a case where this would go wrong: 99% of the data are normally distributed with mean 0.0 and SD 1.0, the other 1% are all 0.0. It is a shame to rely on KDE and the associated assumptions, when most of the questions that actually matter can be answered from the CDF.$\endgroup$
In order to calculate the mode of a continious distribution in python I suggest three options. Scipy have scipy.stats.mode(data)[0] but it's not exact. For the three options after we need to get an approximation of de PDF function of the data. An excellent method is Gaussian Kernel
Density Estimation(WIKIPEDIA). With scipy we can use distribution = scipy.stats.gaussian_kde(data) and get pdf value of the dataset with distribution.pdf(x)[0]. The three methods search the max value of the distribution and get the preimage of the value in the the domian. The first use scipy minimize method, the second use max() native function, and the third SHGO scipy method. Here and imgage of comparison:
And the complete code:
$\begingroup$Although implementation is often mixed with substantive content in questions, we are supposed to be a site for providing information about statistics, machine learning, etc., not code. It can be good to provide code as well, but please elaborate your substantive answer in text for people who don't read this language well enough to recognize & extract the answer from the code.$\endgroup$
But keep in mind it is a rough estimate, and can even be totally wrong if the actual continuous distribution has narrow peaks that are not adequately represented in the sample, or the sample happens to contain sporadic clusters of values. I believe there is no way to quantify uncertainty in the computed estimate without additional information about the actual distribution.
SmoothKernelDistribution function has various options that you can try to adjust to get better results for your specific use cases.
Here are some general solution sketches that also work for high-dimensional distributions:
Train an f-GAN with reverse KL divergence, without giving any random input to the generator (i.e. force it to be deterministic).
Train an f-GAN with reverse KL divergence, move the input distribution to the generator towards a Dirac delta function as training progresses, and add a gradient penalty to the generator loss function.
Train a (differentiable) generative model that can tractably evaluate an approximation of the pdf at any point (I believe that e.g. a VAE, a flow-based model, or an autoregressive model would do). Then use some type of optimization (some flavor of gradient ascent can be used if model inference is differentiable) to find a maximum of that approximation.