在 R 中，应用不基于底层分布的参数建模并使用密度的默认核估计器对 10000 个伽马分布变量的方法：

x <- rgamma(10000, 2, 5)
z <- density(x)
plot(z) # always good to check visually
z$x[z$y==max(z$y)]

返回 0.199，这是估计具有最高密度的 x 的值（密度估计存储为“z$y”）。

假设您从大小为 n 的总样本中制作了一个大小为 b 的直方图，并且最大的 bin 有 k 个条目。然后该 bin 内的平均 PDF 可以估计为 b*k/n。

问题在于，另一个具有较少成员总数的 bin 可能具有较高的点密度。只有当您对 PDF 的变化率有合理的假设时，您才能知道这一点。如果这样做，那么您可以估计第二大 bin 实际包含众数的概率。

根本问题是这样的。通过 Kolmogorov-Smirnov 定理，样本可以很好地了解 CDF，因此可以很好地估计中位数和其他分位数。但是知道 L1 中函数的近似值并不能提供其导数的近似知识。因此，在没有额外假设的情况下，没有样本可以很好地了解 PDF。

为了计算 python 中连续分布的模式，我建议三个选项。Scipy 有scipy.stats.mode(data)[0]但它并不准确。对于这三个选项之后我们需要得到一个近似的de PDF函数的数据。一个很好的方法是高斯核密度估计（维基百科）。使用 scipy，我们可以使用distribution = scipy.stats.gaussian_kde(data)并使用distribution.pdf(x)[0]获取数据集的 pdf 值。三种方法搜索分布的最大值，得到域中值的原像。第一次使用 scipy 最小化方法，第二次使用 max() 本机函数，第三次使用 SHGO scipy 方法。这里和比较的​​图像： 以及完整的代码：

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)

最近我遇到了类似的问题，并在Wolfram Mathematica中提出了以下代码：

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

但请记住，这是一个粗略的估计，如果实际的连续分布具有在样本中没有充分表示的窄峰，或者样本恰好包含零星的值簇，则它甚至可能是完全错误的。我相信如果没有关于实际分布的额外信息，就无法量化计算估计中的不确定性。

SmoothKernelDistribution函数有各种选项，您可以尝试调整这些选项以获得针对特定用例的更好结果。