概率密度函数简称密度函数或密度，是概率统计中的重要概念，用于描述连续型随机变量所服从的概率分布。统计学中许多应用都涉及对密度函数的估计，即根据样本观测数据来估计总体的密度函数。了解随机变量的概率分布可以帮助我们更好地把握数据的特征，有助于进一步的统计分析，如条件分布的估计、非参数回归、可靠性分析等。同时，概率密度函数还可用于图像分类、图像处理^[1]、风险分析^[2-3]、聚类^[4]等相关领域。从某种意义上来说，有了密度函数就获得了随机变量的全部信息。

密度函数的估计通常分为参数密度估计和非参数密度估计。核密度估计(kernel density estimation, KDE)^[5-6]作为一种非参数密度估计，完全从数据出发，利用核函数的叠加平均来近似真实的密度函数。KDE需要恰当地选择核函数和窗宽参数，且窗宽参数的选择对于估计结果至关重要，尤其在高维情形下更是如此。拇指法则^[7-8]是一种简单的最优窗宽选择方法，但在数据集不满足正态分布的假定时会导致结果不准确。作为一种参数化方法，重构方法^[9]是基于插值的函数近似方法，其思想来源于插值器快速的收敛速度，对于充分光滑的函数，甚至可以以指数阶的速度收敛，这种收敛速度远远高于中心极限定理中统计量的收敛速度，因此函数本身与其插值器之间的误差远小于统计误差。对于回归函数估计，函数重构法通过对有限个节点插值近似从而重构出整个函数。

本文关注函数重构法在密度估计中的应用，提出一种重构法中节点的序贯选择方法。该方法操作简单，可以有效缓解由于节点选择不当而造成的估计效果波动，同时根据参数的实际意义，给出密度函数的区间估计的构造方法。该方法可以很自然地应用到其他无监督学习问题中。模拟结果表明，利用序贯方法选择节点构造的密度估计效果优于以往的方法。

1 基于序贯方法的重构密度估计 1.1 重构密度估计

重构方法通过构造一个选定节点的插值器来参数化未知函数，节点处的函数值便是要估计的参数，因此重构方法的参数有着明确解释，这为重构方法的应用提供了便利。假定f为定义在D⊂ $\mathbb{R}^d$上的光滑函数，$\mathcal{X}=\left\{\left(\boldsymbol{x}_1, \boldsymbol{y}_1\right), \cdots, \left(\boldsymbol{x}_n, \boldsymbol{y}_n\right)\right\}$为训练集，定义$\mathcal{I}(\cdot ; \mathcal{A}, \boldsymbol{\gamma})$表示一个插值器，节点集为$\mathcal{A}=\left\{\boldsymbol{a}_1, \cdots, \boldsymbol{a}_m\right\}$，节点处的函数值为γ =(γ₁, …, γ_m)′，从而插值器通过该m个节点，即满足$\mathcal{I}\left(\boldsymbol{a}_j ; \mathcal{A}, \boldsymbol{\gamma}\right)=\gamma_j$。目标函数f可以近似参数化为$\mathcal{I}(\cdot ; \mathcal{A}, \boldsymbol{\gamma})$，从而可以得到f的估计

$ \hat{f}(\boldsymbol{x})=\mathcal{I}(\boldsymbol{x} ; \mathcal{A}, \hat{\boldsymbol{\gamma}}) . $

(1)

其中参数$\boldsymbol{\gamma}=f_{\mathcal{A}}=\left(\mathrm{f}\left(\boldsymbol{a}_1\right), \cdots, f\left(\boldsymbol{a}_m\right)\right)^{\prime}$的估计$\hat{\boldsymbol{\gamma}}$可以通过最小化损失函数L来获得

$ \hat{\boldsymbol{\gamma}}=\arg \min\limits _{\boldsymbol{\gamma} \in \mathbb{R}^m} \sum\limits_{i=1}^n L\left(y_i, \mathcal{I}\left(\boldsymbol{x}_i ; \mathcal{A}, \boldsymbol{\gamma}\right)\right), $

(2)

这种方法称为重构回归法。一种插值器的选择是高斯过程(Gaussian process, GP)^[10]，也被称为Kriging插值器。假定核函数为k(·, ·), g(x)=(g₁(x), …, g_q(x))′为回归项，满足$\boldsymbol{G}_{\mathcal{A}}$=(g(a₁), …, g(a_m))′_m×q列满秩。记$\boldsymbol{r}_{\mathcal{A}}(\boldsymbol{x})=(k(\boldsymbol{x}, $$\left.\left.\boldsymbol{a}_1\right), \cdots, k\left(\boldsymbol{x}, \boldsymbol{a}_m\right)\right)^{\prime}, \boldsymbol{R}_{\mathcal{A}}=\left(k\left(\boldsymbol{a}_i, \boldsymbol{a}_j\right)\right)_{m \times m}$, 从而插值器变为

$ \mathcal{I}(\boldsymbol{x} ; \mathcal{A}, \boldsymbol{\gamma})=\boldsymbol{\gamma}^{\prime} \boldsymbol{b}(\boldsymbol{x}) . $

(3)

其中

$ \boldsymbol{b}(\boldsymbol{x})=\boldsymbol{U} \boldsymbol{g}(\boldsymbol{x})+\boldsymbol{V} \boldsymbol{r}_{\mathcal{A}}(\boldsymbol{x}), $

(4a)

$ \boldsymbol{U}=\boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{G}_{\mathcal{A}}\left(\boldsymbol{G}_{\mathcal{A}}^{\prime} \boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{G}_{\mathcal{A}}\right)^{-1}, $

(4b)

$ \boldsymbol{V}=\left[\boldsymbol{I}_m-\boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{G}_{\mathcal{A}}\left(\boldsymbol{G}_{\mathcal{A}}^{\prime} \boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{G}_{\mathcal{A}}\right)^{-1} \boldsymbol{G}_{\mathcal{A}}^{\prime}\right] \boldsymbol{R}_{\mathcal{A}}^{-1} . $

(4c)

其中I表示单位矩阵，当选取回归项g(x)= 0时，该插值器简化为

$ \mathcal{I}(\boldsymbol{x} ; \mathcal{A}, \boldsymbol{\gamma})=\boldsymbol{\gamma}^{\prime} \boldsymbol{R}_{\mathcal{A}}^{-1} \boldsymbol{r}_{\mathcal{A}}(\boldsymbol{x}), $

(5)

当选择GP时，对于光滑函数，可以以指数阶收敛^[9]。

在函数重构回归法中，目标函数由通过所选节点的插值器重构出来，因此节点集的选择非常重要，Xiong^[9]建议在函数重构回归法中使用m=10d，同时Mu和Xiong^[11]建议选出最小化如下准则

$ \operatorname{cr}(\mathcal{A})=\max\limits _{1 \leqslant i \leqslant j \leqslant m} \sum\limits_{l=1}^d \frac{1}{\left|a_{i l}-a_{j l}\right|} $

(6)

的集合作为节点集$\mathcal{A}$，其中a_j=(a_j1, …, a_jd)′，j=1, …, m。该准则可以很好地平衡空间填充性与低维投影性质。实际操作中可以产生多组集合，从中选取使得上式达到最小的集合作为节点集。同时为选取更好的节点集，Xiong^[9]还提出一种序贯选择节点集的方法。该方法基于一种直观的想法：若在某个点处的损失函数值较大，则表明在该点处的预测结果不够准确，一种最直接的补救措施是将该点纳入到节点集中，即通过以下准则来序贯选择节点

$ \begin{aligned} \boldsymbol{a}_{m+1} & =\arg \max _{\boldsymbol{x}_i \in \mathcal{X} \backslash \mathcal{A}} L\left(y_i, \hat{f}\left(\boldsymbol{x}_i\right)\right) \\ & =\arg \max _{\boldsymbol{x}_i \in \mathcal{X} \backslash \mathcal{A}} L\left(y_i, \hat{\boldsymbol{\gamma}}^{\prime} \boldsymbol{b}\left(\boldsymbol{x}_i\right)\right), \end{aligned} $

(7)

其中: $\hat{\boldsymbol{\gamma}}$为γ的估计，b(x)由式4(a)所定义。

我们将该方法用于密度估计，假设π(x)为未知的密度函数，$\mathcal{X}=\left\{\boldsymbol{x}_1, \cdots, \boldsymbol{x}_n\right\}$为训练集，S⊂ $\mathbb{R}^d$为π的支撑，假定π(x)满足如下形式^[12]

$ \pi(\boldsymbol{x})=\frac{\mathrm{e}^{f(\boldsymbol{x})}}{\int_{\mathrm{S}} \mathrm{e}^{f\left(\boldsymbol{x}^{\prime}\right)} \mathrm{d} \boldsymbol{x}^{\prime}}, $

(8)

其中f(x)为一个未知函数，取指数保证了非负，分母保证了积分为1。运用上述函数重构回归法估计f(x)

$ \hat{f}(\boldsymbol{x})=\hat{\boldsymbol{\gamma}}^{\prime} \boldsymbol{b}(\boldsymbol{x}), $

(9)

其中$\hat{\boldsymbol{\gamma}}$为对应参数的估计，可最小化负对数似然损失得到

$ \hat{\boldsymbol{\gamma}}=\underset{\boldsymbol{\gamma} \in \mathbb{R}^m}{\operatorname{argmin}}-\sum\limits_{i=1}^n \boldsymbol{\gamma}^{\prime} \boldsymbol{b}\left(\boldsymbol{x}_i\right)+n \log \left(\int_S \mathrm{e}^{\boldsymbol{\gamma}^{\prime} \boldsymbol{b}\left(\boldsymbol{x}^{\prime}\right)} \mathrm{d} \boldsymbol{x}^{\prime}\right) . $

(10)

注意该负对数似然损失中包含依赖于未知参数的积分项，在计算时可用伪蒙特卡洛(quasi-Monte Carlo, QMC)^[13]方法近似。若$\mathcal{W}=\left\{\boldsymbol{w}_1, \cdots, \boldsymbol{w}_s\right\}$为调整后的Halton点列，使用求和平均近似积分，从而有

$ \int_S \mathrm{e}^{\boldsymbol{\gamma}^{\prime} \boldsymbol{b}\left(\boldsymbol{x}^{\prime}\right)} \mathrm{d} \boldsymbol{x}^{\prime}=\frac{1}{s} \sum\limits_{t=1}^s \mathrm{e}^{\boldsymbol{\gamma}^{\prime} \boldsymbol{b}\left(\boldsymbol{w}_t\right)}, $

(11)

若选用各向同性高斯核函数，即

$ k(\boldsymbol{u}, \boldsymbol{v})=\exp \left(-\theta \sum\limits_{i=1}^d\left(\mu_i-\nu_i\right)^2\right) . $

(12)

此时负对数似然中包含参数γ和超参数θ，可以采用块坐标下降法^[14]来同时估计，即在第k次迭代中，θ^(k)通过最小化固定γ = γ^(k-1)的负对数似然来获得。同理，γ^(k)可以通过最小化固定θ=θ^(k)的负对数似然来获得，当γ和θ稳定时迭代停止。

1.2 节点的序贯选择算法

密度估计可视为一个无监督学习问题，即没有响应变量y，因此针对回归的节点序贯选择方法在这里不适用。我们考虑如下想法: 若将某个点添加进节点集后能够显著降低损失函数，则应该选择该点作为下一个节点。考虑将节点视为参数，利用贪心算法来进行序贯选择，具体操作如下。

仍假设π为未知的密度函数，$\mathcal{X}=\left\{\boldsymbol{x}_1, \cdots, \right.$$\left.\boldsymbol{x}_n\right\}$为训练集，S⊂$\mathbb{R}^d$为π的支撑。首先最小化$\operatorname{cr}(\mathcal{A})$准则来选取m=7d个节点，重复1 000次，从中选取使得$\operatorname{cr}(\mathcal{A})$最小的一组节点作为初始节点集$\mathcal{A}^{(0)}=\left\{\boldsymbol{a}_1, \cdots, \boldsymbol{a}_m\right\}$。利用初始节点集构造插值器$\hat{f}^{(0)}(\boldsymbol{x})=\hat{\boldsymbol{\gamma}}^{(0)^{\prime}} \boldsymbol{b}^{(0)}(\boldsymbol{x})$及密度估计$\hat{\pi}^{(0)}(\boldsymbol{x})$，其中$\hat{\boldsymbol{\gamma}}^{(0)}=\left(\hat{\gamma}_1, \cdots, \hat{\gamma}_m\right)^{\prime}$为对应参数的估计，即满足$\hat{f}^{(0)}\left(\boldsymbol{a}_j\right)=\hat{\gamma}_j$。下一个节点通过下述方法来选择：遍历集合$\mathcal{D}=\mathcal{X} \backslash \mathcal{A}^{(0)}$，每次从中选取一个样本点$\boldsymbol{d}_i \in \mathcal{D}$作为下一个节点的候选，更新节点集$\mathcal{A}_{\mathrm{i}}^{(1)}=\left\{\boldsymbol{a}_1, \cdots, \boldsymbol{a}_m, \boldsymbol{d}_i\right\}$，更新b_i⁽¹⁾(x)，固定$\hat{\boldsymbol{\gamma}}_1$, …, $\hat{\boldsymbol{\gamma}}_m$的值，记$\boldsymbol{\gamma}_i^{(1)}=\left(\hat{\gamma}_1, \cdots, \hat{\gamma}_m, \gamma_i\right)^{\prime}$，此时$\boldsymbol{\gamma}_i^{(1)}$中只有γ_i为需要估计的未知参数，仍通过最小化负对数似然来得到，其中积分项仍采用QMC方法近似，$\mathcal{W}$和上一节一致，则有

$ \begin{aligned} \hat{\gamma}_i= & \arg \min _{\gamma_i \in \mathbb{R}}-\sum\limits_{j=1}^n \boldsymbol{\gamma}_{\mathrm{i}}^{(1)^{\prime}} \boldsymbol{b}_i^{(1)}\left(\boldsymbol{x}_j\right)+ \\ & n \log \left(\int_S \mathrm{e}^{\boldsymbol{\gamma}_i^{(1)^{\prime}} \boldsymbol{b}_i^{(1)}\left(\boldsymbol{x}^{\prime}\right)} \mathrm{d} \boldsymbol{x}^{\prime}\right) \\ = & \arg \min _{\gamma_i \in \mathbb{R}}-\sum\limits_{j=1}^n \gamma_i b_i\left(\boldsymbol{x}_j\right)+ \\ & n \log \left(\frac{1}{s} \sum\limits_{t=1}^s \mathrm{e}^{\boldsymbol{\gamma}_i^{(1)^{\prime}} \boldsymbol{b}_i^{(1)}\left(\boldsymbol{w}_t\right)}\right) . \end{aligned} $

(13)

若固定高斯核函数中的超参数θ为上一轮的迭代值，则此时该损失函数中只有γ_i一个未知参数，从而变为一维搜索问题，记录此时的损失函数结果L_i。若L_j=min{L₁, …, L_n-m}，则选择a_m+1= d_j，更新节点集$\mathcal{A}^{(1)}=\mathcal{A}_j^{(1)}$以及密度估计$\hat{\pi}^{(1)}(\boldsymbol{x})$。重复上述步骤找出下一个节点，当负对数似然稳定时迭代停止，返回最终的节点集$\mathcal{A}$。

上述算法涉及多次非线性优化问题，一个好的初始值不仅可以加快收敛速度，还有助于提高优化精度，我们充分利用已知信息来构造初始值。由于γ表示光滑函数在节点处的函数值，从而在每一轮迭代过程中，$\hat{\boldsymbol{\gamma}}$只是小幅的波动。因此在γ_i的一维搜索中，可以根据上一轮节点集所构造的密度估计在当前候选节点的估计值作为搜索的初始值，即$\gamma_i=\hat{\pi}^{(0)}\left(\boldsymbol{d}_i\right)$。模拟表明，通过这种初始值的分配，可以有效提高收敛速度和精度。算法过程概括在表 1中。

其中，第3步中的ε是一个很小的常数，控制着停止的时机，ε的选取是实际应用中非常重要的问题。推荐选择$\varepsilon=\min \left\{1, \frac{1}{4} \log \left|L^{(0)}\right|\right\}$，这种选取可以根据损失函数自适应地调整停止时机。该算法能够有效缓解因节点选择不当而造成的估计效果波动，经过序贯选择后的节点集能够使估计精度稳定在一个很好的水平。对于不同数据集，所需要的节点个数也不相同，尤其在高维情形下，一般准则m=10d不一定能够满足估计的需要，根据节点的序贯选择算法可以自适应地选择节点的个数。然而通过这种方法选出的节点集并不能保证最优，因为在每一轮迭代过程中只考虑了局部最优。可以采用束搜索^[15]的思想进行改进，在每一轮迭代过程中选出多个节点作为候选集，遍历候选集产生下一轮的候选集，并淘汰掉损失函数最大的几个组合，通过这种方法可以选出更优的节点集。

2 密度函数的区间估计

实际中除了密度估计本身外，有时还关心密度函数的区间估计。相对于点估计，区间估计可以对不确定性进行量化，从而给出更全面的信息。考虑到重构密度估计的特点，本节利用Bayes思想给出密度函数逐点的区间估计。先验分布是Bayes思想中的重要部分，决定着后验分布的形式以及估计的结果。由于重构方法中的参数γ有着明确的含义，即未知光滑函数在节点处的函数值，而GP假设有限维分布服从多元高斯分布，因此在选择GP作为插值器时，一个合理的先验分布即为

$ \boldsymbol{\gamma} \sim N\left(\boldsymbol{\mu}_\boldsymbol{\gamma}, \boldsymbol{\Sigma}_\boldsymbol{\gamma}\right) . $

(14)

其中: μ_γ为未知函数在节点处的函数值，可取上一节密度估计中的结果，令$\boldsymbol{\mu}_\gamma=\hat{\boldsymbol{\gamma}}$。同时，令$\boldsymbol{\Sigma}_\gamma=$$\left(k\left(\boldsymbol{a}_i, \boldsymbol{a}_j\right)\right)_{m \times m}$，其中k(a_i, a_j)表示高斯核函数，其中核函数中的超参数θ仍取上一节中的结果，因此在此处我们假设它们均为已知，下面给出可信区间的构造过程。

基于γ的先验，计算得到似然函数

$ \pi(\boldsymbol{X} \mid \boldsymbol{\gamma})=\prod\limits_{i=1}^n \pi\left(\boldsymbol{x}_i \mid \boldsymbol{\gamma}\right)=\frac{\mathrm{e}^{\boldsymbol{\gamma}^{\prime} \sum\limits_{i=1}^n \boldsymbol{b}\left(\boldsymbol{x}_\boldsymbol{i}\right)}}{\left[\int_S \mathrm{e}^{\boldsymbol{\gamma}^{\prime} \boldsymbol{b}\left(\boldsymbol{x}^{\prime}\right)} \mathrm{d} \boldsymbol{x}^{\prime}\right]^n}, $

(15)

从而可以推出后验

$ \begin{gathered} {\pi(\boldsymbol{\gamma} \mid \boldsymbol{X}) \propto}\\ {\frac{\exp \left\{\boldsymbol{\gamma}^{\prime} \sum\limits_{i=1}^n \boldsymbol{b}\left(\boldsymbol{x}_i\right)-\frac{1}{2}\left(\boldsymbol{\gamma}-\boldsymbol{\mu}_{\boldsymbol{\gamma}}\right)^{\prime} \boldsymbol{\Sigma}_{\boldsymbol{\gamma}}^{-1}\left(\boldsymbol{\gamma}-\boldsymbol{\mu}_{\boldsymbol{\gamma}}\right)\right\}}{\left[\int_S \mathrm{e}^{\boldsymbol{\gamma}^{\prime} \boldsymbol{b}\left(\boldsymbol{x}^{\prime}\right)} \mathrm{d} \boldsymbol{x}^{\prime}\right]^n\left|\boldsymbol{\Sigma}_{\boldsymbol{\gamma}}\right|^{\frac{1}{2}}}} . \end{gathered} $

(16)

其中积分项仍然采用QMC方法近似计算。需要注意的是在计算过程中积分项$\left[\int_S \mathrm{e}^{\gamma^{\prime} \boldsymbol{b}\left(\boldsymbol{x}^{\prime}\right)} \mathrm{d} \boldsymbol{x}^{\prime}\right]^n$往往会过大导致超出运算范围，因此可以先对后验取对数，基于所求对数后验分布，采用Metropolis抽样算法取得服从后验分布的样本。在维度较高时，普通抽样算法效率较低，可采用哈密尔顿蒙特卡罗算法^[16]提高抽样速度。详细的抽样算法流程如表 2所示。

其中，常数τ_ρ是一个可以调节的参数，控制着抽样过程中的接受率，在实际应用中是一个重要的问题。Sherlock和Roberts^[17]以及Gelman等^[18]对最优接受率进行过讨论，对于随机游走抽样，在建议分布为高斯分布时，一维情形下的最优接受率约为0.439，维数大于1时最优接受率约为0.234，并给出了理论证明，因此在应用时通过调节τ_ρ的值使得接受率接近最优接受率即可。

根据抽样所得后验样本γ₁, …, γ_q代入上一节构造的密度估计中可以得到q个不同的密度估计π₁(x), …, π_q(x)。用样本的分位数来近似总体分位数，从而在给定x₀处的100(1-α)% 可信区间为

$ \left[\pi_{[(\alpha / 2) n]}\left(\boldsymbol{x}_0\right), \pi_{[(1-\alpha / 2) n]}\left(\boldsymbol{x}_0\right)\right], $

(17)

其中: $\pi_{[(\alpha / 2) n]}\left(\boldsymbol{x}_0\right)$为$\left\{\pi_i\left(\boldsymbol{x}_0\right), i=1, \cdots, q\right\}$的次序统计量中第[(α/2)n]个值，[·]表示取整函数。

3 数值模拟 3.1 模拟数据集

本节在模拟数据集上评估我们提出的节点序贯选择算法，并将其与非序贯重构密度估计以及KDE^[5-6]进行比较。其中非序贯重构密度估计中直接利用最小化$\operatorname{cr}(\mathcal{A})$准则选出m=8d个节点，KDE中的窗宽通过最小化渐进均方误差^[7]得到，所用数据集来自如下几个分布^{[12, 19]}：

① X~0.4N(-2, 0.5²)+0.6N(1, 0.5²),

② X~0.7beta(3, 12)+0.3beta(10, 3),

③ X=(X₁, X₂), X₁~2/3beta(1, 2)+1/3beta(10, 10), X₂~U(0, 1),

④ $\boldsymbol{X} \sim 0.4 N\left(\binom{1}{0},\left(\begin{array}{cc} 1 & -1 \\ -1 & 2 \end{array}\right)\right)+0.6 N\left(\binom{0}{1},\right.$

    $\left.\left(\begin{array}{cc} 2 & -1 \\ -1 & 1 \end{array}\right)\right),$

⑤ $\boldsymbol{X}=\left(X_1, X_2\right), X_1 \sim 0.7 t(5)+0.3 t(5),$

    $3\left(X_2-3\right) / 2 \sim 0.7 t(5)+0.3 t(5),$

⑥ $\boldsymbol{X} \sim 0.3 N\left(\binom{\cos 60^{\circ}}{\sin 60^{\circ}}, \boldsymbol{I} / 8\right)+$

    $0.4 N\left(\binom{\cos 120^{\circ}}{\sin 120^{\circ}}, \boldsymbol{I} / 8\right)+0.3 N\left(\binom{\cos 180^{\circ}}{\sin 180^{\circ}}, \boldsymbol{I} / 8\right),$

⑦ $\boldsymbol{X} \sim 0.4 N\left(\left(\begin{array}{l} 1 \\ 0 \\ 1 \end{array}\right),\left(\begin{array}{ccc} 1 / 4 & -1 / 4 & 1 / 4 \\ -1 / 4 & 1 / 2 & -1 / 4 \\ 1 / 4 & -1 / 4 & 1 / 2 \end{array}\right)\right)+$

    $0.6 N\left(\left(\begin{array}{l} 0 \\ 1 \\ 0 \end{array}\right),\left(\begin{array}{ccc} 1 / 2 & -1 / 4 & 1 / 4 \\ -1 / 4 & 1 / 2 & -1 / 4 \\ 1 / 4 & -1 / 4 & 1 / 2 \end{array}\right)\right),$

⑧ $\boldsymbol{X} \sim 0.4 N\left(\left(\begin{array}{l} 1 \\ 0 \\ 1 \\ 0 \end{array}\right), \boldsymbol{I} / 4\right)+0.6 N\left(\left(\begin{array}{l} 0 \\ 1 \\ 0 \\ 1 \end{array}\right), \boldsymbol{I} / 4\right) .$

使用KL散度^[20]衡量3种方法的精度，对于 $\hat{\pi}$和π之间的KL散度定义为

$ \int_{\mathbb{R}} \pi(\boldsymbol{x}) \log \frac{\pi(\boldsymbol{x})}{\hat{\pi}(\boldsymbol{x})} \mathrm{d} \boldsymbol{x} . $

(18)

使用蒙特卡洛算法近似估计该积分，即从分布π(x) 中抽取10 000个独立同分布的样本{s₁, …, s_{10 000} }并计算

$ \frac{1}{10\;000} \sum\limits_{i=1}^{10\;000} \log \frac{\pi\left(\boldsymbol{s}_i\right)}{\hat{\pi}\left(\boldsymbol{s}_i\right)} . $

(19)

将实验重复20次，分别计算3种方法的平均KL散度，结果总结在表 3中。

由表 3可以看出，通过序贯方法选择出的节点集使得估计结果更具稳定性，同时可以使估计精度进一步提升，估计效果优于非序贯方法产生的节点集。

下面选取数据集①中的密度函数来构建置信度为95% 的区间估计。首先基于序贯选择算法得到$\hat{\boldsymbol{\gamma}}$和$\hat{\theta}$，然后利用Metropolis抽样算法得到γ的后验样本，带入密度估计中用样本分位数估计总体分位数。为比较效果，利用Hall^[21]提到的方法，通过构造一个欠光滑的KDE使偏差收敛到0，再利用bootstrap方法构造置信区间。因为这2种方法构造的区间含义不同，这里只用来比较大致的效果。均匀取100个点构造区间估计，用这2种方法重复构造10次，分别取区间上限的平均值和区间下限的平均值作为最后的区间估计，结果如图 1和图 2所示。

图中实线表示真实密度函数图像，上下虚线分别表示区间上限和区间下限。可以看出, 基于重构方法构造的区间估计更好地覆盖了真实密度函数图像，同时有着较小的区间长度，这得益于重构密度估计良好的拟合效果。相反，基于KDE构造的区间较宽，且在尾部有着较大波动。此外，由于KDE的局部特性，其在高维情形估计非常不稳定，此时构造的置信区间也随之变得较不稳定，而基于重构方法构造的区间估计可以很自然地应用于高维情形。

3.2 实例数据

本节展示序贯重构密度估计在地理空间数据中的应用，原始数据集来自Phillips等^[22]，可由Python数据库datasets中的fetch_species_ditributions函数加载。

建立物种的地理分布模型是生物学中的一个重要问题，密度函数作为一种平滑表示，可以更为直观清晰地表示物种的分布。此数据集记录了褐喉树懒和森林小稻鼠2种南美哺乳动物的地理分布，其中褐喉树懒有926个样本点，森林小稻鼠有698个样本点。图 3绘制了2个物种的分布散点图。

可以看出这并不能很好地了解物种的密度，因为有大量的点重叠在了一起。下面将用序贯重构密度估计从分布数据中计算这2个物种的分布密度。密度估计的结果用密度等高线来表示，图 4展示了2个物种分布的密度估计结果。

从密度等高线图中可以很清晰直观地看出2个物种的分布情况，颜色越深的地方即表示物种分布越密集。褐喉树懒主要分布在南美洲与北美洲接壤处以及中部亚马逊河沿岸区域，在哥伦比亚西北部与北美洲接壤处最为集中。森林小稻鼠主要分布在南美洲西北部的沿海地区，在哥伦比亚和厄瓜多尔的沿海地区最为集中。

4 结论

本文提出用于重构密度估计的节点序贯选择方法，由于密度估计可视为无监督学习问题，即没有响应变量y，因此给节点的序贯选择带来了较大困难。通过将节点视作参数，同时基于这样一个想法：若将某个点添加进节点集后能够使得损失函数显著减小，则应考虑选择该点作为下一个节点，从而遍历训练集找出使得损失函数最小的候选点作为下一个节点，直至损失函数趋于稳定。为提升序贯选择的精度和速度，还给出一些初始值的设定方法。最后，根据参数γ的实际含义构造先验分布，并通过Metropolis抽样方法得到后验分布的样本，构造了密度函数逐点的区间估计。这种序贯选择算法可以很自然地应用到其他无监督学习问题，为此类问题的序贯选择算法提供了一个框架。模拟结果表明，通过序贯算法选出节点构造的密度估计效果优于非序贯方法，且具有稳健性，可以很好地缓解由于节点集选择不同而造成的精度波动。在高维情形下，通用准则不一定能够满足估计的需求，序贯方法可以自适应地选择节点个数。