Computing the mode of data sampled from a continuous distribution 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.
 A: 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").
A: 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.
A: 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:
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)

A: 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.
A: 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.
