原文towardsdatascience.com/the-math-behind-kernel-density-estimation-5deca75cba38?sourcecollection_archive---------1-----------------------#2024-09-17探索核密度估计的基础、概念和数学原理https://medium.com/z.nay854?sourcepost_page---byline--5deca75cba38--------------------------------https://towardsdatascience.com/?sourcepost_page---byline--5deca75cba38-------------------------------- Zackary Nay·发表于Towards Data Science ·14 分钟阅读·2024 年 9 月 17 日–https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/7b02965740dfdbe0f77e416a5b6eb7e3.png核密度估计是一个基本的非参数方法是揭示数据隐藏分布的多功能工具。本文深入探讨了估计器的数学基础提供了选择最优带宽的指导并简要触及了核函数的选择以及其他相关话题。第一部分简介假设我给你以下样本数据https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/9d7b5704844ecc01c61ef6440a044011.png分析样本数据的第一个也是最简单的步骤之一是生成一个直方图对于我们的数据我们得到以下图形https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/52fdb96422c868e48af5ae2cfcb9a638.png并不是非常有用我们距离理解底层分布仍然很远。实际上还有一个额外的问题那就是数据很少表现出由直方图产生的那种锐利的矩形结构。核密度估计提供了解决方法它是一种非参数方法用于估计随机变量的概率密度函数。https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/815e4c24ba9ea67563e2a705b1e06d04.pngimportnumpyasnpimportmatplotlib.pyplotaspltfromsklearn.neighborsimportKernelDensity# Generate sample datanp.random.seed(14)uniform_datanp.random.uniform(low1,high7,size20)normal_datanp.random.normal(loc3,scale0.7,size20)# Combine both distributionscombined_datanp.concatenate([uniform_data,normal_data])# Compute histogramhist,bin_edgesnp.histogram(combined_data,bins9,densityTrue)# Compute KDEkdeKernelDensity(kernelgaussian).fit(combined_data[:,None])xnp.linspace(min(combined_data),max(combined_data),1000)[:,None]log_densitykde.score_samples(x)densitynp.exp(log_density)# Plot histogramplt.hist(combined_data,bins9,densityTrue,colorblue,edgecolorblack,labelHistogram)# Plot KDEplt.plot(x,density,colorred,labelKDE)# Add labels and legendplt.ylabel(Density)plt.title(Histogram and KDE)plt.legend()# Show plotplt.show()第二部分推导以下推导灵感来源于 Bruce E. Hansen 的《非参数讲义笔记》2009 年。如果你有兴趣了解更多可以参考他的原始讲义笔记这里。假设我们想从一组数据样本中估计一个概率密度函数 f(t)。一个好的起点是通过使用 经验分布函数 (EDF) 来估计累积分布函数 F(t)。设 X1, …, Xn 为独立同分布的实值随机变量其共同的累积分布函数为 F(t)。EDF 定义为https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/124f9bd1a658b206a5279dfd27173ddd.png然后根据大数法则当 n 趋近于无穷大时EDF 几乎必然收敛于 F(t)。现在EDF 是一个阶梯函数可能如下所示https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/b8c9f6aac2ad35cb70daede61f52be29.pngimportnumpyasnpimportmatplotlib.pyplotaspltfromscipy.statsimportnorm# Generate sample datanp.random.seed(14)datanp.random.normal(loc0,scale1,size40)# Sort the datadata_sortednp.sort(data)# Compute ECDF valuesecdf_ynp.arange(1,len(data_sorted)1)/len(data_sorted)# Generate x values for the normal CDFxnp.linspace(-4,4,1000)cdf_ynorm.cdf(x)# Create the plotplt.figure(figsize(6,4))plt.step(data_sorted,ecdf_y,wherepost,colorblue,labelECDF)plt.plot(x,cdf_y,colorgray,labelNormal CDF)plt.plot(data_sorted,np.zeros_like(data_sorted),|,colorblack,labelData points)# Label axesplt.xlabel(X)plt.ylabel(Cumulative Probability)# Add gridplt.grid(True)# Set limitsplt.xlim([-4,4])plt.ylim([0,1])# Add legendplt.legend()# Show plotplt.show()因此如果我们尝试通过取经验分布函数EDF的导数来找到 f(t) 的估计器我们将得到一个 狄拉克 delta 函数 的加权和这并没有太大帮助。相反让我们考虑使用估计器的 t二点中心差分公式作为导数的近似。对于一个小的 h0我们得到https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/59ae71a5b9101204a61466fda298a3b7.png现在定义函数 k(u) 如下https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/39d48ca6e7d2b0c4db071cd0eaddc43c.png然后我们得到https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/130318df46d4d42879bfb8b4ef0d52e2.png这是核密度估计器的特例其中 k 是均匀核函数。更一般地核函数是一个从实数到实数的非负函数满足以下条件https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/d490e3f1e63f3b07459eb3f97598cde8.png我们假设本文讨论的所有核函数都是对称的因此我们有 k(-u) k(u)。核的矩它为核函数的形状和行为提供了洞察定义如下https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/3a624958e05b0beb451fdc6b1e39836a.png最后核的阶数定义为第一个非零矩。我们只能通过改变 h 值带宽或核函数来最小化核密度估计器的误差。带宽参数对结果估计的影响远大于核函数但选择起来也更加困难。为了演示 h 值的影响考虑以下两个核密度估计。使用高斯核来估计从标准正态分布生成的样本两个估计器之间唯一的区别就是所选择的 h 值。https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/8202c10fcb27e1b851f069889558d274.pngimportnumpyasnpimportmatplotlib.pyplotaspltfromscipy.statsimportgaussian_kde# Generate sample datanp.random.seed(14)datanp.random.normal(loc0,scale1,size100)# Define the bandwidthsbandwidths[0.1,0.3]# Plot the histogram and KDE for each bandwidthplt.figure(figsize(12,8))plt.hist(data,bins30,densityTrue,colorgray,alpha0.3,labelHistogram)xnp.linspace(-5,5,1000)forbwinbandwidths:kdegaussian_kde(data,bw_methodbw)plt.plot(x,kde(x),labelfBandwidth {bw})# Add labels and titleplt.title(Impact of Bandwidth Selection on KDE)plt.xlabel(Value)plt.ylabel(Density)plt.legend()plt.show()差异相当显著。现在让我们看看在保持带宽不变的情况下改变核函数的影响。https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/c99738f2f355e32578811f64b2dd879c.pngimportnumpyasnpimportmatplotlib.pyplotaspltfromsklearn.neighborsimportKernelDensity# Generate sample datanp.random.seed(14)datanp.random.normal(loc0,scale1,size100)[:,np.newaxis]# reshape for sklearn# Intialize a constant bandwidthbandwidth0.6# Define different kernel functionskernels[gaussian,epanechnikov,exponential,linear]# Plot the histogram (transparent) and KDE for each kernelplt.figure(figsize(12,8))# Plot the histogramplt.hist(data,bins30,densityTrue,colorgray,alpha0.3,labelHistogram)# Plot KDE for each kernel functionxnp.linspace(-5,5,1000)[:,np.newaxis]forkernelinkernels:kdeKernelDensity(bandwidthbandwidth,kernelkernel)kde.fit(data)log_densitykde.score_samples(x)plt.plot(x[:,0],np.exp(log_density),labelfKernel {kernel})plt.title(Impact of Different Kernel Functions on KDE)plt.xlabel(Value)plt.ylabel(Density)plt.legend()plt.show()虽然在尾部视觉上有很大的差异但不同核函数的估计器整体形状相似。因此我将主要集中于为估计器找到最优带宽。现在让我们探讨一些核密度估计器的性质包括其偏差和方差。第三部分核密度估计器的性质我们需要利用的第一个事实是估计器在实数线上的积分为 1。为了证明这一事实我们将使用 u 代换https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/e33cf2e2049c29aca18e8a093ca66c13.png运用 u 代换我们得到以下结果https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/9f3a57706737f576377668a492d07d83.png现在我们可以找到估计密度的均值https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/ed503405029561f4a7976cd6a585ae03.png因此估计密度的均值就是样本均值。现在让我们找到估计量的第二阶矩。https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/bfe471ec7996859cc904c593b9a23c86.png然后我们可以找到估计密度的方差https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/2c822ab90eee6dff0793fa3b7c7ddaa8.png为了找到最优带宽和核函数我们将最小化估计量的均方误差。为此我们首先需要找到估计量的偏差和方差。随机变量 X 的期望值具有概率密度函数 f可以通过以下公式计算https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/1675a80c94a4c62e844235a9ee73b25f.png因此https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/3f155d13144fd99e7f129ce6610dc563.png其中 z 是一个虚拟的积分变量。然后我们可以使用变量代换得到以下结果https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/c0a469d13dbaa87efedfbe3e45fac73f.png因此我们得到以下结果https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/33a5697549c9c8fc7aef8cc5cc5ac786.png然而大多数情况下上述积分无法解析求解。因此我们将通过使用泰勒展开来近似 f(x hu)。作为提醒f(x) 关于 a 的泰勒展开式为https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/e21d33426eddb14757b89477692c5e9c.png因此假设 f(x) 对 v1 次项可微。对于 f(xhu)展开式为https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/b0e29e10c63fa4e7688400f1b90e4c1e.png然后对于 v 阶核函数我们可以将表达式展开到 v’项https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/e062240689da239fe38625a7e2a75c84.png其中最后一项是余项当 h 趋近于 0 时该项的消失速度比 h 快。现在假设 k 是 v 阶核函数我们得到https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/f7854ca857959f4e956f84b66bff38be.png因此https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/9618d12cb67889aa94553603e126194a.png因此估计量的偏差为https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/aeab30eaee8173e839d4b93f05452fbf.png方差的上界可以通过以下计算获得 [Chen 2].https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/4082e2eee9c96f468b9e2bd8e068cdae.png其中https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/349e7cdc90db769123596f25b306f7f0.png这个项也被称为粗糙度可以表示为 R(k)。最后我们可以得到均方误差的表达式https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/9c1d3c9c4559f6512c1fd4194a06d4ea.pngAMSE 是渐近均方误差Asymptotic Mean Squared Error的缩写。我们可以最小化渐近均积分平方误差该误差定义如下https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/02ada6a1f661f96e4fbc20bc99bbf909.png以找到能带来最小误差的带宽Silverman, 1986https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/c473fd6806cfa8781ba6b8a6a67c3ebb.png然后通过将方程设置为 0 并简化得到一个核函数的最优带宽为https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/b95e30c3665889767077d44b94e3e699.png更为熟悉的表达式适用于二阶核其中最小化 AMISE 的带宽为https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/53599fe5eeb322d7a7b1003ac5a7276d.png然而这个解法可能令人失望因为我们需要知道我们正在估计的分布才能找到最优的带宽。实际上如果我们一开始就使用核密度估计器我们就不会知道这个分布。第四部分带宽选择尽管无法找到最小化均积分平方误差的带宽但仍有几种方法可以在不知道基础分布的情况下选择带宽。然而值得注意的是较大的 h 值会导致估计量的方差较小但偏差较大而较小的带宽会产生一个粗略的估计偏差较小Eidous 等2010。一些用于找到带宽的方法包括使用1: 解方程规则2: 最小二乘交叉验证3: 偏差交叉验证4: 直接代入法5: 对比法值得注意的是根据样本大小和数据的偏斜程度选择带宽的“最佳”方法会有所不同。Python 中的大多数软件包允许你使用 Silverman 提出的法则你可以直接代入某个分布通常是正态分布来表示 f然后根据你选择的核函数计算带宽。这一过程被称为Silverman 经验法则它为带宽提供了一个相对简单的估计。然而它往往会高估带宽导致较平滑的估计方差较低但偏差较大。Silverman 经验法则对于双峰分布的表现也不好对于这种情况有更精确的技术可供选择。如果假设数据服从均值为 0 的正态分布使用样本标准差并应用 Epanechnikov 核下文讨论你可以通过以下方程使用经验法则选择带宽https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/15398a17b0be9cacc757582fe73205c7.pngEidous 等人发现与我列出的其他方法相比对比法表现最佳。然而这种方法有其缺点因为它增加了需要选择的参数数量。交叉验证方法是带宽选择的另一种不错选择因为它通常会导致较小的偏差但较大的方差Heidenreich 等2013。它在寻找倾向于低估最佳带宽的带宽选择器时非常有用。为了避免这篇文章过长我将只介绍普通最小二乘交叉验证方法。该方法对于较为弯曲的密度和样本量大约为 50 到 100 的情况效果较好。如果样本量非常小或非常大这篇论文是一个很好的资源可以帮助你找到另一种选择带宽的方法。如 Heidenreich 所指出的“选择哪个带宽选择器确实有很大不同不仅在数值上而且对密度估计的质量也有影响”。为了快速复习当我们创建模型时我们会保留一部分样本作为验证集以查看模型在未训练的样本上的表现。K 折交叉验证通过将数据集分成 K 部分并训练模型 K 次从而最小化可能选择错误测试集的影响每一部分都会轮流作为验证集一次。留一法交叉验证Leave One Out Cross Validation是 K 折交叉验证的极端形式对于样本量为 n 的情况我们训练模型 n 次每次留出一个样本。这里有一篇文章详细探讨了这种方法在训练机器学习算法中的应用。让我们回到 AMISE 表达式。我们将不再研究渐近平均积分平方误差而是最小化平均积分平方误差MISE。首先我们可以展开平方https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/8d07064b850466ae87e1c1924a01c17d.png由于我们无法控制最后一项的期望值因此我们可以尽力最小化前两项并舍弃第三项。然后由于 X 是 E[X] 的无偏估计量我们可以直接求得第一项https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/5e82bb063a7a2eccf4342d6c403d000a.png然后k 与 k 的卷积定义如下https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/57fa4edb38fda6db5b996091bd9dd4df.png因此https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/bba7957ebdfc234384c04e6ebd32959b.png因此对于第一项我们得到https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/dfba3c5255b22e2d8c8b1b0a13996df6.png接下来我们可以通过蒙特卡洛方法来逼近第二项。首先正如前面讨论的密度函数等于累积分布函数的导数我们可以通过经验分布函数来逼近它。然后积分可以通过估计量在样本上的平均值来逼近。https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/f96f9260cac17a62a41b9d40be4bb682.png然后最小二乘交叉验证选择器LSCV被定义为如下https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/d4fc33d8d1ce23ecd2bd8b275c7012fc.pnghttps://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/f372f058f5adfa4cc5fd14285e1ffcbb.png我们得到最终定义的选择器如下https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/2d69b80d31ae868cc3d74514a1728103.png最优带宽是最小化 LSCV 的 h 值定义如下https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/3088842523b4847cd1a05dd2477eb061.pngLSCV(h)函数可能有多个局部最小值因此你找到的最优带宽可能对选择的区间非常敏感。绘制该函数并直观地检查全局最小值所在位置是很有用的。第五部分最优内核选择如果我们使用的是二阶内核这很常见则内核选择比带宽选择要简单得多。在最小化 AMISE 的标准下Epanechnikov**内核是最优的内核。**完整的证明可以在Muler的论文中找到。https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/cfe9062c4a6e3b26a7b94a474d190507.png还有一些内核与 Epanechnikov 内核一样高效这些也在 Muler 的论文中有所提及。然而我认为你不必过于担心内核函数的选择带宽的选择要重要得多。第六部分进一步的主题与结论自适应带宽提出的一种改进核密度估计器的方法是通过自适应带宽。自适应带宽会根据每个数据点的密度调整带宽当样本数据的密度较低时增加带宽密度较高时减少带宽。虽然理论上很有前景但在一维情况下自适应带宽的表现却非常差Terrel, Scott 1992。虽然它在高维空间中可能表现更好但在一维情况下我认为坚持使用常数带宽是最好的选择。多变量核密度估计核密度估计器还可以扩展到更高维度其中内核是径向基函数或是多个内核函数的乘积。这种方法确实存在维度灾难的问题随着维度的增加所需数据点的数量会呈指数增长。它计算成本高昂对于高维数据分析并不是一个理想的方法。非参数多变量密度估计仍然是一个非常活跃的领域Masked Autoregressive Flow似乎是一个全新且有前景的方法。现实世界应用核密度估计在各个学科中有着广泛的应用。首先它已被证明能改善机器学习算法例如在灵活的朴素贝叶斯分类器中。它还被用于估算交通事故密度。该链接的论文使用 KDE 帮助构建模型指示日本不同城市的交通事故风险。另一个有趣的应用是在地震学y 中KDE 被用于建模不同地点的地震风险。结论核密度估计器是数据分析师工具箱中的一个优秀补充。虽然直方图无疑是分析数据的一种良好方式它不依赖任何假设但核密度估计器为单变量数据提供了一个稳固的替代方案。对于高维数据或当计算时间是一个问题时我建议考虑其他方法。尽管如此KDE 仍然是数据分析中直观、强大且多用途的工具。除非另有注明所有图片均为作者提供。