时间戳

首页 > 解决方案 > 如何使用 sklearn 加快核密度估计？

python - 如何使用 sklearn 加快核密度估计？

问题描述

我正在尝试使用sklearn.neighbors.KernelDensity估计概率密度函数（PDF） 。但是，我不知道用于带宽的最佳值。因此，我使用sklearn.model_selection.GridSearchCV来计算最佳带宽（我从阅读这篇文章中得到了想法）。

我的代码真的很慢，特别是如果我将数据点的数量增加到 1000 以上。我希望数据点的数量为 10^6。你知道我怎样才能加快代码速度吗？或者你知道我可以估计最佳带宽的另一种方法吗？

我附上了一个示例代码。这里的 PDF 是高斯的。但是，实际上，我不知道 PDF，需要计算一个估计值。

import matplotlib.pyplot as plt
import numpy as np
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV, LeaveOneOut

mu = 0
sigma = 1
num_markers = int(1e3)
marker_positions = np.random.normal(mu, sigma, num_markers)

x_min = -3
x_max = 3
nx = 101
x = np.linspace(x_min, x_max, nx)

bandwidths = 10 ** np.linspace(-1, 0, 5)
grid = GridSearchCV(KernelDensity(kernel='gaussian'), \
                    {'bandwidth': bandwidths}, \
                    cv=LeaveOneOut(), \
                    n_jobs = -1, \
                    verbose = 1)
grid.fit(marker_positions[:, None])

kde = grid.best_estimator_.fit(marker_positions[:, None])
log_density = kde.score_samples(x[:, None])

fig, ax = plt.subplots()
ax.plot(x, np.exp(-x**2/2) / np.sqrt(2 * np.pi), label = 'Exact')
ax.plot(x, np.exp(log_density), label = 'Estimate')
ax.set_ylabel('Probability density')
ax.set_xlabel('x')
ax.legend()
plt.savefig('probability_density_estimate.png')

这是输出：

Fitting 1000 folds for each of 5 candidates, totalling 5000 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 16 concurrent workers.
[Parallel(n_jobs=-1)]: Done  18 tasks      | elapsed:    0.9s
[Parallel(n_jobs=-1)]: Done 257 tasks      | elapsed:    2.2s
[Parallel(n_jobs=-1)]: Done 757 tasks      | elapsed:    4.8s
[Parallel(n_jobs=-1)]: Done 1457 tasks      | elapsed:    8.7s
[Parallel(n_jobs=-1)]: Done 2357 tasks      | elapsed:   14.6s
[Parallel(n_jobs=-1)]: Done 3457 tasks      | elapsed:   21.4s
[Parallel(n_jobs=-1)]: Done 4757 tasks      | elapsed:   29.4s
[Parallel(n_jobs=-1)]: Done 4969 out of 5000 | elapsed:   30.7s remaining:    0.2s
[Parallel(n_jobs=-1)]: Done 5000 out of 5000 | elapsed:   30.8s finished

和图表： 在此处输入图像描述

标签： pythonnumpyscikit-learncross-validation

解决方案

有多种方法可以选择带宽。交叉验证是其中之一，但对于更快的方法，您可以参考一些经验法则。SciPy 的高斯核密度估计实现 ( scipy.stats.gaussian_kde) 内置了两个这样的经验法则：

1. 斯科特法则： n**(-1. / (d + 4))

2. 西尔弗曼法则： (n * (d + 2) / 4.)**(-1. / (d + 4)).

在这两者中，d是数据的维度，n是您拥有的数据点的数量。对于加权数据集，您可以将n这些方程替换为：

neff = sum(weights)^2 / sum(weights^2)

所有这些都直接来自链接页面；那里有更多细节，以及参考链接。希望这些经验法则之一对您有用，或者至少链接页面可能能够引导您走上一条富有成效的道路！

推荐阅读