获取SciPy的gaussian_kde函数使用的带宽

Question

我正在使用SciPy的stats.gaussian_kde函数，从一组点的数据中生成一个核密度估计（kde）函数。

这是我代码的一个简单示例：

import numpy as np
from scipy import stats

def random_data(N):
    # Generate some random data.
    return np.random.uniform(0., 10., N)

# Data lists.
x_data = random_data(100)
y_data = random_data(100)

# Obtain the gaussian kernel.
kernel = stats.gaussian_kde(np.vstack([x_data, y_data]))

因为我没有手动设置带宽（通过bw_method这个参数），所以这个函数默认使用Scott的规则（可以查看函数的说明）。我需要的是获取这个由stats.gaussian_kde函数自动设置的带宽值。

我尝试过使用：

print kernel.set_bandwidth()

但它总是返回None，而不是一个浮点数。

Answer 1

我想再补充一点，之前@nivniv的回答已经非常棒了。

我在尝试从头开始实现scipy.stats.gaussian_kde这个功能，但我用的样本数量非常少（这可能不太现实），因为我想更好地理解KDE是怎么工作的。起初，我的实现和scipy的结果不一致，但后来我用自由度为1来计算标准差，这才让结果对上了。

# Implementation from scratch
import numpy as np
np.random.seed(12)
Xgrid = np.linspace(-6, 6, 200)
sample = np.random.rand(3) # only 3 samples

scottFactor = len(sample)**(-1/5)
bw = scottFactor * sample.std(ddof=1)
pde_manual = 0.
for s in sample:
    _in = (s - Xgrid) / bw
    pde_manual = pde_manual + _kernel(_in) / bw
pde_manual = pde_manual / len(sample)

为了让你明白我的意思，这里是用scipy构建KDE的代码，

# Scipy
from scipy.stats import gaussian_kde
k_scipy = gaussian_kde(sample, bw_method="scott")
pde_scipy = np.reshape(k_scipy(Xgrid).T, Xgrid.shape)

然后是绘制它们的代码……

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(11, 5), sharey=True)
ax[0].plot(Xgrid, pde_manual)
ax[0].set_title("Implementation from scratch")

ax[1].plot(Xgrid, pde_scipy)
ax[1].set_title("scipy.stats.gaussian_kde")

另外，我还想展示一下如果不设置ddof=1会是什么样子（根据NumPy的文档，ddof=0是默认值）……

不过，如果样本数量多了，ddof的设置就没那么重要了。

Answer 2

我看到这个老问题，因为我也想知道Scipy中的gaussian_kde使用了多少带宽。我想补充一下之前的回答，Scipy的kde.py代码中使用了协方差因子，具体是这样的：

self.covariance = self._data_covariance * self.factor**2

所以，完整的核协方差是样本协方差乘以所谓的协方差因子的平方（也叫Scott因子），这个因子可以通过kde.factor或者kde.covariance_factor()来获取。

Answer 3

我明白了，这一行是：

kernel.covariance_factor()

来自 scipy.stats.gaussian_kde.covariance_factor：

这个函数计算一个系数（kde.factor），它会乘以数据的协方差矩阵，以得到核的协方差矩阵。默认情况下使用的是 scotts_factor。你可以通过子类重写这个方法来提供不同的计算方式，或者通过调用 kde.set_bandwidth 来设置它。

你可以检查一下，使用这个带宽值得到的核和使用默认带宽生成的核是等价的。要做到这一点，你需要用 covariance_factor() 得到一个新的核，然后在一个随机点上比较它的值和原始核的值：

kernel = stats.gaussian_kde(np.vstack([x_data, y_data]))
print kernel([0.5, 1.3])

bw = kernel.covariance_factor()    
kernel2 = stats.gaussian_kde(np.vstack([x_data, y_data]), bw_method=bw)
print kernel2([0.5, 1.3])

Answer 4

简短回答

带宽是通过 kernel.covariance_factor() 乘以你使用的样本的标准差 来计算的。

（这是针对一维样本的情况，默认情况下是根据Scott的经验法则来计算的。）

举个例子：

from scipy.stats import gaussian_kde
sample = np.random.normal(0., 2., 100)
kde = gaussian_kde(sample)
f = kde.covariance_factor()
bw = f * sample.std()

你得到的概率密度函数（pdf）是这样的：

from pylab import plot
x_grid = np.linspace(-6, 6, 200)
plot(x_grid, kde.evaluate(x_grid))

你可以这样检查，如果你使用一个新的函数来创建一个核密度估计（kde），比如说，使用sklearn：

from sklearn.neighbors import KernelDensity
def kde_sklearn(x, x_grid, bandwidth):
    kde_skl = KernelDensity(bandwidth=bandwidth)
    kde_skl.fit(x[:, np.newaxis])
    # score_samples() returns the log-likelihood of the samples
    log_pdf = kde_skl.score_samples(x_grid[:, np.newaxis])
    pdf = np.exp(log_pdf)
    return pdf

现在使用上面的相同代码，你会得到：

plot(x_grid, kde_sklearn(sample, x_grid, f))

plot(x_grid, kde_sklearn(sample, x_grid, bw))