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

Posted 2023-02-14

【中文标题】如何在核密度估计中找到局部最大值？

【英文标题】：How to find Local maxima in Kernel Density Estimation?

【发布时间】：2015-09-20 18:38:30

【问题描述】：

我正在尝试使用内核密度估计器 (KDE) 制作一个过滤器（以去除异常值和噪声）。我在我的 3D (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)

【参考方案1】：

您将需要使用称为Mean Shift 的算法。它是一种通过查找 KDE 的模式（又名 f(x) 的最大值）来工作的聚类算法。请注意，为您的 KDE 设置的带宽会影响模式的数量及其位置。由于您使用的是 python，因此scikit-learn 中有一个实现。

【讨论】：

感谢您的想法。我听从了您的建议，并将均值偏移应用于我的密度值。但我不确定如何获得局部最大值。它给了我 6 个集群 :( 。这是Source Code，我做得对吗？

簇中心应该包含最大值，因为“中心”没有多大意义，因为簇形状可能非常不规则。

【参考方案2】：

这是一个简短的函数，演示了如何估计最大值。注意：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(0, 10, no_samples)

      probs = kde.evaluate(samples)

      maxima_index = probs.argmax()

      maxima = samples[maxima_index]

      return maxima

【参考方案3】：

你可以使用 scipy.optimize。

一维数据示例：

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 中的其他方法将其扩展到多变量数据。