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

Question

我正在尝试使用the scipy.stats.gaussian_kde class 来平滑一些使用纬度和经度信息收集的离散数据，因此它最终显示为有点类似于等高线图，其中高密度是峰值和低密度是山谷。

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

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

我认为您将内核密度估计与插值或内核回归混为一谈。如果您有更大的点样本，KDE 会估计点的分布。

我不确定您想要哪种插值，但 scipy.interpolate 中的样条线或 rbf 会更合适。

如果您想要一维内核回归，那么您可以在 scikits.statsmodels 中找到具有多个不同内核的版本。

更新：这是一个例子（如果这是你想要的）

>>> data = 2 + 2*np.random.randn(2, 100)
>>> kde = stats.gaussian_kde(data)
>>> kde.evaluate(np.array([[1,2,3],[1,2,3]]))
array([ 0.02573917,  0.02470436,  0.03084282])

gaussian_kde 在行中具有变量，在列中具有观察值，因此与统计数据中的通常方向相反。在您的示例中，所有三个点都在一条线上，因此具有完美的相关性。那就是，我猜，奇异矩阵的原因。

调整阵列方向并添加一个小噪点，示例有效，但看起来仍然非常集中，例如您在 (3,3) 附近没有任何采样点：

>>> data = np.array([[1.1, 1.1],
              [1.2, 1.2],
              [1.3, 1.3]]).T
>>> data = data + 0.01*np.random.randn(2,3)
>>> kde = stats.gaussian_kde(data)
>>> kde.evaluate(np.array([[1,2,3],[1,2,3]]))
array([  7.70204299e+000,   1.96813149e-044,   1.45796523e-251])

Answer 2

This example 似乎就是您要找的东西：

import numpy as np
import scipy.stats as stats
from matplotlib.pyplot import imshow

# Create some dummy data
rvs = np.append(stats.norm.rvs(loc=2,scale=1,size=(2000,1)),
                stats.norm.rvs(loc=0,scale=3,size=(2000,1)),
                axis=1)

kde = stats.kde.gaussian_kde(rvs.T)

# Regular grid to evaluate kde upon
x_flat = np.r_[rvs[:,0].min():rvs[:,0].max():128j]
y_flat = np.r_[rvs[:,1].min():rvs[:,1].max():128j]
x,y = np.meshgrid(x_flat,y_flat)
grid_coords = np.append(x.reshape(-1,1),y.reshape(-1,1),axis=1)

z = kde(grid_coords.T)
z = z.reshape(128,128)

imshow(z,aspect=x_flat.ptp()/y_flat.ptp())

显然，轴需要修复。

您还可以使用

绘制数据的散点图

scatter(rvs[:,0],rvs[:,1])

Answer 3

最佳答案中发布的示例对我不起作用。我不得不稍微调整一下，它现在可以工作了：

import numpy as np
import scipy.stats as stats
from matplotlib import pyplot as plt

# Create some dummy data
rvs = np.append(stats.norm.rvs(loc=2,scale=1,size=(2000,1)),
                stats.norm.rvs(loc=0,scale=3,size=(2000,1)),
                axis=1)

kde = stats.kde.gaussian_kde(rvs.T)

# Regular grid to evaluate kde upon
x_flat = np.r_[rvs[:,0].min():rvs[:,0].max():128j]
y_flat = np.r_[rvs[:,1].min():rvs[:,1].max():128j]
x,y = np.meshgrid(x_flat,y_flat)
grid_coords = np.append(x.reshape(-1,1),y.reshape(-1,1),axis=1)

z = kde(grid_coords.T)
z = z.reshape(128,128)

plt.imshow(z,aspect=x_flat.ptp()/y_flat.ptp())
plt.show()

Answer 4

我发现很难理解 SciPy 手册中关于 gaussian_kde 如何处理 2D 数据的描述。这是一个旨在补充 @endolith 示例的解释。我用 cmets 将代码分成几个步骤来解释不太直观的部分。

首先，导入：

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 数组，gxy 的 [i,j]-th 元素包含对应“X”和“Y”值的 2 元素列表：gxy[i, j] 的值为 @ 987654340@.

我们必须在每个二维网格点上调用dens()（或dens.pdf()，这是同一件事）。 NumPy 为此目的提供了一个非常优雅的函数：

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

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

out [返回值，L.D.] 的形状与 arr 的形状相同，除了沿 轴维度。此轴已删除，并替换为新尺寸 等于func1d的返回值的形状。所以如果 func1d 返回 标量的维度将比 arr 少一个。

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

z = z.reshape(128, 128)

之后我们可以生成图像：

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

这是图片。 （请注意，我也实现了 @endolith 的版本，并且得到了一张与这个没有区别的图像。）