scipy gaussian_kde 根据使用的方法产生不同的结果（权重与无权重）

Question

我有一系列坐标要应用 KDE，并且一直在使用 scipy.stats.gaussian_kde 来执行此操作。这里的问题是这个函数需要一组离散的坐标，然后它会执行密度估计。

当我希望记录我的数据时，这会导致问题（对于坐标特别少的集合，并且使用未触及的数据提供的信息很少）。可以想象，如果您必须处理离散数量的点，如果 2 个点出现 18 次，另外 24 次，取 18 和 24 的对数将使它们相同，因为它们必须按顺序四舍五入到最接近的整数保持离散。

作为解决此问题的方法，我一直在 scipy.stats.gaussian_kde 函数中使用 weights 参数。不是每个点出现的次数等于其密度的数组，而是每个点出现一次，而是按其密度加权。所以现在，使用之前的示例，密度为 18 和 24 的 2 个点将不会相同，因为这些密度可以是连续的。

这很有效，并产生了看起来不错的估计值，但是使用这两种不同的方法，它们都会产生具有微小差异的图表。如果我只使用一种方法，我会很高兴地保持无知，但现在我已经使用了两种方法，我无法确定估计是否合理。

这两种方法产生不同结果是否有原因？

请参阅下面的一些重现该问题的示例代码：

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

np.random.seed(1)

discrete_points = np.random.randint(0,10,size=(2,400))

continuous_points = [[],[]]
continuous_weights = []

recorded_points = []
for i in range(discrete_points.shape[1]):
    p = discrete_points[:,i]
    if tuple(p) in recorded_points:
        continuous_weights[recorded_points.index(tuple(p))] += 1
    else:
        continuous_points[0].append(p[0])
        continuous_points[1].append(p[1])
        continuous_weights.append(1)
        recorded_points.append(tuple(p))

resolution = 1

kde = gaussian_kde(discrete_points)
x, y = discrete_points
# https://www.oreilly.com/library/view/python-data-science/9781491912126/ch04.html
x_step = int((max(x)-min(x))/resolution)
y_step = int((max(y)-min(y))/resolution)
xgrid = np.linspace(min(x), max(x), x_step+1)
ygrid = np.linspace(min(y), max(y), y_step+1)
Xgrid, Ygrid = np.meshgrid(xgrid, ygrid)
Z = kde.evaluate(np.vstack([Xgrid.ravel(), Ygrid.ravel()]))
Zgrid = Z.reshape(Xgrid.shape)

ext = [min(x), max(x), min(y), max(y)]
earth = plt.cm.gist_earth_r

plt.imshow(Zgrid,
    origin='lower', aspect='auto',
    extent=ext,
    alpha=0.8,
    cmap=earth)

plt.title("Discrete method (no weights)")
plt.savefig("noweights.png")

kde = gaussian_kde(continuous_points, weights=continuous_weights)
x, y = continuous_points
# https://www.oreilly.com/library/view/python-data-science/9781491912126/ch04.html
x_step = int((max(x)-min(x))/resolution)
y_step = int((max(y)-min(y))/resolution)
xgrid = np.linspace(min(x), max(x), x_step+1)
ygrid = np.linspace(min(y), max(y), y_step+1)
Xgrid, Ygrid = np.meshgrid(xgrid, ygrid)
Z = kde.evaluate(np.vstack([Xgrid.ravel(), Ygrid.ravel()]))
Zgrid = Z.reshape(Xgrid.shape)

ext = [min(x), max(x), min(y), max(y)]
earth = plt.cm.gist_earth_r

plt.imshow(Zgrid,
    origin='lower', aspect='auto',
    extent=ext,
    alpha=0.8,
    cmap=earth)

plt.title("Continuous method (weights)")
plt.savefig("weights.png")

这会产生以下图：

和

Answer 1

kde 的一个重要方面是使用的带宽。 Scipy 的 gaussian_kde 使用 "Scott's factor" 作为带宽的猜测。

特别是，gaussian_kde 使用 n**(-1./(d+4))，其中 d 是维度（在本例中为 2），n

• 非加权版本的数据点数
• 加权版本的“有效数据点数”；计算为neff = sum(weights)^2 / sum(weights^2)

以帖子n = 400 和neff = sum(continuous_weights)**2 / sum([w**2 for w in continuous_weights]) = 84.0336 为例。

要获得相同的结果，应在两种情况下使用相同的带宽。可以显式设置为gaussian_kde(..., bw_method=bandwidth。

bandwidth = discrete_points.shape[1]**(-1./(2+4))

# kde without weights
kde = gaussian_kde(discrete_points, bw_method=bandwidth)

# kde for the weighted points
kde = gaussian_kde(continuous_points, weights=continuous_weights, bw_method=bandwidth)

如果您计划创建多个绘图，您可能希望对所有绘图使用相同的带宽，而与点数或权重无关。您可能想试验分辨率和带宽。更高的带宽可以在更远的距离上平滑所有内容，更小的带宽更忠实于给定的数据。