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

jqu*_*404 5 python numpy machine-learning kernel-density scikit-learn

我有像下面这样的数据集,我试图找出具有最佳带宽的核密度估计。

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]])
Run Code Online (Sandbox Code Playgroud)

但我不知道如何处理它。还怎么找的矩阵。

更新

我尝试了 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()
Run Code Online (Sandbox Code Playgroud)

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

J R*_*ape 5

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

\n\n
    \n
  1. 继续 scikit-learn
  2. \n
  3. 使用不同的库。例如,如果您感兴趣的内核是高斯 - 那么您可以使用scipy.gaussian_kde,这可以说更容易理解/应用。这个问题中有一个很好的例子来说明这种技术
  4. \n
  5. 从第一原则开始你自己的。这非常困难,我不推荐它
  6. \n
\n\n

这篇博文详细介绍了内核密度估计 (KDE) 的各种库实现的相对优点。

\n\n
\n\n

我将向您展示什么(在我看来 - 是的,这有点基于意见)是最简单的方法,我认为在您的情况下是选项 2。

\n\n

注意此方法使用链接文档中所述的经验法则来确定带宽。使用的确切规则是斯科特规则。您提到的 \xce\xa3 矩阵让我认为经验法则带宽选择对您来说是可以的,但您也谈论了最佳带宽,并且您提供的示例使用交叉验证来确定最佳带宽。因此,如果此方法不适合您的目的 - 请在评论中告诉我。

\n\n
import numpy as np\nfrom scipy import stats\ndata = np.array([[1, 4, 3], [2, .6, 1.2], [2, 1, 1.2],\n         [2, 0.5, 1.4], [5, .5, 0], [0, 0, 0],\n         [1, 4, 3], [5, .5, 0], [2, .5, 1.2]])\n\ndata = data.T #The KDE takes N vectors of length K for K data points\n              #rather than K vectors of length N\n\nkde = stats.gaussian_kde(data)\n\n# You now have your kde!!  Interpreting it / visualising it can be difficult with 3D data\n# You might like to try 2D data first - then you can plot the resulting estimated pdf\n# as the height in the third dimension, making visualisation easier.\n\n# Here is the basic way to evaluate the estimated pdf on a regular n-dimensional mesh\n# Create a regular N-dimensional grid with (arbitrary) 20 points in each dimension\nminima = data.T.min(axis=0)\nmaxima = data.T.max(axis=0)\nspace = [np.linspace(mini,maxi,20) for mini, maxi in zip(minima,maxima)]\ngrid = np.meshgrid(*space)\n\n#Turn the grid into N-dimensional coordinates for each point\n#Note - coords will get very large as N increases...\ncoords = np.vstack(map(np.ravel, grid))\n\n#Evaluate the KD estimated pdf at each coordinate\ndensity = kde(coords)\n\n#Do what you like with the density values here..\n#plot them, output them, use them elsewhere...\n
Run Code Online (Sandbox Code Playgroud)\n\n
\n\n

警告

\n\n

这可能会产生可怕的结果,具体取决于您的具体问题。显然要记住的事情是:

\n\n
    \n
  1. 随着维度数量的增加,您想要的观察数据点的数量将不得不呈指数级增长 - 您的 3 维 9 个点的样本数据非常稀疏 - 尽管我认为这些点表明实际上您还有更多。

  2. \n
  3. 正如正文中提到的 - 以特定方式选择带宽 - 这可能会导致估计的 pdf 过度(或可以想象但不太可能不足)平滑。

  4. \n
\n

  • 您还需要查看“kde.covariance”。`kde.factor` 乘以 `kde.covariance` 以获得我认为您想要看到的内核协方差矩阵或带宽(我认为您在上面称之为 Σ )。[文档页面](http://docs.scipy.org/doc/scipy-0.15.1/reference/ generated/scipy.stats.gaussian_kde.html)底部有详细说明 (2认同)