Python：正值的核密度估计

Question

我想获得正数据点的核密度估计。使用 Python Scipy Stats 包，我想出了以下代码。

def get_pdf(data):
    a = np.array(data)
    ag = st.gaussian_kde(a)
    x = np.linspace(0, max(data), max(data))
    y = ag(x)
    return x, y

这对于大多数数据集都非常有效，但是对于“所有正数”数据点会给出错误的结果。为确保它正常工作，我使用数值积分来计算这条曲线下的面积。

def trapezoidal_2(ag, a, b, n):
    h = np.float(b - a) / n
    s = 0.0
    s += ag(a)[0]/2.0
    for i in range(1, n):
        s += ag(a + i*h)[0]
    s += ag(b)[0]/2.0
    return s * h

由于数据分布在区域 (0, int(max(data))) 中，因此在执行以下行时，我们应该得到一个接近 1 的值。

b = 1
data = st.pareto.rvs(b, size=10000)
data = list(data)

a = np.array(data)
ag = st.gaussian_kde(a)
trapezoidal_2(ag, 0, int(max(data)), int(max(data))*2)

但当我测试时它给出的值接近 0.5。

但是当我从 -100 积分到 max(data) 时，它提供了一个接近 1 的值。

trapezoidal_2(ag, -100, int(max(data)), int(max(data))*2+200)

原因是，为小于 0 的值定义了 ag (KDE)，即使原始数据集仅包含正值。

那么我怎样才能得到一个只考虑正值的核密度估计，使得区域中曲线下的面积（o，max（数据））接近 1？

Answer 1

在进行核密度估计时，带宽的选择非常重要。我认为斯科特法则和西尔弗曼法则适用于类似于高斯分布的分布。但是，它们不适用于帕累托分布。

引用自doc：

带宽选择会强烈影响从 KDE（比内核的实际形状更重要）。带宽选择 可以通过“经验法则”、交叉验证、“插件”来完成 方法”或其他方式；请参阅 [3]、[4] 了解评论。gaussian_kde 使用经验法则，默认为 Scott's Rule。

尝试不同的带宽值，例如：

import numpy as np
import matplotlib.pyplot as plt

from scipy import stats

b = 1

sample = stats.pareto.rvs(b, size=3000)
kde_sample_scott = stats.gaussian_kde(sample, bw_method='scott')
kde_sample_scalar = stats.gaussian_kde(sample, bw_method=1e-3)


# Compute the integrale:
print('integrale scott:', kde_sample_scott.integrate_box_1d(0, np.inf))
print('integrale scalar:', kde_sample_scalar.integrate_box_1d(0, np.inf))

# Graph:
x_span = np.logspace(-2, 1, 550)
plt.plot(x_span, stats.pareto.pdf(x_span, b), label='theoretical pdf')
plt.plot(x_span, kde_sample_scott(x_span), label="estimated pdf 'scott'")
plt.plot(x_span, kde_sample_scalar(x_span), label="estimated pdf 'scalar'")
plt.xlabel('X'); plt.legend();

给予：

integrale scott: 0.5572130540733236
integrale scalar: 0.9999999999968957

和：

我们看到使用Scott方法的kde是错误的。