如何在 multivariate/3D 中实现核密度估计

Question

我有如下来自at 的数据集，我正在尝试找出具有最佳带宽的内核密度估计。

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

任何人都可以帮助我将其扩展到多变量/在这种情况下是 3D 数据吗？

Answer 1

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

继续使用 scikit-learn
使用不同的库。例如，如果您感兴趣的内核是高斯 - 那么您可以使用 scipy.gaussian_kde which is arguably easier to understand / apply. There is a very good example of this technique in this question.
从第一原则开始你自己的。这很难，我不推荐

This blog post 详细介绍了内核密度估计 (KDE) 的各种库实现的相对优点。

我将向您展示最简单的方法（在我看来 - 是的，这有点基于意见），我认为这是您的情况的选项 2。

注意此方法使用链接文档中描述的经验法则来确定带宽。使用的确切规则是 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...

警告

这可能会产生糟糕的结果，具体取决于您的具体问题。要记住的事情显然是：

随着维度数量的增加，您想要的观察数据点数量将呈指数级增长——3 维中 9 个点的样本数据非常稀疏——尽管我假设点表示实际上你还有很多。
如正文中所述 - 以特定方式选择带宽 - 这可能导致估计的 pdf 过度（或可以想象但不太可能不足）平滑。

如何在 multivariate/3D 中实现核密度估计

How to implement Kernel density estimation in multivariate/3D

python

numpy

machine-learning

kernel-density

scikit-learn