python - sklearn 中的 2D KDE 带宽与 scipy 中的带宽之间的关系

Question

我正在尝试比较 sklearn.neighbors.KernelDensity 的性能与 scipy.stats.gaussian_kde对于二维数组。

来自 this article我看到带宽 (bw) 在每个函数中的处理方式不同。这篇文章给出了在 scipy 中设置正确 bw 的方法，因此它将等同于 sklearn 中使用的方法。基本上它将 bw 除以样本标准差。结果是这样的:

# For sklearn
bw = 0.15

# For scipy
bw = 0.15/x.std(ddof=1)

其中 x 是我用来获取 KDE 的示例数组。这在 1D 中工作得很好，但我不能让它在 2D 中工作。

这是我得到的MWE:

import numpy as np
from scipy import stats
from sklearn.neighbors import KernelDensity

# Generate random data.
n = 1000
m1, m2 = np.random.normal(0.2, 0.2, size=n), np.random.normal(0.2, 0.2, size=n)
# Define limits.
xmin, xmax = min(m1), max(m1)
ymin, ymax = min(m2), max(m2)
# Format data.
x, y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([x.ravel(), y.ravel()])
values = np.vstack([m1, m2])

# Define some point to evaluate the KDEs.
x1, y1 = 0.5, 0.5

# -------------------------------------------------------
# Perform a kernel density estimate on the data using scipy.
kernel = stats.gaussian_kde(values, bw_method=0.15/np.asarray(values).std(ddof=1))
# Get KDE value for the point.
iso1 = kernel((x1,y1))
print 'iso1 = ', iso[0]

# -------------------------------------------------------
# Perform a kernel density estimate on the data using sklearn.
kernel_sk = KernelDensity(kernel='gaussian', bandwidth=0.15).fit(zip(*values))
# Get KDE value for the point.
iso2 = kernel_sk.score_samples([[x1, y1]])
print 'iso2 = ', np.exp(iso2[0])

(iso2 呈现为指数，因为 sklearn 返回对数值)

iso1 和 iso2 得到的结果不同，我不知道应该如何影响带宽(在任一函数中)以使它们相等(正如他们应该的那样)。

---

添加

我在 sklearn 聊天中(由 ep)被建议我应该在使用 scipy< 计算内核之前缩放 (x,y) 中的值 以获得与 sklearn 可比较的结果。

这就是我所做的:

# Scale values.
x_val_sca = np.asarray(values[0])/np.asarray(values).std(axis=1)[0]
y_val_sca = np.asarray(values[1])/np.asarray(values).std(axis=1)[1]
values = [x_val_sca, y_val_sca]
kernel = stats.gaussian_kde(values, bw_method=bw_value)

即:我在使用 scipy 获取内核之前缩放了两个维度，同时保留了在 sklearn 中获取内核的行。

这给出了更好的结果，但获得的内核仍然存在差异:

其中红点是代码中的(x1,y1)点。因此可以看出，密度估计的形状仍然存在差异，尽管差异很小。也许这是可以实现的最好结果？

Answer 1

几年后，我尝试了这个，并认为我可以在不需要重新缩放数据的情况下使用它。带宽值确实需要一些缩放:

# For sklearn
bw = 0.15

# For scipy
bw = 0.15/x.std(ddof=1)

两个 KDE 对同一点的评估并不完全相等。例如，这是对 (x1, y1) 点的评估:

iso1 =  0.00984751705005  # Scipy
iso2 =  0.00989788224787  # Sklearn

但我想它已经足够接近了。

这是 2D 情况下的 MWE 和输出，据我所知，它们看起来几乎完全相同:

import numpy as np
from scipy import stats
from sklearn.neighbors import KernelDensity
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

# Generate random data.
n = 1000
m1, m2 = np.random.normal(-3., 3., size=n), np.random.normal(-3., 3., size=n)
# Define limits.
xmin, xmax = min(m1), max(m1)
ymin, ymax = min(m2), max(m2)
ext_range = [xmin, xmax, ymin, ymax]
# Format data.
x, y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([x.ravel(), y.ravel()])
values = np.vstack([m1, m2])

# Define some point to evaluate the KDEs.
x1, y1 = 0.5, 0.5
# Bandwidth value.
bw = 0.15

# -------------------------------------------------------
# Perform a kernel density estimate on the data using scipy.
# **Bandwidth needs to be scaled to match Sklearn results**
kernel = stats.gaussian_kde(
    values, bw_method=bw/np.asarray(values).std(ddof=1))
# Get KDE value for the point.
iso1 = kernel((x1, y1))
print 'iso1 = ', iso1[0]

# -------------------------------------------------------
# Perform a kernel density estimate on the data using sklearn.
kernel_sk = KernelDensity(kernel='gaussian', bandwidth=bw).fit(zip(*values))
# Get KDE value for the point. Use exponential since sklearn returns the
# log values
iso2 = np.exp(kernel_sk.score_samples([[x1, y1]]))
print 'iso2 = ', iso2[0]


# Plot
fig = plt.figure(figsize=(10, 10))
gs = gridspec.GridSpec(1, 2)

# Scipy
plt.subplot(gs[0])
plt.title("Scipy", x=0.5, y=0.92, fontsize=10)
# Evaluate kernel in grid positions.
k_pos = kernel(positions)
kde = np.reshape(k_pos.T, x.shape)
plt.imshow(np.rot90(kde), cmap=plt.cm.YlOrBr, extent=ext_range)
plt.contour(x, y, kde, 5, colors='k', linewidths=0.6)

# Sklearn
plt.subplot(gs[1])
plt.title("Sklearn", x=0.5, y=0.92, fontsize=10)
# Evaluate kernel in grid positions.
k_pos2 = np.exp(kernel_sk.score_samples(zip(*positions)))
kde2 = np.reshape(k_pos2.T, x.shape)
plt.imshow(np.rot90(kde2), cmap=plt.cm.YlOrBr, extent=ext_range)
plt.contour(x, y, kde2, 5, colors='k', linewidths=0.6)

fig.tight_layout()
plt.savefig('KDEs', dpi=300, bbox_inches='tight')