使用scipy.stats.gaussian_kde处理二维数据

Question

我正在尝试使用{{link1：scipy.stats.gaussian_kde类}}来平滑一些带有纬度和经度信息的离散数据，以便最终显示为类似于等高线地图的形式，其中高密度是峰值，低密度是山谷。

我很难将二维数据集放入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]。我想使用x轴和y轴宽度为1来进行核密度估计，从1到3进行估计。

创建高斯核密度时，一直出现错误。

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

查看gaussian_kde的源代码后，我意识到我的数据集解释与计算维数方式完全不同，但我找不到任何示例代码展示该模块如何处理多维数据。有人能提供一些使用gaussian_kde处理多维数据的示例方法吗？

- jet

---

尝试使用不全在一行的数据进行测试。我不确定它是否应该因此而失败，或者这是一个错误。 - endolith

Answer 1

这个例子似乎符合您的要求：

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

- endolith

---

当你说轴需要修正时，你指的是什么？因为我正在处理一些数据，但不知何故它会在最小值和最大值之下和之上返回一些多余的数据。 - Srivatsan

@Srivatsan：我想我只是意味着它应该有一个更加正方形的长宽比。 - endolith

Answer 2

我认为你把核密度估计和插值或核回归混淆了。如果你有更多的点样本，KDE可以估计点的分布。

我不确定你需要哪种插值方法，但是在scipy.interpolate中的样条插值或径向基函数插值会更合适。

如果你想要一维核回归，那么你可以在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])

- Josef

---

我不是统计学家，但我对KDE和核回归的阅读以及Jet提到的“等高线图”让我认为KDE是指的内容。 - endolith

Answer 3

我发现理解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)

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

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

所以，“X”数据在第一行xy [0,:] ，“Y”数据在第二行xy [1,:] ， xy.shape 为（2,2000）。现在创建 gaussian_kde 对象：

dens = st.gaussian_kde(xy)

我们将在二维网格上评估估计的二维密度概率分布函数。在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" 是一个三维数组，其中的 "[i,j]" 元素包含相应的 "X" 和 "Y" 值的二元列表："gxy[i, j]" 的值为 "[gx[i], gy[j]]"。

我们必须在每个二维网格点上调用 "dens()"（或者等效的 "dens.pdf()"）函数。NumPy 有一个非常优雅的函数可以实现这个目的。

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

将可调用的`dens`（也可以是`dens.pdf`）应用于三维数组`gxy`的第三个轴`axis=2`，并将值作为二维数组返回。唯一的问题是`z`的形状将是`(128, 128, 1)`而不是我期望的`(128, 128)`。请注意，文档中写道：

输出的形状与输入的形状相同，除了沿着轴的维度。该轴被删除，并替换为func1d返回值的形状相同的新维度。因此，如果func1d返回标量，则out的维数比arr少一个。

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

z = z.reshape(128, 128)

之后，我们可以生成图像：

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

这是图片。（请注意，我也实现了@endolith的版本，并获得了一张与此相似的图片。）

- András Aszódi

Answer 4

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

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

- Abhishek