如何从内核密度估计中获取内核（最好是 sklearn.neighbors）？

Posted 2023-03-12

【中文标题】如何从内核密度估计中获取内核（最好是 sklearn.neighbors）？

【英文标题】：How to get Kernels from kernel density estimation (preferrably sklearn.neighbors)?

【发布时间】：2017-12-31 10:08:09

【问题描述】：

我目前正在对时间序列数据集进行一些季节性估计。

我得到的是数据集中可能出现的频率/周期的数据集。因此，这有点嘈杂（例如，有一些周期为 [100, 98, 101, 102]，实际上应该是“相同的”）。

为了估计尖锐周期，我尝试通过核密度估计（kde、sklearn.neighbors.KernelDensity）来估计峰值，如下所示：

import numpy as np
from sklearn.neighbors import KernelDensity
from scipy import signal
import matplotlib.pyplot as plt

X1 = np.random.randint(1, 4, 20)
X2 = np.random.randint(10, 13, 200)
X = np.concatenate((X1, X2), axis=0)
# the peaks schould be at 2 and 11!

bw = 1

kde = KernelDensity(kernel='gaussian', bandwidth=bw).fit(X.reshape(-1, 1))
estimator = np.linspace(0, 15, 100)
kde_est = np.exp(kde.score_samples(estimator.reshape(-1, 1)))

plt.plot(estimator, kde_est)

peaks_pos = signal.argrelextrema(kde_est, np.greater)[0]

print(estimator[peaks_pos])
# the peaks are at around 2 and 11!

另外，我想知道这个估计的内核是什么样子的。对于高斯情况，应该有一组 /mu 和 /sigma 应该可用于所有 [默认] 40 个内核。我可以访问这些信息吗？ 我在文档或 kde 属性的详细信息中找不到线索。但我很确定，这应该可以在这里找到。

为了澄清，我为什么需要这个：

在以下示例中，两个峰靠得太近而无法找到，但我确信内核会显示出来。

X1 = np.random.randint(1, 4, 20)
X2 = np.random.randint(5, 8, 200)
X = np.concatenate((X1, X2), axis=0)
# the peaks schould be at 2 and 6!

bw = 1

kde = KernelDensity(kernel='gaussian', bandwidth=bw).fit(X.reshape(-1, 1))
estimator = np.linspace(0, 15, 100)
kde_est = np.exp(kde.score_samples(estimator.reshape(-1, 1)))

plt.plot(estimator, kde_est)

peaks_pos = signal.argrelextrema(kde_est, np.greater)[0]

print(estimator[peaks_pos])
# the peaks are at around 6 and sometimes 2!

【参考方案1】：

我相信在核密度估计中找不到您要查找的内容。 KDE 中的所有内核都具有完全相同的形状（标准差）并且以数据点为中心（因此均值由 X 中的值决定）。

为了防止正态分布与模糊峰值的接近，您可以做的是调整带宽（如果您的第二个样本，我设法通过使用 0.7 的带宽来获得非常一致的 2 个峰值。有代数方法可以做到这一点（参见：***），或者您可以使用交叉验证为您的样本选择最佳带宽（参见：blog on the subject)。

但是，如果您想将数据集拆分为由具有各种形状（权重、均值和协方差）的正态分布描述的不同分量，您可能需要使用高斯混合建模。您可以在下面找到一个示例。为了确定组件的最佳数量，有多种方法，例如轮廓标准或 akaike 信息标准（内置于 scikitlearn）。由于我们知道示例中有 2 个正态分布，因此我没有实施这样的标准，但您可以在互联网上轻松找到更多信息。

X1 = np.random.randint(1, 4, 20)
X2 = np.random.randint(5, 8, 200)
X = np.concatenate((X1, X2), axis=0)
# the peaks schould be at 2 and 6!

components = 2

gmm = GaussianMixture(n_components = components).fit(X.reshape(-1,1))

#you can now directly get the means from the gaussian mixture models components,
#skipping the score_samples and signal.argrelextrema steps.
print gmm.means_
#the means are around 2 and 6!


#your original method of getting the peaks:
estimator = np.linspace(0, 15, 100)
gmm_est = np.exp(gmm.score_samples(estimator.reshape(-1,1)))

plt.hist(X,normed=True)
plt.plot(estimator,gmm_est,linewidth=5,color='black',alpha=0.7)


peaks_pos = signal.argrelextrema(gmm_est, np.greater)[0]

print(estimator[peaks_pos])


#plotting the separate components:
for n,weight in enumerate(gmm.weights_):
    plt.plot(estimator,weight*stats.norm.pdf(estimator,gmm.means_[n][0],np.sqrt(gmm.covariances_[n][0][0])))
plt.show()

image of results

【讨论】：

感谢您的精彩评论。我将尝试您提出的高斯去混合方案。对于带宽：我已经有了一个粗略的估计。交叉验证不是一个选项，因为这会花费很多时间（=> 应该能够在 3D 堆栈上快速运行）。

比交叉验证更快的带宽估计也可以在 statsmodels 包中找到（Silverman 或 Scott 的经验法则）仅供参考。 statsmodels.org/stable/generated/…

是的，这就是我实现的（不使用这个包）:)

我在我的“更复杂”的系统上尝试了他的，但它不够稳定。谢谢你的建议。 scikitlearn.mixture.GaussianMixture 符合我的要求，但我的信号太弱了。