核密度估计（KDE）（二）

核密度估计（KDE）（二）

1. 自适应带宽选择

在前文中，KDE 算法要求各样本点共用同一带宽，即认为各样本点的影响大小、影响范围是一样的. 但在实际生活中，这一假设并不总是成立.

以某一大型区域的便利店服务研究为例：在城市中，人口密集、各类设施丰富，便利店分布较为密集. 但每家便利店主要服务对象仅为周边一小片区域的居民，影响范围小；而在乡村中，人口相对分散，便利店数量较少，一家便利店往往要服务整个村子甚至邻近几个村落的居民，服务范围就大得多. 在这种情况下，位于高密度区域的样本点的影响范围应该更小；位于低密度区域的样本点的影响范围应该更大. 于是我们给出如下自适应修正模型：

对于某一初始带宽 $h$，计算其密度估计函数：

在此基础上，引入带宽修正公式：

其中 $F = \left[ \prod\limits_{i = 1}^{n} \hat{f}_h \left( x^{(i)} \right) \right]^{1 / n}$，最终得到修正后的密度估计函数：

对于其它问题，可运用此类思想对带宽进行自适应修正.

2. 多元 KDE 算法

2.1 多元 KDE 算法引入

基于一元 KDE 的基本原理，我们可以将 KDE 推广到 $d \space (d \in \mathbb{N}_+)$ 元情况.

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

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

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

核函数	$d$ 元核公式
$d$ 元高斯核	$K(\boldsymbol{u}) = \frac{1}{(2 \pi)^{d / 2}} \exp(-\frac{1}{2} \boldsymbol{u}^\mathrm{T}\boldsymbol{u})$
$d$ 元 Epanechnikov 核	$K(\boldsymbol{u}) = \begin{cases} \frac{d + 2}{2 V_d}(1 - \boldsymbol{u}^\mathrm{T}\boldsymbol{u}), & \text{if } \vert \boldsymbol{u} \vert < 1 \\ 0, & \text{otherwise}\end{cases}$
$d$ 元均匀核	$K(\boldsymbol{u}) = \begin{cases} \frac{1}{V_d}, & \text{if } \vert \boldsymbol{u} \vert < 1 \\ 0, & \text{otherwise}\end{cases}$

其中 $V_d = \frac{\pi^{d/2}}{\Gamma(\frac{d}{2} + 1)}$ 是 $d$ 维单位球体积.

根据一元 KDE 的基本思想，我们很容易得到下述公式：

其中 $\boldsymbol{H} = h^2 \boldsymbol{I}_d$，$\boldsymbol{I}_d$ 是 $d$ 阶单位矩阵，$K(\boldsymbol{u}) = \frac{1}{V_d} \mathbb{1} \left( \parallel \boldsymbol{u} \parallel_2 \leq 1 \right)$ 是 $d$ 元均匀核. 容易验证：

若将上式中的 $K(x)$ 替换为其它核函数，得到的 $\hat{p}(\boldsymbol{x})$ 仍为概率密度函数. 下面我们将对 $\boldsymbol{H}$ 进行进一步的讨论：

首先，样本的不同维度的数据规模不同，适用的带宽也不同，因此，我们可以对 $\boldsymbol{H}$ 进行扩展：

这样的 $\boldsymbol{H}$ 是十分容易参与计算的，其几何意义也十分明显：样本点的各维度彼此独立，各自享用自己的带宽. 如果各维度并不彼此独立，则可考虑使用样本的协方差矩阵构建带宽矩阵：

其中参数 $\varepsilon > 0$ 起微小扰动作用，避免协方差矩阵 $\boldsymbol{\Sigma}$ 半正定.

2.2 多元 KDE 算法

$\color{blue}{\textbf{Algorithm: Multivariate Kernel Density Estimation}}$ 设 $\mathcal{S} = (\boldsymbol{x}^{(1)}, \boldsymbol{x}^{(2)}, ..., \boldsymbol{x}^{(n)})$ 是从总体中抽取的一组独立同分布（i.i.d.）样本数据，$\boldsymbol{x}^{(i)} \in \mathbb{R}^d$，$K(\boldsymbol{x})$ 是给定的 $d$ 元核函数，$\boldsymbol{H} \in \mathrm{SPD}(\mathbb{R}, d)$ 为给定（$d$ 阶实正定）带宽矩阵，则核密度估计给出的概率密度估计函数为：

在这里，我们继续引入积核（Product Kernel）的概念：

$\color{blue}{\textbf{Def 2.2 积核: }}$ 设 $K_1, K_2, \dots, K_d$ 是 $d$ 个一元核函数，定义多元核函数 $K: \mathbb{R}^d \to [0, +\infty)$ 满足：

其中 $\boldsymbol{u} = (u_1, u_2, \dots, u_d)^\mathrm{T} \in \mathbb{R}^d$，则称 $K$ 为基于 $K_1, K_2, \dots, K_d$ 的积核.

显然，$d$ 元 Gauss 核是 $d$ 个一元 Gauss 核函数的积核. 使用积核的一个好处是其支持并行计算；此外，在需要调整带宽、增加估计点时，积核方法的计算复用性显著高于非积核方法.

我们接下来讨论多元 KDE 算法的渐进性质，将原来的三条假设推广至多元：

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

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

3. 记 $\boldsymbol{H} = \boldsymbol{H}_n$，假设 $n \rightarrow \infty$ 且 $\boldsymbol{H} \rightarrow 0^+$ 时满足 $n\det \boldsymbol{H}^{1/2} \rightarrow +\infty$.

利用多元 Taylor 定理及一元 KDE 渐进性质的分析思路，我们可以得到（不做证明）：

其中 $\mathrm{vec}(\cdot)$ 表示将目标 $n \times m$ 矩阵转换为 $nm$ 阶向量，$\mathrm{H}(\cdot)$ 表示 Hessian 矩阵.

2.3 带宽选择

带宽选择的基本方法已经在前文叙述，以下仅给出多元情况下的带宽选择公式. 首先我们推广 $\mathrm{MSE}$ 和 $\mathrm{MISE}$：

它们的渐进形式为：

其中用到了一个定理：

$\color{blue}{\textbf{Thm 2.4: }}$ $A, B \in \mathbb{R}^{n \times n}$，$\mathrm{vec}(A)^{\mathrm{T}} \mathrm{vec}(B) = \mathrm{tr}\left( A^{\mathrm{T}} B \right)$.

$\color{brown}{\textbf{Proof: }}$ 设 $A = \left[ a_{ij} \right]$，$B = \left[ b_{ij} \right]$，$C = A^\mathrm{T}B = \left[ c_{ij} \right]$有：

于是：

于是，渐进意义下的最优带宽为：

与一元 KDE 不同，通过上式无法求出 $\boldsymbol{H}_{\mathrm{asym-opt}}$ 的表达式. 但是特别地，若 $\boldsymbol{H} = h^2 \boldsymbol{I}_d$，则：

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

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

尽管 $(34)$ 依赖于大量化简，但其揭示了一个重要的规律：当 $d$ 增大时，最优带宽 $h_{\mathrm{asym-opt}}$ 也会增大. 此外，我们发现 $(34)$ 的计算也面临着和一元 KDE 同样的问题：起计算依赖于关于真实密度函数 $f$ 的泛函. 我们首先采取和 Silverman 方法一致的思想，用正态分布函数来近似 $f$，得到（不作具体计算）：

在前文中，我们还提到了插入法. 对于多元 KDE，其插入法推导思路与一元 KDE 无异，但各阶段最优带宽 $\boldsymbol{G}$ 没有明确的公式，不适合实际计算，这里不再展开. 接下来我们讨论交叉验证法，交叉验证法可以很好地适应多元情况.

留一交叉验证法可以表述为：

有偏交叉验证法可以表述为：

2.4 多元 KDE 算法下的导数估计

设 $f: \mathbb{R}^d \rightarrow \mathbb{R}$ 是通过 KDE 算法生成的概率密度函数，我们想讨论其导数. 对于一阶和二阶导数，我们分别有：

需要指出的是，$(29)$ 式是矩阵形式，显然不利于推广到高维，我们利用 Kronecker 积将其转换为向量形式：

我们将上述符号记为 $\nabla^{\otimes 2} f(\boldsymbol{x})$. 一般地，我们给出 $f(\boldsymbol{x})$ 的 $r$ 阶导数：

$\color{blue}{\textbf{Lem 2.3 多元 Taylor 定理: }}$ 设 $f: \mathbb{R}^p \to \mathbb{R}$，$\boldsymbol{x} \in \mathbb{R}^p$。若 $f \in C^r[B(\boldsymbol{x}, \delta)]$，则对任意 $\vert \boldsymbol{h} \vert < \delta$，有：

证明略. 下面给出 $f$ 的导数估计量：

其中：

特别地，取 $r = 1$ 时，有：

若定义缩放一阶梯度核函数：

则一阶导数估计量可写作：

对于 $r > 1$ 的讨论我们不多赘述，可根据实际情况进行讨论. 注意到 $\boldsymbol{H}^{-1/2}$ 的存在，这说明了 $(37)$ 的带宽选择不同于 $(8)$ 的带宽选择，我们给出 $r$ 阶导数下 $\boldsymbol{H}_{\mathrm{NS}, r}$ 的表达式（对 $\mathrm{AMSE}$ 的推广、Silverman 方法的应用不作叙述）：

当 $r = 0$ 时，$\boldsymbol{H}_{\mathrm{NS}, 0}$ 退化为 $\boldsymbol{H}_{\mathrm{NS}}$. 随着 $r$ 的增大，$\boldsymbol{H}_{\mathrm{NS},r}$ 也会增大，这表明若使用密度估计的最优带宽来估计导数的最优带宽，会导致导数估计出现“欠平滑”现象.

3. 线特征的核密度估计（KDE-Linear, KDE-L）

线特征（如道路、河流、轨迹、等高线等）是地理信息科学（Geographic Information Science, GIS）、计算机视觉、模式识别等领域中常见的空间数据类型. 与点特征的核密度估计（KDE）不同，线特征的核密度估计需考虑线的连续性、长度、方向等几何属性，核心是通过“核函数”将线特征的空间分布“平滑化”，最终生成反映线特征空间集聚程度的密度表面，用于识别线特征的密集区域（如城市路网密集区）或稀疏区域.

在这种意义下，我们不再需要保障最终的密度估计函数 $\hat{f}$ 满足 $\int_{\mathbb{R}^d} \hat{f}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} = 1$，而仅使用 $\hat{f}(\boldsymbol{x})$ 的大小来判断 $\boldsymbol{x}$ 处的密度高低.

设 $l$ 是一个线特征，为确定 $l$ 的概率影响，我们可以采取采样点的方法. 作划分 $T: t_0, t_1, t_2, \dots, t_n$，其中 $t_0$ 是 $l$ 的起点，$t_n$ 是 $l$ 的终点，$t_0, t_1, \dots, t_n$ 依次位于 $l$ 上，$\overline{t_{i - 1}, t_i}$ 表示 $t_{i - 1}$ 和 $t_i$ 点之间的线段，$\parallel T \parallel = \max\limits_{i = 1, 2, \dots, n} \vert \overline{t_{i - 1}, t_i} \vert$. 任取 $s_i \in \overline{t_{i - 1}, t_i}$，考虑使用 $\mathcal{S} = \{ s_1, s_2, \dots, s_n \}$ 来估计 $l$ 的概率影响，对于任意点 $x$：

其中 $K(\cdot)$ 是选定核函数. 随着划分 $T$ 越来越精细，$\hat{f}_T(x)$ 越接近真实的密度影响：

忽略掉系数 $(\vert \overline{t_0, t_n} \vert)^{-1}$，得到：

值得指出的是，忽略掉系数 $(\vert \overline{t_0, t_n} \vert)^{-1}$ 是有意义的：一方面，由于不需要进行归一化，忽略该系数不会产生错误；另一方面，由于不同线特征的长度可能不同，若保留该系数，则较长的线特征产生的影响更广，导致其影响强度降低，使得不同线特征之间影响不匹配.

对于带宽 $h$ 的选择，经验法（即：前文所提的直观法）是尤为重要的. 这是因为此类问题通常关注可视化的结果，得到如道路路网密度图等图像，因此视觉上的最优性是选取带宽的第一要素；对于部分问题，还要兼顾其计算复杂性和精度. 此外，自适应的带宽选择方法也会被使用.

4. 网络核密度估计（Network KDE, NKDE）

4.1 网络密度核估计算法

网络核密度估计是传统核密度估计在网络空间中的拓展，核心是解决“空间点事件在网络结构（如道路、河流、管线等线性网络）上的密度分布估计”问题，而非传统 KDE 针对的二维/三维欧氏空间.

在现实场景中，许多事件（如交通事故、店铺分布、犯罪案件）的发生和扩散严格依赖于网络（例如交通事故仅会沿道路分布，而非在道路外的区域），传统 KDE 直接用欧氏距离计算密度会忽略网络拓扑结构（如道路的连通性、弯曲度），导致估计结果失真；而 NKDE 通过将 “欧氏距离” 替换为 “网络距离”，并结合网络拓扑特性调整核函数，实现对网络上点事件密度的精准刻画. 同样地，网络和密度估计也不需要进行归一化.

设 $G = (V, E)$ 是有向连通图，其中 $V$ 是顶点集，$E$ 是边集. 我们用 $d_+(v)$ 表示顶点 $v$ 的出度，$\mathrm{out}(v)$ 表示离开顶点 $v$ 的边的集合. 显然满足：$\vert \mathrm{out}(v) \vert = d_+(v)$.

下面，我们来定义网络：记 $\mathcal{N} = (G, \omega)$ 称为一个网络，其中 $G = (V, E)$ 是有向连通图，$\omega$ 是定义在 $E$ 上的函数，用于指明 $E$ 中边 $e$ 的权重（对应现实意义中的：道路权重、河流流量、管道大小等）. 称 $x \in \bigcup E$ 为 $\mathcal{N}$ 中的点，对于 $\mathcal{N}$ 中的任一点 $x_1, x_2$，记 $x_1$ 到 $x_2$ 的最短路径长度为 $\mathrm{dm}(x_1, x_2)$，经过的顶点的全体构成的集合为 $\mathrm{Vm}(x_1, x_2)$.

设 $\mathcal{S} = (x^{(1)}, x^{(2)}, ..., x^{(n)})$ 是 $\mathcal{N}$ 中的样本，$K$ 是取定一元核函数. 对于任意样本点 $x^{(i)}$，不妨设 $x^{(i)} \in e_i$，对于 $e_i$ 上任意点 $x$，$x^{(i)}$ 在 $x$ 处的影响为：

此时的影响衰减函数即为一元 KDE 的影响衰减函数. 下面考虑 $x^{(i)}$ 对位于其它边上的点的影响：设 $e_j$ 的起点通过顶点 $v$ 连接到 $e_i$，$x \in e_j$，则可以给出 $x^{(i)}$ 在 $x$ 处的影响：

其直观的实际意义是：$x^{(j)}$ 的影响随着 $x$ 移动到 $e_j$ 的终点 $v$ 处而发生了分裂，平均分为 $d_+$ 份，沿着 $\mathrm{out}(v)$ 中的每条边继续延申. 如果考虑到 $\omega$ 的作用（即：影响的分裂不是均分，而是依权重分裂），则有：

据此，对于一般 $\mathcal{N}$ 中的点 $x$，$x^{(i)}$ 对 $x$ 的影响为：

（说明：上式实际表示的是 $x^{(i)}$ 对其下游点的影响，若要表示对上游点的影响，只需将 $\mathrm{Vm}(x^{(i)}, x)$ 改为 $\mathrm{Vm}(x, x^{(i)})$，下同）

其中 $e_v$ 表示从 $x^{(i)}$ 到 $x$ 的最短路径中以 $v \space \left(v \in \mathrm{\mathrm{Vm}(x^{(i)}, x)} \right)$ 为起点的边. 尽管这一式子十分复杂，在实际应用中（无论是选取有紧支集的核函数还是 Gauss 核函数），只需考虑带宽范围内的点即可.

最后，基于全部的样本，我们可以给出 $\mathcal{N}$ 上的样本密度估计函数：

4.2 基于线特征的网络核密度估计算法

在某些应用网络密度核估计算法的场景中，样本点不再是点特征，而是线特征，如：出租车的载客点分布（说明：实际应用中，无法获取出租车载客点的实际位置，而是通过连续采样，粗略地认为载客点位于某两采样点间. 而这两个采样点间的道路即为一个线特征）. 下面我们将 NKDE 推广至线特征：

先给出线性事件的定义：记 $\mathscr{l} = ( s, t, p )$，其中 $s \in e_s, \space t \in e_t, \space \exist e_s, e_p \in E$，$p$ 是线性事件的其它额外信息（如起止时间、ID 编号等），不在数学讨论范围中.

设 $\mathcal{S_\mathrm{linear}} = (\mathscr{l}^{(1)}, \mathscr{l}^{(2)}, ..., \mathscr{l}^{(n)})$ 是网络 $\mathcal{N}$ 中的线特征样本，每个特征的原始影响衰减函数为：

联合上述分析，得到受网络约束的影响衰减函数：

最后，基于全部的样本，我们可以给出 $\mathcal{N}$ 上的线特征样本密度估计函数：

最后，我们指出本节所述仅是一种推广的思想，对于实际问题需要进行具体分析，调整式子，从而给出更优的解答. 例如，在前文所述出租车载客点分布的问题中，我们有时需要分别讨论道路两侧的载客点. 这就需要我们在求平均值时只考虑“同向”线性事件的影响，而忽略“反向”线性事件的影响.