scipy - 将 scipy.stats.gaussian_kde 与二维数据一起使用

Question

我正在尝试使用 the scipy.stats.gaussian_kde class将一些用经纬度信息收集的离散数据进行平滑处理，使其最终显示为有点类似于等高线图，其中高密度为峰，低密度为谷。

我很难将二维数据集放入 gaussian_kde类(class)。我一直在尝试弄清楚它是如何处理一维数据的，所以我认为二维应该是这样的:

from scipy import stats
from numpy import array
data = array([[1.1, 1.1],
              [1.2, 1.2],
              [1.3, 1.3]])
kde = stats.gaussian_kde(data)
kde.evaluate([1,2,3],[1,2,3])

这就是说我在 [1.1, 1.1], [1.2, 1.2], [1.3, 1.3] 有 3 分.我想在 x 和 y 轴上使用 1 的宽度使用 1 到 3 进行内核密度估计。

创建 gaussian_kde 时，它​​一直给我这个错误:

raise LinAlgError("singular matrix")
numpy.linalg.linalg.LinAlgError: singular matrix

查看 gaussian_kde 的源代码，我意识到我思考数据集的含义的方式与计算维度的方式完全不同，但是我找不到任何示例代码来显示多维数据如何与模块一起使用。有人可以帮我提供一些使用示例 gaussian_kde多维数据？

Answer 1

我发现很难理解 SciPy 手册中关于 gaussian_kde 的描述适用于二维数据。这是一个旨在补充@endolith 示例的解释。我将代码分成几个带有注释的步骤来解释不太直观的部分。
一、进口:

import numpy as np
import scipy.stats as st
from matplotlib.pyplot import imshow, show

创建一些虚拟数据:这些是“X”和“Y”点坐标的一维数组。

np.random.seed(142)  # for reproducibility
x = st.norm.rvs(loc=2, scale=1, size=2000)
y = st.norm.rvs(loc=0, scale=3, size=2000)

对于二维密度估计， gaussian_kde必须使用包含“X”和“Y”数据集的两行数组初始化对象。在 NumPy 术语中，我们“垂直堆叠它们”:

xy = np.vstack((x, y))

所以“X”数据在第一行 xy[0,:]和“Y”数据在第二行 xy[1,:]和 xy.shape是 (2, 2000) .现在创建 gaussian_kde目的:

dens = st.gaussian_kde(xy)

我们将在二维网格上评估估计的二维密度 PDF。在 NumPy 中创建这样一个网格的方法不止一种。我在这里展示了一种不同于(但在功能上等同于)@endolith 的方法的方法:

gx, gy = np.mgrid[x.min():x.max():128j, y.min():y.max():128j]
gxy = np.dstack((gx, gy)) # shape is (128, 128, 2)

gxy是一个 3-D 数组， [i,j] - gxy 的第一个元素包含对应“X”和“Y”值的 2 元素列表: gxy[i, j]的值为 [ gx[i], gy[j] ] .
我们必须调用 dens() (或 dens.pdf() 是同一件事)在每个二维网格点上。为此，NumPy 有一个非常优雅的函数:

z = np.apply_along_axis(dens, 2, gxy)

换句话说，可调用的 dens (也可能是 dens.pdf)沿 axis=2 调用(第三轴)在 3-D 数组中 gxy并且值应作为二维数组返回。唯一的小故障是 z 的形状将是 (128,128,1)而不是 (128,128)我所期望的。请注意 documentation说:

The shape of out [the return value, L.D.] is identical to the shape of arr, except along theaxis dimension. This axis is removed, and replaced with new dimensionsequal to the shape of the return value of func1d. So if func1d returnsa scalar out will have one fewer dimensions than arr.

最有可能 dens()返回一个 1 长的元组，而不是我希望的标量。我没有进一步调查这个问题，因为这很容易解决:

z = z.reshape(128, 128)

之后我们可以生成图像:

imshow(z, aspect=gx.ptp() / gy.ptp())
show()  # needed if you try this in PyCharm

这是图像。 (请注意，我也实现了 @endolith 的版本，并且得到了一张与这个无法区分的图像。)