使用带有2维数据的scipy.stats.gaussian_kde

Question

我想使用的scipy.stats.gaussian_kde类来平滑经度和纬度信息收集了一些离散数据,所以它显示为有点类似于在最后一个等高线图,其中高密度的峰值密度和低密度的山谷.

我很难将二维数据集放入gaussian_kde课堂.我已经玩弄了解它如何与1维数据一起工作,所以我认为2维将是这样的:

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多维数据吗？

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

Answer 2

我认为您将内核密度估计与插值或内核回归混合在一起。如果您有较大的点样本，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])