如何在多元/三维中实现核密度估计

Question

如何在多元/三维中实现核密度估计

pythonnumpymachine-learningscikit-learnkernel-density

5

我有以下格式的数据集，正在尝试找到最佳带宽的核密度估计。

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 function
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)

# optimal bandwidth selection
from sklearn.grid_search import GridSearchCV
grid = GridSearchCV(KernelDensity(), {'bandwidth': np.linspace(.1, 1.0, 30)}, cv=20)
grid.fit(x)
bw = grid.best_params_

# pdf using kde
pdf = kde_sklearn(x, x_grid, bw)
ax.plot(x_grid, pdf, label='bw={}'.format(bw))
ax.legend(loc='best')
plt.show()

有人能帮我将这个扩展到多变量/在这种情况下是三维数据吗？

- jquery404

我想知道是否可以帮忙，但需要更多了解。我能看到每个数据点将有三个值，但是您编写的方式是将这些三元组进一步分组为三个一组。您的输入数据为什么要分两次组？另外，我想再次确认您所说的Σ矩阵的含义。我假设您指的是估计的数据协方差-因此可以使用Σ^(-1/2)的经验法则带宽？如果是这样，您是打算将其作为带宽优化的起点，还是代替优化？ - J Richard Snape

我的回答有帮到您吗？如果没有，欢迎留下评论，因为我可能能够根据您的需求进行调整。 - J Richard Snape

@JRichardSnape 你说得对，我把数据分组错了，实际上在我的代码中它就像你实现的那样，但是当我把代码复制到SO时搞砸了。是的，我所说的Σ是协方差矩阵。 - jquery404

好的。但我仍然不确定我的下面的回答是否有帮助 - 这是否满足了您的需求？还是您的问题还有其他方面需要解决？如果您想输出协方差矩阵，我可以添加它。 - J Richard Snape

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- J Richard Snape · Accepted Answer

有趣的问题。你有几个选择：

继续使用scikit-learn
使用不同的库。例如，如果你感兴趣的内核是高斯内核，则可以使用scipy.gaussian_kde，这可能更容易理解/应用。这里有一个非常好的例子在这个问题中。
从头开始自己编写。这非常困难，我不建议这样做。

这篇博客文章详细介绍了各种库实现核密度估计（KDE）的相对优点。

我将向您展示我认为是最简单的方法（在我看来，这有点基于个人意见），我认为这是您情况下的第二种选择。请注意，此方法使用拇指规则来确定带宽，如链接文档中所述。使用的确切规则是Scott的规则。您提到的Σ矩阵使我想到拇指规则带宽选择对您来说是可以接受的，但您也谈到了最佳带宽，并且您提供的示例使用交叉验证来确定最佳带宽。因此，如果此方法不适合您的目的，请在评论中告诉我。

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 #The KDE takes N vectors of length K for K data points
              #rather than K vectors of length N

kde = stats.gaussian_kde(data)

# You now have your kde!!  Interpreting it / visualising it can be difficult with 3D data
# You might like to try 2D data first - then you can plot the resulting estimated pdf
# as the height in the third dimension, making visualisation easier.

# Here is the basic way to evaluate the estimated pdf on a regular n-dimensional mesh
# Create a regular N-dimensional grid with (arbitrary) 20 points in each dimension
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)

#Turn the grid into N-dimensional coordinates for each point
#Note - coords will get very large as N increases...
coords = np.vstack(map(np.ravel, grid))

#Evaluate the KD estimated pdf at each coordinate
density = kde(coords)

#Do what you like with the density values here..
#plot them, output them, use them elsewhere...

注意

这可能会导致糟糕的结果，具体取决于您的问题。需要考虑以下几点：

随着维度的增加，您想要观察到的数据点数量将呈指数级增长 - 您在三个维度中拥有9个样本数据点相对较少 - 虽然我假设点实际上更多。
如主体部分所述 - 带宽是以特定方式选择的 - 这可能导致估计的pdf过度平滑（或者极少情况下会低估平滑）。