核密度估计（KDE）（一）

核密度估计（KDE）（一）

1. 算法简介

1.1 核密度估计算法介绍

设 $\mathcal{S} = (x^{(1)}, x^{(2)}, ..., x^{(n)})$ 是从总体中抽取的一组独立同分布（i.i.d.）样本数据，在统计分析中，我们通常需要通过这组有限样本，推断出总体的概率分布特征，而核密度估计正是实现这一目标的非参数方法.

1.2 核密度估计算法基本思路

核密度估计的核心思想是：每个样本点都会对其周围区域的概率密度产生一定的 “影响”，将所有样本点的贡献叠加起来，就能得到整个数据空间的概率密度分布，具体可拆解为以下两个步骤：

1. 为每个样本点分配一个“影响衰减函数”：它用于描述单个样本点对周围区域密度的“影响范围”和“影响强度”. 距离样本点越近的位置受到的影响越强，函数值越大；反之受到的影响越弱，函数值越小.

2. 叠加各样本的“影响衰减函数”：将每个样本点对应的缩放核函数在整个数据空间中进行叠加，叠加后的结果就是核密度估计得到的总体概率密度函数.

上述“影响衰减函数”通常用核函数来定义，因此本算法称为“核密度估计”. 从直观上理解，核密度估计相当于用无数个 “小的概率密度峰”（每个样本点对应的核函数）共同构建出一个 “整体的概率密度山”，这个 “山” 的形状完全由样本点的分布和核函数的特性决定，既能避免参数估计对分布假设的依赖，又能有效平滑数据中的噪声.

2. 核函数

2.1 核函数的定义

在介绍核密度估计前我们先引入对核函数的介绍.

$\color{blue}{\textbf{Def 2.1 核函数: }}$ 一个核函数 $K: \mathbb{R} \to [0, +\infty)$ 是一个实值函数，满足以下三条性质：

1. 非负性：$K(u) \geq 0$ 对任意的 $u \in \mathbb{R}$；
2. 归一性：$\int_{\mathbb{R}} K(u) \mathrm{d}u = 1$；
3. 对称性：$K(-u) = K(u)$ 对任意的 $u \in \mathbb{R}$.

非负性确保核函数对周围区域的密度贡献是正向的；归一性保证单个样本点对整体密度的总贡献为 $1$；对称性则确保样本点对其左侧、右侧区域的密度贡献对称，避免引入方向偏差.

以下给出了常见的核函数及其公式：

核函数	公式
高斯核（Gaussian）	$K(u) = \frac{1}{\sqrt{2 \pi}} \exp(-\frac{u^2}{2})$
Epanechnikov 核	$K(u) = \begin{cases} \frac{3}{4}(1 - u^2), & \text{if } \vert u \vert < 1 \\ 0, & \text{otherwise}\end{cases}$
均匀核（Uniform / Rectangular）	$K(u) = \begin{cases} \frac{1}{2}, & \text{if } \vert u \vert < 1 \\ 0, & \text{otherwise}\end{cases}$
三角核（Triangular）	$K(u) = \begin{cases} 1 - \vert u \vert, & \text{if } \vert u \vert < 1 \\ 0, & \text{otherwise}\end{cases}$
双权核（Biweight / Quartic）	$K(u) = \begin{cases}\frac{15}{16}(1 - u^2)^2, & \text{if } \vert u \vert < 1 \\ 0, & \text{otherwise}\end{cases}$
三权核（Triweight）	$K(u) = \begin{cases} \frac{35}{32}(1 - u^2)^3, & \text{if } \vert u \vert < 1 \\ 0, & \text{otherwise}\end{cases}$

2.2 缩放核函数的定义与应用

在上述核函数中存在带宽的概念. 以 Epanechnikov 核为例，其默认带宽（bandwidth）为 $h = 1$，即对 $\vert u \vert < 1$ 的区域产生影响，对 $\vert u \vert \geq 1$ 的区域不产生影响. 这在实际应用中会产生“尺度不匹配”的问题，即样本的分布范围可能与带宽不符.

例如：某身高样本的相邻间隔约为 $5$（单位：厘米，下同），而 Epanechnikov 核函数的带宽默认为 $1$，导致核函数的影响范围无法覆盖相邻样本，进而使密度曲线"断裂"；为解决这一问题，我们可以缩放带宽至 $h = 5$，将核函数影响范围扩展到 $\pm 5$，从而有效衔接相邻样本的贡献.

$\color{blue}{\textbf{Def 2.2 缩放核函数: }}$ 设 $K: \mathbb{R} \to [0, +\infty)$ 是一给定的核函数，对带宽 $h > 0$，定义其缩放核函数：

显然，$\int_\mathbb{R} K_h(u) \mathrm{d}u = \frac{1}{h} \int_\mathbb{R} K(\frac{u}{h}) \mathrm{d}u = 1$，即缩放后的核函数不影响核函数的基本性质.

3. 核密度估计算法

3.1 从直方图估计到核密度估计

对于该类问题，最值观的解法是“直方图法”. 例如，设有 $n$ 个 $1$ 维样本点 $x^{(1)}, x^{(2)}, \dots, x^{(n)}$，$x^{(i)} \in [a, b]$，其分布服从概率密度函数 $p(x)$. 作 $[a, b]$ 的等距分划 $T: a = a_0 < a_1 < a_2 < \dots < a_n = b$，$a_i = a_0 + \frac{i}{n}$，$\parallel T \parallel = a_k - a_{k - 1} = \frac{b - a}{n}$. 则可以估计密度函数：

记 $p_k = \mathbb{P} \left[ x^{(i)} \in [a_{k - 1}, a_k) \right]$，若 $p \in C\left( \mathbb{R} \right)$，则：

计算 $\hat{p}$ 的期望为：

计算 $\hat{p}$ 的方差为：

当 $\parallel T \parallel \rightarrow 0^+$ 时，$\xi_k \rightarrow x, \space x \in [a_{k - 1}, a_k)$，有：

即：$\hat{p}$ 是 $p$ 的无偏估计. 但注意到此时 $\mathrm{Var} \left[ \hat{p}(x) \right]$ 会增大，若要减小方差，还需保证 $n \parallel T \parallel \rightarrow +\infty$，即 $n \gg \frac{1}{h}$.

显然，该方法存在以下缺陷：

1. 区间的位置和宽度会严重影响估计结果：例如，改变区间的端点或改变分划 $T$ 都会对结果产生影响；
2. 密度函数不连续：$\hat{p}(x)$ 不连续，无法反映数据分布的平滑性；
3. 难以推广到高维空间：当数据维度增加时，区间数量会呈指数级增长；且大多数区间内无样本点，致使其密度估计为 $0$.

核密度估计正是为解决这些问题而提出的：它将直方图的 “区间贡献” 替换为 “样本点的核函数贡献”，用连续的核函数替代离散的区间，从而实现平滑、连续的密度估计.

3.2 核密度估计算法推导

首先我们解决第一个问题：摆脱对 $T$ 的依赖，其核心思想是不再将 $x$ 与 $T$ 作关联，而是直接计算 $x$ 附近的点的密度. 设总体的随机变量为 $X$，其概率分布函数为 $F(x) = P(X \leq x)$，概率密度函数为 $f(x)$，则有：

根据大数定律，当样本空间足够大时，可以用频率近似概率，得到：

将 $(8)$ 带入 $(7)$，得：

其中 $K(u) = \frac{1}{2} \mathbb{1} ( \vert u \vert \leq 1)$，即均匀核函数. 下面我们验证 $f(x)$ 是概率密度函数，事实上：

为解决第二个问题，只需将不连续的均匀核函数替换为其它连续核函数即可. 根据归一性，替换后得到的 $\hat{f}(x)$ 仍为概率密度函数. 在此基础上，问题三也得以解决，我们称这样的方法为核密度估计.

3.3 核密度估计算法

$\color{blue}{\textbf{Algorithm: Kernel Density Estimation}}$ 设 $\left( x^{(1)}, x^{(2)}, \dots, x^{(n)} \right)$ 是一组给定的独立同分布的简单随机样本，$K(x)$ 是给定的核函数，$h > 0$ 为给定带宽，则核密度估计给出的概率密度估计函数为：

其中核函数 $K(x)$ 即为样本点的“影响衰减函数”，不同样本点的影响衰减函数彼此相同. 通过选取不同的核函数和带宽，可以改变样本点的影响衰减效应，从而改变概率密度估计. 本文剩余内容将介绍如何选取合适的核函数与带宽，但在此之前，我们先讨论核密度估计的期望、方差及其渐进性质. 为此，我们给出三条假设：

1. 假设 $f \in C^2(\mathbb{R})$ 且其二阶导数 $\nabla^2f$ 平方可积，记为 $R(\nabla^2 f) = \int_\mathbb{R} \left [ \nabla^2f(x) \right]^2 \mathrm{d}x$；

2. 假设核函数 $K$ 二阶矩存在且平方可积，分别记为$\mu_2(K) = \int_{\mathbb{R}} u^2 K(u) \mathrm{d}u$，$R(K) = \int_{\mathbb{R}} K^2(u) \mathrm{d}u$；

3. 记 $h = h_n$，假设 $n \rightarrow \infty$ 且 $h \rightarrow 0^+$ 时满足 $nh \rightarrow +\infty$.

计算 $\hat{f}$ 的期望为：

对 $f(x - hu)$ 作二阶 Taylor 展开：

带入到 $\mathbb{E}\left[ \hat{f}(x) \right]$ 中得到渐进等式：

注意到 $f(x - hu) = f(x) + O(hu)$，计算 $\hat{f}$ 的方差为：

其中用到了如下渐进等式：

$\color{blue}{\textbf{Prop 3.1: }}$ $\int_{\mathbb{R}} K^2(u) O(hu) \mathrm{d}u = O(h)$.

$\color{brown}{\textbf{Proof: }}$ 对于 $O(z)$，存在 $C > 0$，使得 $\vert O(z) \vert < C \vert z \vert$，于是有：

最后我们给出 $\mathrm{Bias} \left[ \hat{f}(x) \right]$ 的渐进等式：

可以发现，$\mathrm{Bias} \left[ \hat{f}(x) \right]$ 与 $h^2$ 和 $\nabla^2 f(x)$ 有关. 在 $f(x)$ 为凸函数的区域，$\nabla^2 f(x) > 0$，有 $\mathrm{Bias} \left[ \hat{f}(x) \right] > 0$，即核密度估计会高估 $f$；反之在 $f(x)$ 为凹函数的区域，$\nabla^2 f(x) < 0$，有 $\mathrm{Bias} \left[ \hat{f}(x) \right] < 0$，即核密度估计会低估 $f$.

4. 带宽的选择

4.1 带宽评价指标

要选取合适的带宽，首先要给出带宽的评价指标. 我们发现当带宽 $h$ 增大时，影响衰减函数图像越“扁”、曲线越平滑，此时方差减小、偏差增大；反之当带宽 $h$ 减小时，影响衰减函数图像越“尖”、曲线越陡峭，此时方差增大、偏差减小. 因此，要评价带宽，必须给出对方差和偏差的综合评价指标.

$\color{blue}{\textbf{Def 4.1 均方误差: }}$ 函数 $\mathrm{MSE} \left[ \hat{f}(x) \right] = \mathbb{E}\left[ \left( \hat{f}(x) - f(x) \right)^2 \right]$ 称为均方误差（Mean Squared Error）.

$\color{blue}{\textbf{Prop 4.2: }}$ $\mathrm{MSE} \left[ \hat{f}(x) \right] = \mathrm{Var} \left[ \hat{f}(x) \right] + \mathrm{Bias}^2 \left[ \hat{f}(x) \right]$.

$\color{brown}{\textbf{Proof: }}$ 易知 $\mathbb{E} \left[ \hat{f}(x) - \mathbb{E} \left( \hat{f}(x) \right) \right] = 0$，则有：

$\color{blue}{\textbf{Def 4.3 均方积分误差: }}$ 积分值 $\mathrm{MISE} \left[ \hat{f}(x) \right] = \int_\mathbb{R} \mathrm{MSE} \left[ \hat{f}(x) \right] \mathrm{d}x$ 称为均方积分误差（Mean Integrated Squared Error）.

最优带宽 $h_{\mathrm{opt}}$ 理论上是使均方积分误差最小的带宽，即：

4.2 渐进最优带宽

首先给出 KDE 中 $\mathrm{MSE}$ 和 $\mathrm{MISE}$ 的渐进形式：

$\color{blue}{\textbf{Prop 4.4: }}$ $\hat{f}(x) \xrightarrow{2} f(x)$.

$\color{brown}{\textbf{Proof: }}$ 只需注意到 $\mathrm{MSE} \left[ \hat{f}(x) \right] = O(\frac{1}{nh}) + O(h^4) = o(1)$. $\square$

$\color{blue}{\textbf{Cor 4.5: }}$ $\hat{f}(x) \xrightarrow{2, \mathbb{P}, d} f(x)$.

下面讨论渐进意义下求解最优带宽，即：

对 $\mathrm{AMISE} \left[ \hat{f}(x) \right]$ 关于 $h$ 求导，得：

令 $\frac{\mathrm{d}}{\mathrm{d}h} \mathrm{AMISE} \left[ \hat{f}(x) \right] = 0$，有：

在实际应用中，由于不清楚 $f$ 的真实面貌，无法直接求解 $h_{\mathrm{asym-opt}}$. 下面给出一些常用方法.

4.3 Silverman 方法

Silverman 方法假定 $f$ 服从正态分布 $\mathcal{N}(\mu, \sigma^2)$，即：

易知：

求解 $R(\nabla^2 f)$：

记 $C = \sqrt[5]{\frac{8 \sqrt{\pi} R(K)}{3 \mu_2^2(K)}}$，于是有：

最后，我们对 $\sigma$ 进行估计. 最直接的方法是使用样本的标准差 $s$ 估计 $\sigma$：

这是因为 $\mathbb{E} \left[ s^2 \right] = \sigma^2$. 然而，如果样本点中包含极端值，会导致 $\sigma$ 被严重高估. 为解决这一问题，我们考虑变换 $X = \mu + \sigma Z, \space Z \sim \mathcal{N}(0, 1)$，再将样本点进行排序，得到其上四分位数 $Q_{0.75}$ 和下四分位数 $Q_{0.25}$. 根据正态分布的性质，可得：

其中 $\Phi = \frac{1}{\sqrt{2 \pi}} \exp(-\frac{1}{2}x^2)$ 是标准正态函数. 两式相减，得：

尽管 $\hat{\sigma}_Q$ 相较于 $\hat{\sigma}_s$ 更能避免极端值的影响，但其之运用到了 $50\%$ 的数据，信息利用不完整. 在数据符合正态分布时，$\hat{\sigma}_s$ 是最优估计，$\hat{\sigma}_Q$ 会略大；而数据中出现极端值时，$\hat{\sigma}_s$ 会大幅增大，而 $\hat{\sigma}_Q$ 不受影响. 基于此，我们给出如下对 $\sigma$ 进行估计的经验公式：

于是得到 Silverman 方法的最终带宽公式：

若选择 Gauss 核，则 $C_\mathrm{G} \approx 1.06$（这也是著名的经验带宽公式）；若选择 Epanechnikov 核，则 $C_\mathrm{E} \approx 2.34$；若选择均匀核，则 $C_{\mathrm{U}} \approx 1/36$. 具体计算过程略.

4.3 插入法（Plug-in）

$\color{blue}{\textbf{Def 4.6 } \boldsymbol{k} \textbf{ 阶核函数: }}$ 设 $K$ 是一个核函数，$k \in \mathbb{N}_+$. 若对任意的 $r = 1, 2, \dots, k - 1$ 满足 $\int_{\mathbb{R}} u^rK(u) \mathrm{d}u = 0$，且 $\int_{\mathbb{R}} u^kK(u) \mathrm{d}u > 0$，则称 $K$ 是一个 $k$ 阶核函数.

$\color{blue}{\textbf{Prop 4.7: }}$ $K$ 是 $k$ 阶核函数，则 $k$ 是偶数.

$\color{brown}{\textbf{Proof: }}$ 只需注意到 $r < k$ 是奇数时，$u^r K(u)$ 是奇函数即可. $\square$

下面，我们介绍一个重要的推导思路，引理 $4.8$ 是其最简单的形式. 后续我们将利用其思想推导各等式，不作单独证明.

$\color{blue}{\textbf{Lem 4.8: }}$ $s \in \mathbb{N}_+$，$f \in C^{2s}(\mathbb{R})$，$\lim\limits_{x \rightarrow \infty} \nabla^r f(x) = 0, \space \forall r = 1, 2, \dots, 2s$，则：

$\color{brown}{\textbf{Proof: }}$

记：

基于引理 $4.8$，我们可以将对 $R(\nabla^2 f)$ 的计算转化为对 $\psi_4$ 的计算. 下面我们给出 $\psi_r$ 的一种非参数估计方式：

其中 $g$ 是带宽，$L$ 是核函数. 注意这里的带宽和核函数可以与 KDE 的带宽与核函数不同. 于是问题转化为求解 $g$ 使得 $\mathrm{AMSE} \left[ \hat{\psi}_r(g) \right]$ 最小. 其中：

利用和前文一样的思想，先给出 $\mathbb{E} \left[ \hat{\psi}_r(g) \right]$：

设 $L$ 是 $k$ 阶核函数，其中：

注意到 $\frac{n - 1}{n} \sim 1$，得到：

故偏差为：

接下来，我们推导 $\mathrm{Var} \left[ \hat{\psi}_r(g) \right]$，记 $\nabla^r L_g(x^{(i)} - x^{(j)}) = L_{ij}$，对 $\mathrm{Var} \left[ \hat{\psi}_r(g) \right]$ 进行展开：

对 $\mathrm{Cov}(L_{ij}, L_{kl})$ 进行分类讨论：

1. 若 $i = j$ 或 $k = l$，此时 $\mathrm{Cov}(L_{ij}, L_{kl})$ 退化为 $\mathrm{Cov}\left( L_g(0), L_g(0) \right) = 0$.
2. 若 $i \not= j$ 且 $k \not= l$，则：
   1. 若 $\vert \{i, j\} \cap \{k, l\} \vert = 2$，这样的情况有 $2 A_n^2 = 2n(n - 1)$ 种，根据核函数的对称性，有：$\mathrm{Cov}(L_{ij}, L_{kl}) = \mathrm{Var}(L_{ij}, L_{ij})$.
   2. 若 $\vert \{i, j\} \cap \{k, l\} \vert = 1$，不妨设 $j = l$，这样的情况有 $4A_n^3 = 4n(n - 1)(n - 2)$ 种，根据核函数的对称性，有：$\mathrm{Cov}(L_{ij}, L_{kl}) = \mathrm{Var}(L_{ij}, L_{jk})$.
   3. 若 $\vert \{i, j\} \cap \{k, l\} \vert = 0$，$L_{ij}$ 与 $L_{jk}$ 独立，此时有 $\mathrm{Cov}(L_{ij}, L_{kl}) = 0$.

于是，得到 $\mathrm{Var} \left[ \hat{\psi}_r(g) \right]$ 的展开式：

其中：

最后，注意到 $\mathbb{E}(L_{1, 2}) = \psi_r + o(1)$，有：

从而得到渐进均方误差：

这一式子十分复杂，难以通过求偏导并令导数为 $0$ 的方法求解理论最优带宽，但是我们可以计算其近似最优带宽：通过分析发现，当 $g$ 增大时，$g^{2k}$ 成为主导项，$\mathrm{AMSE} \left[ \hat{\psi}_r(g) \right]$ 显著增高；$g$ 减小时，$\frac{1}{g^{2r + 2}}$ 成主导项，$\mathrm{AMSE} \left[ \hat{\psi}_r(g) \right]$ 也会显著增高. 因此，减小 $\mathrm{AMSE} \left[ \hat{\psi}_r(g) \right]$ 的关键在于选取一个平衡的 $g$，使得 $\mathrm{AMSE} \left[ \hat{\psi}_r(g) \right]$ 既不会太大也不会太小. 我们根据这样的思想，忽略不占据主导位置的 $\frac{2}{n^2 g^{2r + 1}} \psi_0 R(\nabla^r L)$ 项和不含 $g$ 的 $\frac{4}{n} \left[ \int_{\mathbb{R}} \left[ \nabla^r f(x) \right]^2 f(x) \mathrm{d}x - \psi_r^2 \right]$ 项，令：

得到：

将 $g_{\mathrm{approx-opt}}$ 带入到 $(51)$ 中，仅保留 $n$ 最高次数项（主导项，其他项省略）：

一般来讲，我们使用的核函数是 $2$ 阶核函数，以下也仅以 $k = 2$ 为例进行讨论. 带入 $k = 2$，得：

我们先回到问题本身，我们的目标是计算 KDE 的最优带宽，即：

将 $\psi_4$ 替换为核估计量 $\hat{\psi}_4(g)$，便得到插入法的最优带宽计算公式：

然而，$\hat{\psi}_4(g)$ 的计算依赖于 $g_{\mathrm{approx-opt}, 4}$，根据公式，$g_{\mathrm{approx-opt}, 4}$ 的计算又依赖于 $\psi_6$. 依次类推，最终会导致无限的循环：对 $\psi_r$ 的估计总是依赖于 $\psi_{r + 2}$. 解决这一问题的方法为：选取一个上限 $\mathscr{l} \in \mathbb{N}_+$，在计算到 $\psi_{2\mathscr{l} + 4}$（或者更一般地，计算到 $\psi_{k\mathscr{l} + 4}$）时停止，转而用另一种方法进行估计（推导过程将放在本节末尾处）：

例如，我们选取 $\mathscr{l} = 2$，先估计 $\hat{\psi}_8 = \hat{\psi}_{\mathrm{NS}, 8}$，利用 $\hat{\psi}_8$ 计算 $g_{\mathrm{approx-opt}, 6}$，进而得到 $\hat{\psi}_6$. 利用 $\hat{\psi}_6$ 我们可以计算 $g_{\mathrm{approx-opt}, 4}$，进而得到 $\hat{\psi}_4$，最终得到 $h_{\mathrm{pluge-in}, 2}$. 我们将上述方法综合记为（在 $k = 2$ 意义下）：

以上方法称为 $\mathscr{l}$ 阶段插入法. 目前尚未存在一个客观选择 $\mathscr{l}$ 的标准方法，一般而言，取 $\mathscr{l} = 2$ 即可.

现在来推导式子 $(58)$，假设 $f$ 服从正态分布 $\mathcal{N}(\mu, \sigma^2)$，则：

其中 $\mathrm{He}_r (z)$ 是概率论的埃尔米特多项式（Hermite Polynomial, HP），其相关知识这里不作叙述.

记标准正态分布函数为 $\Phi$，对式子 $(64)$ 进行变量代换，得：

利用 HP 的生成函数：

一方面：

另一方面：

比较式子 $(67)$ 和 $(68)$ 可得：

于是当 $r$ 是偶数时：

4.4 留一交叉验证法（Leave-One-Out CV, LOO-CV）

留一交叉验证法同样使用 $\mathrm{MISE}$ 作为最优带宽的选择标准. 不同的是，留一交叉验证法构建了 $\mathrm{LOOCV}$ 函数，

首先，回顾 $\mathrm{MISE}$ 的公式：

其中 $\int_{\mathbb{R}} f^2(x) \mathrm{d}x$ 与 $h$ 无关. 因此最小化 $\mathrm{MISE}$ 等价于最小化 $\mathbb{E}\left[ \int_{\mathbb{R}} \hat{f}^2 (x) \mathrm{d}x \right] - 2 \mathbb{E}\left[ \int_{\mathbb{R}} \hat{f}(x) f(x) \mathrm{d}x \right]$. 由于 $f(x)$ 未知，无法直接求解. 下面采用 Monte-Carlo 算法的思想，对其进行近似.

对于每个 $j = 1, 2, \dots, n$，作训练集 $\mathcal{S}_{-j} = \mathcal{S} \setminus \{x^{(j)}\}$，利用 $\mathcal{S}_{-j}$ 作函数 $\hat{f}_{h, -j}(x)$：

显然有：

定义 $\mathrm{LOOCV}$ 函数：

容易验证：

因此最小化 $\mathrm{LOOCV}(h)$ 等价于最小化 $\mathrm{MISE}(\hat{f}(x))$. 即：

下面补充 $\mathrm{LOOCV}$ 中 $\int_{\mathbb{R}} \hat{f}^2 (x) \mathrm{d}x$ 的计算方法：

例如对于 Gauss 核，可直接利用其卷积公式 $(K_{\mathrm{G}} * K_{\mathrm{G}} )(t) = \frac{1}{2\sqrt{\pi}}e^{-t^2/4}$ 快速计算.

4.5 有偏交叉验证法（Biased CV, BCV）

有偏交叉验证的核心思想是用 $R(\nabla^2 \hat{f})$ 来估计 $R(\nabla^2 f)$. 事实上，前者是后者的有偏估计，我们先计算 $R(\nabla^2 \hat{f})$ 的期望：

1. 若 $i = j$，$\mathbb{E} \left[ \left( \nabla^2 K_h \ast \nabla^2K_h \right) \left( x^{(i)} - x^{(j)} \right) \right] = \mathbb{E} \left[ h^{-5}R(\nabla^2 K) \right] = h^{-5}R(\nabla^2 K)$；

2. 若 $i \not= j$，有：

   $$\begin{align} \mathbb{E} \left[ \left( \nabla^2 K_h \ast \nabla^2K_h \right) \left( x^{(i)} - x^{(j)} \right) \right] &= \iint_{\mathbb{R^2}} \left( \nabla^2 K_h \ast \nabla^2K_h \right) \left( s - t \right) f(s) f(t) \mathrm{d}s \mathrm{d}t \nonumber \\ &= \int_{\mathbb{R}} \left( \nabla^2 K_h \ast \nabla^2K_h \right) \left( u \right) \left( f \ast f \right) \left( u \right) \mathrm{d}u. \end{align}$$

   带入，得：

   $$\begin{align} \mathbb{E} \left[ R(\nabla^2 \hat{f}) \right] &= \frac{1}{nh^5} R(\nabla^2 K) + \frac{n - 1}{n} \int_{\mathbb{R}} \left( \nabla^2 K_h \ast \nabla^2K_h \right) \left( u \right) \left( f \ast f \right) \left( u \right) \mathrm{d}u . \end{align}$$

   下面证明 $\lim\limits_{h \rightarrow 0^+} \int_{\mathbb{R}} \left( \nabla^2 K_h \ast \nabla^2K_h \right) \left( u \right) \left( f \ast f \right) \left( u \right) \mathrm{d}u = R(\nabla^2 f)$，根据 Parseval 定理及卷积定理，有：

   $$\begin{align} &\quad \space \int_{\mathbb{R}} \left( \nabla^2 K_h \ast \nabla^2K_h \right) \left( u \right) \left( f \ast f \right) \left( u \right) \mathrm{d}u \nonumber \\ &= \int_{\mathbb{R}} \mathcal{F} \left[ \nabla^2 K_h \ast \nabla^2K_h \right] (t) \cdot \mathcal{F} \left[ f \ast f \right] (t) \mathrm{d}t \nonumber \\ &= \int_{\mathbb{R}} \left[ \mathcal{F} \left[ \nabla^2 K_h\right] (t) \right]^2 \cdot \left[ \mathcal{F} \left[ f \right] (t) \right]^2 \mathrm{d}t \nonumber \\ &= \int_{\mathbb{R}} t^4\left[ \mathcal{F} \left[ K \right] (ht) \right]^2 \cdot \left[ \mathcal{F} \left[ f \right] (t) \right]^2 \mathrm{d}t. \end{align}$$

   由 $\mathcal{F} [K] (ht) \rightarrow \mathcal{F} [K](0) = 1$ 得：

   $$\begin{align} &\quad \space \int_{\mathbb{R}} \left( \nabla^2 K_h \ast \nabla^2K_h \right) \left( u \right) \left( f \ast f \right) \left( u \right) \mathrm{d}u \nonumber \\ &= \int_{\mathbb{R}} t^4 \left[ \mathcal{F} \left[ f \right] (t) \right]^2 \mathrm{d}t \nonumber \\ &= \int_{\mathbb{R}} \left[ \mathcal{F} \left[ \nabla^2 f \right] (t) \right]^2 \mathrm{d}t \nonumber \\ &= \int_{\mathbb{R}} \left[ \nabla^2 f(t) \right]^2 \mathrm{d}t \nonumber \\ &= R(\nabla^2 f). \end{align}$$

   于是：

   $$\begin{align} \mathbb{E} \left[ R(\nabla^2 \hat{f}) \right] &= \frac{1}{nh^5} R(\nabla^2 K) + R(\nabla^2 f) + o(\frac{1}{n}). \end{align}$$

   因此，$R(\nabla^2 f)$ 的一个良好的无偏估计为：

   $$\begin{equation} \widetilde{R(\nabla^2 f)} = R(\nabla^2 \hat{f}) - \frac{1}{nh^5} R(\nabla^2 K). \end{equation}$$

   类似 Silverman 方法，给出有偏交叉验证函数：

   $$\begin{align} \mathrm{BCV}(h) &= \frac{1}{4}h^4 \mu_2^2(K) \widetilde{R(\nabla^2 f)} + \frac{R(K)}{nh}, \\ h_{\mathrm{BCV}} &= \arg \mathop{\min}\limits_{h} \mathrm{BCV}\left( h \right). \end{align}$$

   容易得出：

   $$\begin{equation} h_{\mathrm{BCV}} = \left(\frac{R(K)}{\mu_2^2(K) [R(\nabla^2 \hat{f}) - \frac{1}{nh^5} R(\nabla^2 K)] n} \right)^{1/5} \end{equation}$$

4.6 直观法

直观法是核密度估计中带宽选择的 “经验驱动型方法”，核心思路是通过可视化密度曲线形态和主观判断平滑性与细节的平衡选择带宽，无需复杂的数学推导或迭代计算，适用于对估计精度要求不高、数据结构简单或需快速探索的场景.

直观法一般给出一系列可选择的带宽 $\mathcal{H} = \{h_k \mid k = 1, 2, \dots, m\}$，生成多组密度函数，进而选取最优带宽（此处的最优指最适合展示或计算的带宽，非理论最优带宽）. 例如，在处理小样本问题过程中，Silverman 方法和 LOO-CV 法都有较大误差. 此时可以使用直观法，确定较为合适的带宽. 此外，在数据分析初期，也可通过直观法生成多组带宽的密度曲线，辅助分析数据分布趋势.

5. 核函数的选择

5.1 Epanechnikov 核：理论最优核

在讨论最优核之前，我们先对核函数补充两个限制：

1. 要求任意核函数 $K$ 必须有紧支集，即 $\mathrm{supp}(K) = \overline{K^{-1}(\mathbb{R} \setminus \{0\})}$ 有界，不妨设 $\mathrm{supp}(K) = [-1, 1]$.
2. 核函数连续.

不符合上述两个限制的核包含 Gauss 核及均匀核，我们将在后面进行讨论. 而在选择最优带宽下，满足上述两个限制的最优核函数是 Epanechnikov 核. 我们将继续在最小化 $\mathrm{AMISE}$ 下讨论.

在前文中我们计算了 $\mathrm{AMISE} \left[ \hat{f}(x) \right]$ 的最小值：

这表明在最优化带宽下最小化 $\mathrm{AMISE}$ 等价于最小化 $R^2(K) \mu_2(K)$. 接下来我们使用变分法求解.

由于核函数 $K$ 具有紧支集，故对 $K$ 的各种积分的界限都可以改为 $[-1, 1]$，仍记为原符号. 设 $\mu_2(K) = A$，在 $\mu_2(K) = A$ 及 $\int_{[-1, 1]} K(u) \mathrm{d}u$ 的约束下，最小化 $R(K)$，可以建立 Lagrange 函数：

对 $\mathcal{L}$ 关于 $K$ 求变分：

令其泛函导数为 $0$，解得：

考虑到约束 $K(-1) = K(1) = 0$，得：

再考虑到 $K$ 的归一性，有：

上式中 $K(u)$ 即为 Epanechnikov 核.

5.2 Gauss 核：实践最优核

Gauss 核因适配多数实践场景需求，被广泛视为‘实践最优核”.

首先我们讨论 Gauss 核的有效区域. 虽然 Gauss 核不具有紧支集，但是对于 $h > 0$：

1. 当 $\vert x \vert > 3h$ 时，$\exp\left(-\frac{x^2}{2h^2}\right) < \exp\left(-\frac{9}{2}\right) \approx 0.011$，贡献已不足 1%；
2. 当 $\vert x \vert > 5h$ 时，$\exp\left(-\frac{x^2}{2h^2}\right) < \exp\left(-\frac{25}{2}\right) \approx 3 \times 10^{-6}$，贡献几乎可忽略.

即：Gauss 核的实际有效区域是 $[-5h, 5h]$，在实际计算中完全可将有效区域外样本的权重视为 $0$，计算效率与紧支集核几乎无差异. Gauss 核的广泛应用还离不开以下三个个关键特质：

1. 无线光滑性：$K_{\mathrm{G}} \in C^{\infty}(\mathbb{R})$，其生成的密度函数 $\hat{f}$ 也继承了这一优良数学性质，能完美匹配真实密度的光滑特性. 而其它核（如：Epanechnikov 核）在需要高阶导数条件下（如：峰值分析、梯度分析）表现不佳.
2. 抗干扰性：由于 Gauss 核的权重是连续平滑的，且对远离待估计点的 “异常值样本” 赋予极低的权重，因此在数据存在少量噪声或异常值时，Gauss 核的密度估计结果不易被干扰. 而其它核若恰好将异常值纳入支集内，会直接影响局部密度估计，抗干扰性相对较弱.
3. 高维稳定性：高维数据中，“距离” 的定义会因维度灾难导致样本分布稀疏，紧支集核可能因支集内样本过少而无法有效估计密度（甚至支集内无样本，导致估计值为 $0$）. 而 Gauss 核的指数衰减特性能通过带宽调整，灵活适配高维数据的稀疏性——即使局部样本少，也能通过平滑的权重分配维持密度估计的连续性，避免出现 “零密度区域” 的不合理结果.

此外，在其它分析中，Gauss 核还具有严格正定性等优良性质，是核分析中的常用函数.

5.3 均匀核

均匀核既非理论最优核，也没有 Gauss 核的优良数学性质，反而缺点重重. 但是均匀核在特定情况下有其作用. 例如，在认为“某点周围 $K$ 个样本的类别完全同等重要”时，可以使用均匀核. 此外，在一些低计算成本场景下（如：嵌入式开发、实时数据处理），均匀核因其计算简单而受青睐.

总而言之，不同核有其各自的优缺点，实际应用中要根据应用场景的特点，选择合适的核函数.

---

更新日志

[2025.09.03]：

1. 调整了学习顺序，优化叙述逻辑；
2. 添加“插入法”、“有偏交叉验证法”.