核密度估计是一种估计随机变量概率密度函数的非参数方式。本示例使用一维数据来展示核密度估计的基本概念。首先，来看直方图在可视化一维数据点密度时遇到的问题。直方图可以想象为在规则网格上每个点上方堆叠一个单位“块”。然而，如上图所示，这些块的网格选择可能导致对密度分布形状的截然不同的理解。如果将每个块的中心放在它所代表的点上，得到了左下角图所示的估计。这是一个使用“顶帽”核的核密度估计。这个想法可以推广到其他核形状：第一张图的右下角图显示了相同分布上的高斯核密度估计。
Scikit-learn库通过KernelDensity估计器实现了高效的核密度估计，可以使用Ball Tree或KD Tree结构。本示例的第二张图展示了可用的核。第三张图比较了一维100个样本分布的核密度估计。尽管本示例使用的是一维分布，但核密度估计也很容易且高效地扩展到更高维度。
以下是使用Python的matplotlib和scikit-learn库实现一维核密度估计的代码示例。首先，导入必要的库，然后生成一些随机数据来模拟一维分布。接着，使用KernelDensity类来拟合这些数据，并使用不同的核函数来估计密度。最后，将结果可视化出来。
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
from sklearn.neighbors import KernelDensity
# 设置随机种子以获得可重复的结果
np.random.seed(1)
# 生成模拟数据
N = 100
X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)), np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
# 绘制直方图和核密度估计
fig, ax = plt.subplots(2, 2, sharex=True, sharey=True)
fig.subplots_adjust(hspace=0.05, wspace=0.05)
# 绘制直方图
ax[0, 0].hist(X[:, 0], bins=np.linspace(-5, 10, 10), fc="#AAAAFF", density=True)
ax[0, 0].text(-3.5, 0.31, "直方图")
# 绘制偏移的直方图
ax[0, 1].hist(X[:, 0], bins=np.linspace(-5, 10, 10) + 0.75, fc="#AAAAFF", density=True)
ax[0, 1].text(-3.5, 0.31, "直方图，偏移的bins")
# 绘制顶帽核密度估计
kde = KernelDensity(kernel="tophat", bandwidth=0.75).fit(X)
log_dens = kde.score_samples(X_plot)
ax[1, 0].fill(X_plot[:, 0], np.exp(log_dens), fc="#AAAAFF")
ax[1, 0].text(-3.5, 0.31, "顶帽核密度")
# 绘制高斯核密度估计
kde = KernelDensity(kernel="gaussian", bandwidth=0.75).fit(X)
log_dens = kde.score_samples(X_plot)
ax[1, 1].fill(X_plot[:, 0], np.exp(log_dens), fc="#AAAAFF")
ax[1, 1].text(-3.5, 0.31, "高斯核密度")
# 绘制所有可用的核
X_plot = np.linspace(-6, 6, 1000)[:, None]
X_src = np.zeros((1, 1))
fig, ax = plt.subplots(2, 3, sharex=True, sharey=True)
fig.subplots_adjust(left=0.05, right=0.95, hspace=0.05, wspace=0.05)
def format_func(x, loc):
if x == 0:
return "0"
elif x == 1:
return "h"
elif x == -1:
return "-h"
else:
return "%ih" % x
for i, kernel in enumerate(["gaussian", "tophat", "epanechnikov", "exponential", "linear", "cosine"]):
axi = ax.ravel()[i]
log_dens = KernelDensity(kernel=kernel).fit(X_src).score_samples(X_plot)
axi.fill(X_plot[:, 0], np.exp(log_dens), "-k", fc="#AAAAFF")
axi.text(-2.6, 0.95, kernel)
axi.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
axi.xaxis.set_major_locator(plt.MultipleLocator(1))
axi.yaxis.set_major_locator(plt.NullLocator())
axi.set_ylim(0, 1.05)
axi.set_xlim(-2.9, 2.9)
ax[0, 1].set_title("可用核")
# 绘制一维密度示例
N = 100
np.random.seed(1)
X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)), np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
true_dens = 0.3 * norm(0, 1).pdf(X_plot[:, 0]) + 0.7 * norm(5, 1).pdf(X_plot[:, 0])
fig, ax = plt.subplots()
ax.fill(X_plot[:, 0], true_dens, fc="black", alpha=0.2, label="输入分布")
colors = ["navy", "cornflowerblue", "darkorange"]
kernels = ["gaussian", "tophat", "epanechnikov"]
lw = 2
for color, kernel in zip(colors, kernels):
kde = KernelDensity(kernel=kernel, bandwidth=0.5).fit(X)
log_dens = kde.score_samples(X_plot)
ax.plot(X_plot[:, 0], np.exp(log_dens), color=color, lw=lw, linestyle="-", label="kernel = '%s'" % kernel)
ax.text(6, 0.38, "N=%d points" % N)
ax.legend(loc="upper left")
ax.plot(X[:, 0], -0.005 - 0.01 * np.random.random(X.shape[0]), "+k")
ax.set_xlim(-4, 9)
ax.set_ylim(-0.02, 0.4)
plt.show()
上述代码首先设置了随机种子以确保结果的可重复性。然后，生成了一些模拟数据，这些数据由两个正态分布混合而成。接着，使用KernelDensity类来拟合这些数据，并使用不同的核函数来估计密度。最后，将结果可视化出来，包括直方图、顶帽核和高斯核密度估计，以及所有可用核的展示。