16
$\begingroup$

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.

$\endgroup$
17
  • 8
    $\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$
    – Macro
    Dec 18, 2011 at 0:50
  • 1
    $\begingroup$ @Macro - that's helpful. You can then read my question as, 'What are the best methods to determine the (peak) density?' $\endgroup$
    – keflavich
    Dec 18, 2011 at 0:52
  • 1
    $\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$
    – Macro
    Dec 18, 2011 at 0:53
  • 2
    $\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$ Dec 18, 2011 at 1:53
  • 3
    $\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$
    – whuber
    Dec 18, 2011 at 16:53

5 Answers 5

5
$\begingroup$

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

$\endgroup$
2
  • 5
    $\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$ Feb 17, 2012 at 12:02
  • $\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$
    – keflavich
    Oct 5, 2021 at 13:36
5
$\begingroup$

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.

$\endgroup$
2
  • $\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$
    – wnoise
    Feb 9, 2022 at 18:44
  • 1
    $\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$ Feb 11, 2022 at 12:27
2
$\begingroup$

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: enter image description here And the complete code:

import numpy as np
import scipy.stats
import scipy.optimize
import matplotlib as mpl
import matplotlib.pyplot as plt
import random

mpl.style.use("ggplot")

def plot_histogram(data, distribution, modes):
    plt.figure(figsize=(8, 4))
    plt.hist(data, density=True, ec='white')
    plt.title('HISTOGRAM')
    plt.xlabel('Values')
    plt.ylabel('Frequencies')
    
    
    x_plot = np.linspace(min(data), max(data), 1000)
    y_plot = distribution.pdf(x_plot)
    plt.plot(x_plot, y_plot, linewidth=4, label="PDF KDE")
    
    for name, mode in modes.items():
        plt.axvline(mode, linewidth=2, label=name+": "+str(mode)[:7], color=(random.uniform(0, 1), random.uniform(0, 1), random.uniform(0, 1)))
    
    plt.legend(title='DISTRIBUTIONS', bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.show()
    

## SCIPY MODE
def calc_scipy_mode(data):
    return scipy.stats.mode(data)[0]

## METHOD 1: MAXIMIZE PDF SCIPY MINIMIZE
def calc_minimize_mode(data, distribution):
    def objective(x):
        return 1/distribution.pdf(x)[0]
    
    bnds = [(min(data), max(data))]
    solution = scipy.optimize.minimize(objective, [1], bounds = bnds)
    
    return solution.x[0]

## METHOD 2: MAXIMIZE PDF AND GET PREIMAGE
def calc_max_pdf_mode(data, distribution):
    x_domain = np.linspace(min(data), max(data), 1000)
    y_pdf = distribution.pdf(x_domain)
    i = np.argmax(y_pdf)
    return x_domain[i]

## METHOD 3: ## METHOD 3: MAXIMIZE PDF SCIPY SHGO
def calc_shgo_mode(data, distribution):
    def objective(x):
        return 1/distribution.pdf(x)[0]
    
    bnds = [[min(data), max(data)]]
    solution = scipy.optimize.shgo(objective, bounds= bnds, n=100*len(data))
    return solution.x[0]



def calculate_mode(data):
    ## KDE
    distribution = scipy.stats.gaussian_kde(data)
    
    scipy_mode = calc_scipy_mode(data)[0]
    minimize_mode = calc_minimize_mode(data, distribution)
    max_pdf_mode = calc_max_pdf_mode(data, distribution)
    shgo_mode = calc_shgo_mode(data, distribution)
    
    modes = {
        "scipy_mode": scipy_mode, 
        "minimize_mode": minimize_mode, 
        "max_pdf_mode": max_pdf_mode,
        "shgo_mode": shgo_mode
    }
    plot_histogram(data, distribution, modes)
    
    
if __name__ == "__main__":
  
    data = [d1, d2, d3, ..........]
    
    calculate_mode(data)
$\endgroup$
3
  • 1
    $\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$ Aug 8, 2021 at 0:55
  • $\begingroup$ what if there are 2 modes? $\endgroup$
    – Adrian
    Jun 7, 2022 at 6:01
  • $\begingroup$ Find the preimage of the greatest pdf and then the preimage of second greatest pdf $\endgroup$ Jun 8, 2022 at 3:56
1
$\begingroup$

Recently I faced a similar problem, and came up with this code in Wolfram Mathematica:

ModeEstimate[data_?VectorQ] := 
  MaximalBy[data, PDF[SmoothKernelDistribution[data]], 1][[1]];

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.

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.