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 is the bandwidth used. Scipy's gaussian_kde uses "Scott's factor" 的一个重要方面作为对带宽的猜测。

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

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

例子中的post 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)

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

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

scipy gaussian_kde produces different results depending on method used (weights vs no weights)

python

numpy

matplotlib

scipy

kernel-density