7

我正在尝试使用该类来平滑一些使用纬度和经度信息收集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])

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

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

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

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

4

4 回答 4

7

这个例子似乎是你正在寻找的：

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

在此处输入图像描述

于 2011-05-25T14:55:52.200 回答
4

我认为您将内核密度估计与插值或内核回归混为一谈。如果您有更大的点样本，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])
于 2010-11-09T00:28:43.233 回答
1

我发现很难理解 SciPy 手册中关于如何gaussian_kde处理 2D 数据的描述。这是一个旨在补充@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（也可以）在 3-D 数组中dens.pdf沿axis=2（第三个轴）被调用，gxy并且值应该作为 2-D 数组返回。唯一的问题是zwill的形状(128,128,1)不是(128,128)我所期望的。请注意，文档说：

out [返回值，LD] 的形状与 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 的版本，并且得到了一张与这个无法区分的图像。）

上述命令的输出

于 2020-10-23T12:13:17.690 回答
0

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

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()
于 2017-09-19T16:11:58.100 回答