我有一个如下格式的数据集,并且我想找出最佳带宽的核密度估计。
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]])
但我不知道该如何着手解决这个问题。另外,如何找到Σ矩阵呢?
更新
我尝试使用scikit-learn工具包中的KDE函数来找出一元(1D)kde,
# kde 函数def kde_sklearn(x, x_grid, bandwidth): kde = KernelDensity(kernel='gaussian', bandwidth=bandwidth).fit(x) log_pdf = kde.score_samples(x_grid[:, np.newaxis]) return np.exp(log_pdf)# 最佳带宽选择from sklearn.grid_search import GridSearchCVgrid = GridSearchCV(KernelDensity(), {'bandwidth': np.linspace(.1, 1.0, 30)}, cv=20)grid.fit(x)bw = grid.best_params_# 使用 kde 计算 pdfpdf = kde_sklearn(x, x_grid, bw)ax.plot(x_grid, pdf, label='bw={}'.format(bw))ax.legend(loc='best')plt.show()
有谁能帮我将其扩展到多变量/在这种情况下是3D数据吗?
回答:
有趣的问题。你有几种选择:
- 继续使用scikit-learn
- 使用不同的库。例如,如果你感兴趣的核是高斯核,那么你可以使用scipy.gaussian_kde,它可能更容易理解和应用。这项技术的一个很好的例子可以在这个问题中找到。
- 从基本原理开始自己实现。这非常困难,我不推荐这样做
这篇博客文章详细介绍了各种库实现核密度估计(KDE)的相对优点。
我将向你展示(在我看来——是的,这有点基于意见)最简单的方法,我认为在你的情况下是选项2。
注意 这种方法使用链接文档中描述的经验法则来确定带宽。使用的具体规则是Scott的规则。你提到Σ矩阵让我认为经验法则带宽选择对你来说是可以的,但你也谈到了最佳带宽,并且你展示的例子使用交叉验证来确定最佳带宽。因此,如果这种方法不适合你的目的,请在评论中告诉我。
import numpy as npfrom scipy import statsdata = 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需要N个长度为K的向量来表示K个数据点 #而不是K个长度为N的向量kde = stats.gaussian_kde(data)# 现在你有了你的kde!! 解释它/可视化它在3D数据中可能很困难# 你可能想先尝试2D数据——然后你可以将结果的估计pdf作为第三维度的高度进行绘制,使可视化更容易。# 这是在一个规则的n维网格上评估估计pdf的基本方法# 创建一个规则的N维网格,每个维度有(任意)20个点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)# 将网格转换为每个点的N维坐标# 注意 - 随着N的增加,coords会变得非常大...coords = np.vstack(list(map(np.ravel, grid)))# 在每个坐标处评估KD估计的pdfdensity = kde(coords)# 你可以在这里对密度值做任何你想做的事情..# 绘制它们,输出它们,在其他地方使用它们...
警告
这可能会根据你的具体问题产生糟糕的结果。需要注意的显然是:
-
随着维度的增加,你想要的观测数据点的数量将必须以指数级增加——你在3维中的9个数据点样本相当稀疏——尽管我假设这些点实际上表示你有更多的数据点。
-
正如正文中提到的——带宽是以特定方式选择的——这可能会导致估计的pdf过度(或不太可能但可能不足)平滑。