如何在核密度估计中找到局部极大值？

Question

我正在尝试使用核密度估计器(KDE)制作一个过滤器(用于去除离群值和噪声)。我在我的三维(d=3)数据点中应用了KDE，并得到了概率密度函数(PDF) f(x)。现在，正如我们所知道的那样，密度估计f(x)的局部极大值定义了数据点簇的中心。因此，我的想法是定义适当的f(x)，以确定这些簇。

我的问题是如何以及什么方法最适合找到f(x)中的局部极大值。如果有人能够提供一些示例代码/思路，我将不胜感激。

以下是查找3D数据中给出f(x)的KDE的代码。

import numpy as np
from scipy import stats

data = np.array([[1, 4, 3], [2, .6, 1.2], [2, 1, 1.2],
         [2, 0.5, 1.4], [5, .5, 0], [0, 0, 0],
         [1, 4, 3], [5, .5, 0], [2, .5, 1.2]])
data = data.T 
kde = stats.gaussian_kde(data)
minima = data.T.min(axis=0)
maxima = data.T.max(axis=0)
space = [np.linspace(mini,maxi,20) for mini, maxi in zip(minima,maxima)]
grid = np.meshgrid(*space)
coords = np.vstack(map(np.ravel, grid))
#Evaluate the KD estimated pdf at each coordinate
density = kde(coords)

- jquery404

Answer 1

这是一个简短的函数，演示了如何估计最大值。注意：采样数量no_samples越高，最大值的估计越准确。

from scipy.stats import gaussian_kde
import numpy as np

 def estimate_maxima(data):
    kde = gaussian_kde(data)
    no_samples = 10
    samples = np.linspace(min(data), max(data), no_samples)
    probs = kde.evaluate(samples)
    maxima_index = probs.argmax()
    maxima = samples[maxima_index]
    
    return maxima

- brandon musa

---

样例假设数据分布在0-10之间，应该是： samples = np.linspace(min(data), max(data), no_samples) - darc

Answer 2

你需要使用一种名为Mean Shift的算法。它是一种聚类算法，通过找到KDE的模式（也称f(x)的最大值）来工作。请注意，为KDE设置的带宽将影响模式的数量和位置。由于你正在使用Python，scikit-learn中有一个实现。

- Raff.Edward

---

谢谢你的建议。我按照你的建议，对我的密度值应用了meanshift算法。但是我不确定如何获取局部极大值。它给了我6个聚类:(。这是源代码，我做得对吗？ - jquery404

聚类中心“应该”包含极大值，因为“中心”并没有太多意义，因为聚类形状可以非常不规则。 - Raff.Edward

在中等规模的数据（1e4）上非常慢。 - lovetl2002

是的，Python相对较慢 - 默认的scikit实现也是如此。它有一个简单的分箱选项可以帮助提速，并且还有其他更快的均值漂移方法。 - Raff.Edward

Answer 3

你可以使用scipy.optimize。

1D数据示例：

import numpy as np
from scipy import optimize
from scipy import stats


# Generate some random data
shape, loc, scale = .5, 3, 10
n = 1000
data = np.sort(stats.lognorm.rvs(shape, loc, scale, size=n))

kernel = stats.gaussian_kde(data)
# Minimize the negative instead of maximizing
# Depending on the shape of your data, you might want to set some bounds
opt = optimize.minimize_scalar(lambda x: -kernel(x))
opt

     fun: array([-0.08363781])
    nfev: 21
     nit: 14
 success: True
       x: array([10.77361776])

这个分布的实际模式位于

mode = scale/np.exp(shape**2) + loc
mode
10.788007830714049

绘制结果：

import matplotlib.pyplot as plt

data_es = np.linspace(0, data.max(), 201)  # x-axis points
ecdf = (np.arange(n) + 1)/n  # empirical CDF

fig, axes = plt.subplots(2, 1, sharex=True, dpi=300, figsize=(6,7))
axes[0].hist(x, bins=30, density=True, alpha=.5, rwidth=.9)  # histogram
axes[0].plot(data_es, kernel.pdf(data_es), 'C0')  # estimated PDF
axes[0].plot(data_es, stats.lognorm.pdf(data_es, shape, loc, scale), 'k--', alpha=.5)  # true PDF
axes[0].plot(opt.x, kernel.pdf(opt.x), 'C0.')  # estimated mode
axes[0].plot(mode, stats.lognorm.pdf(mode, shape, loc, scale), 'k.', alpha=.5)  # true mode

axes[1].plot(np.sort(data), ecdf)  # estimated CDF
axes[1].plot(opt.x, np.interp(opt.x, np.sort(data), ecdf), 'C0.')  #estimated mode
axes[1].plot(data_es, stats.lognorm.cdf(data_es, shape, loc, scale), 'k--', alpha=.5)  # true CDF
axes[1].plot(mode, stats.lognorm.cdf(mode, shape, loc, scale), 'k.', alpha=.5)  # true mode

fig.tight_layout()

从图中可以看出，估计的模式非常贴切。我认为可以使用scipy.optimize中的其他方法将其扩展到多变量数据。

- Aubergine