核密度估计基础-Part I

核平滑方法理论-I

@(机器学习)[MachineLearning, Econometrics]

0. Introduction

核密度估计是一种非参数估计方法，在机器学习领域，是一种非监督性学习方法。用于从给定分布的样本重建总体的分布函数。

优点：

• 非参数：假设少，不假设样本服从任何分布

缺点：

• 计算量：比起参数估计，非参数估计运算量大很多

1. 核密度估计(Kernel Density Estimation)

1.1 单变量(Univariable)密度估计

1.1.1 单变量的核密度估计

定理 1.1: 均匀核估计量
$$\hat{f}(x)=\frac{1}{nh}\sum _{i=1}^{n}k\left(\frac{X_i-x}{h}\right)$$
要$\hat{f}$是$f$的一致估计量，只要核函数$k(\cdot )$满足

1. 归一化, $\int k(v)dv=1$
2. 对称性, $k(v)=k(-v)$
3. 二阶矩有限, $\int v^{2}k(v)dv<\infty$

并且可证明$\hat{f}$有一个渐进正态分布，也就是说$\hat{f}(x)$统计量服从中心极限定理。

注意到
\begin{align*}
\hat{f}(x)&=\frac{1}{nh}\sum_{i=1}^{n}k\left(\frac{X_i-x}{h}\right)\
&=(\hat{F}\star k_h)(x)
\end{align*}
其中$\star$表示卷积，$\hat{F}(x)={\sum }_{i=1}^n\delta (x-X_i)$，是一堆针刺。这也就是说，通过总体的密度分布$\hat{f}$是利用宽度为$h$的核函数$k_h$平滑了针板函数$\hat{F}(x)$得到的。

均方误差(Mean Square Error)分析
\begin{align*}
MSE[\hat{f}(x)]&\equiv E\left[\left(\hat{f}-f(x)\right)^2\right]\
&=var\left[\hat{f}(x)\right]+\left[E\left(\hat{f}(x)\right)-f(x)\right]^2\
&=var\left[\hat{f}(x)\right]+bias^2\left[\hat{f}(x)\right]
\end{align*}
可以利用Taylor展开方便的证明$\hat{f}(x)$具有均方误差一致收敛速度满足下面定理

定理 1.2：设三阶可微概率密度$f(x)$有一组i.i.d.的$n$个观测值{Xn}\{X_n\}{Xn​}。核函数$k(\cdot )$满足归一性，对称性和二阶矩存在，且当$n\to \infty$时，有$h$宏观无穷小$h\to 0$，微观无穷大$nh\to \infty$。则对于$x\in \text{supp}(X)$
\begin{align*}
MSE[\hat{f}(x)]&=bias^2[\hat{f}(x)]+var[\hat{f}(x)] \
&=\frac{h^{4}{4}\left[\kappa_2f}{(2)}(x)\right]^2+\frac{\kappa f(x)}{nh}+o(h^4+(nh){-1})
\end{align*}
其中$\kappa =\int k^2(x)dx$，${\kappa }_2=\int x^{2}k(x)dx$由核函数性质决定。并且$|f^{(1)}(x)|<\infty$，$\int |x^{3}k(x)|dx<\infty$。

因此$\hat{f}(x)$在均方误差意义下一致收敛于$f(x)$。

更近一步，如果将MSE作为判据，为了使MSE最小($\frac{dMSE(\hat{f}(x))}{dh}=0$)，应该选取的核宽度为
$$h_{opt}=c(x)n^{-1/5}$$
其中c(x)={κ f(x)(κ2 f(2)(x))2}1/5c(x)=\{\frac{\kappa\,f(x)}{\left(\kappa_2\,f^{(2)}(x)\right)^2}\}^{1/5}c(x)={(κ2​f(2)(x))2κf(x)​}1/5

注意到上面的窗口宽度随着$x$变化的函数，如果希望使用固定窗口宽度，我们选择固定核宽度的积分均方误差作为评判标准，即估计密度函数和总体密度之间的期望希尔伯特距离
\begin{align*}
IMSE[\hat{f}(x)]&\equiv \int E\left[\left(\hat{f}(x)-f(x)\right)^2\right],dx\
&=\frac{1}{4}h^4\kappa_22\int \left[,f^{{(2)}(x),\right]}2,dx+\frac{\kappa}{nh}+o(h^4+(nh){-1})
\end{align*}
在这个意义下，可以求得是IMSE最小的优化$h_{opt}$
$$h_{opt}=c_0{n}^{-1/5}$$
其中$c_0={\kappa }_2^{-2/5}{\kappa }^{1/5}{\left\{\int [f^{(2)}(x){]}^{2}dx\right\}}^{-1/5}>0$。

1.1.2 窗宽选择

1. 插入法(plug-in methods)
   为了求出在IMSE条件下最有的窗宽，需要确定常数$c_0$中的$\int {\left[f^{(2)}(x)\right]}^{2}dx$。由于$f$是未知的，所以这个量无法事先知道。如果选择一个$h$初始的“试验值”(pilot value)，然后将这个值代入$h_{opt}$的计算式求出的优化$h$，则这种方法称为“插入法”(plug-in methods)。
   Silverman(1986)提出假定$f$是一个以${\sigma }^2$为方差的正态分布，则其二阶导可确定，$\int {\left[f^{(2)}(x)\right]}^{2}dx=\frac{3}{8\sqrt{\pi }{\sigma }^5}$，代入优化窗宽，可以得到试验窗宽估计
   $$h_{pilot}=(4\pi {)}^{-1/10}[(3/8)\sqrt{\pi }{]}^{-1/5}\sigma n^{-1/5}\approx 1.06\sigma n^{-1/5}$$
   用此试验值进一步迭代计算$\int {\left[{\hat{f}}^{(2)}(x)\right]}^{2}dx$，定出最终的优化结果$h_{opt}$。
   Silverman还提出一种更加稳健的分散程度度量，就是用min⁡{σ,q1/4/1.34}\min\{\sigma, q_{1/4}/1.34\}min{σ,q1/4​/1.34}来代替$\sigma$，其中$q_{1/4}$表示四分位矩。

2. 交错鉴定法
   交错鉴定法是一种完全由数据驱动的方法，其核心在于用一部分样本拟合模型来检验另一部分样本的拟合程度。通过不断改变训练集合测试集，来评价模型的好坏。当每次都只留一个样本作为检验对象，其他样本均做训练集时，所得到的估计量称为去一估计量(leave-one-out estimator)。
   通过这种方法，我们可以来估计$\hat{f}$和$f$的希尔伯特距离，并以距离作为判据来选择窗宽，这种方法称为最小二乘交叉检验。
   \begin{align*}
   L(\hat{f}, f)&=\int[\hat{f}(x)-f(x)]^2,dx\
   &=\int\hat{f}(x)^2dx-2\int\hat{f}(x)f(x)dx+\int f(x)^2dx\
   \end{align*}
   其中第三项和$\hat{f}$无关，视为常数
   $$\int f(x{)}^{2}dx=C$$
   第二项采用去一估计量估计，即
   $$\int \hat{f}(x)f(x)dx=E_X\left[\hat{f}(X)\right]=\frac{1}{n}\sum _{i=1}^n{\hat{f}}_{-i}(X_i)+O(n^{-1/2})$$
   其中$E_x[\cdot ]$是对$x$求期望，用来区别对观测量$X_i$求期望。在$X_i$处的去一估计量${\hat{f}}_{-i}(X_i)$定义为
   $${\hat{f}}_{-i}(X_i)=\frac{1}{(n-1)h}\sum _{j\ne i}^{n}k\left(\frac{X_i-X_j}{h}\right)$$
   表示用除了$X_i$这个观测量外的其他观测量来估计$X_i$处的密度函数。
   第一项直接代入$\hat{f}(x)$的估计式，可以得到
   \begin{align*}
   \int\hat{f}(x)^{2dx&=\int\left[\frac{1}{nh}\sum_{i=1}}nk\left(\frac{X_i-x}{h}\right)\right]^2dx\
   %&=\frac{1}{n^2h2}\sum_{i=1}^n\sum_{j=1}n\int k\left(\frac{X_i-x}{h}\right)k\left(\frac{X_j-x}{h}\right)dx\
   %&=\frac{1}{n^2h2}\sum_{i=1}^n\sum_{j=1}n\int k\left(\frac{x}{h}\right)k\left(\frac{x+X_i-X_j}{h}\right)\cdot h, d\left(\frac{x}{h}\right)\
   &=\frac{1}{n^{2h}\sum_{i=1}}n\sum_{j=1}^n\bar{k}\left(\frac{X_i-X_j}{h}\right)
   \end{align*}
   其中$\bar{k}(t)=\int k(x)k(t-x)dx$是$k(\cdot )$的重卷积核(two-fold convolution)，一般是两个独立同分布的随机变量之和的分布。可证明，$\bar{k}(\cdot )$也是偶函数。

定理 1.3 总体分布函数为$f(x)$，通过去一核估计交叉检验得到的估计量$\hat{f}$的积分平方误差$CV$为
$$C{V}_f(h)=\frac{1}{n^{2}h}\sum _{i=1}^n\sum _{j=1}^n\bar{k}\left(\frac{X_i-X_j}{h}\right)-\frac{2}{n(n-1)h}\sum _{i=1}^n\sum _{j\ne i}^{n}k\left(\frac{X_i-X_j}{h}\right)+C$$
其中$\bar{k}(t)=\int k(x)k(t-x)dx$是$k(\cdot )$的重卷积核。

可以通过成熟的数值算法对$C{V}_f(h)$进行优化求解得到使交叉检验$C{V}_f$最小的核宽度$h$。
将$C{V}_f(h)$的首项提出，并使首项最小，会发现得到的最优解退化为IMSE最优解的情形。

除了最小二乘方法，还可以使用最概然交叉检验。根据玻尔兹曼熵定义，这种方法以最大化去一核最概然函数的对数为标准来选取$h$，即
$$\mathcal{L}=k\ln L=k\sum _{i=1}^n\ln \left[{\hat{f}}_{-i}(X_i)\right]$$
其中$k$为玻尔兹曼常数。这种方法受到尾部行为影响严重，对厚尾分布会引起不一致的结果，因此最概然交错检验不太流行。

1.2 单变量累计分布函数

1.2.1 累计分布函数的核估计

为了得到平滑的CDF估计量，我们从核函数出发，将密度分布函数估计进行积分
$$\hat{F}(x)={\int }_{-\infty }^x\hat{f}(x)dx=\frac{1}{n}\sum _{i=1}^{n}G\left(\frac{x-X_i}{h}\right)$$
其中$G(x)={\int }_{-\infty }^{x}k(x)dx$是核的累计分布函数。其均方误差有下面定理给出

定理 1.4：总体的累计分布函数$F(x)$二阶连续可微，且二阶倒数Holder连续，$k(x)$为对称的核函数，$G(x)={\int }_{\infty }^{x}k(x)$为核积分函数。则当$n\to \infty$时，
\begin{align*}
MSE[\hat{F}]&=bias[\hat{F}]^2+var[\hat{F}] \
& = \left{ \frac{1}{2}\kappa_2h^2F{(2)}(x) + o\left(h^{2\right)\right}}2 \
& + \left{\frac{1}{n}F(x)[1-F(x)]-\frac{1}{n}\alpha_0f(x)h+o\left(\frac{h}{n}\right)\right}\
&=c_0(x)n^{{-1}-c_1(x)hn}{-1}+c_2(x)h^4+o(h4+hn^{-1})
\end{align*}
其中系数项为
\begin{align*}
c_0(x)&=F(x)[1-F(x)]\
c_1(x)&=\alpha_0f(x)\
c_2(x)&=\left[\frac{\kappa_2}{2}F^{{(2)}(x)\right]}2\
\alpha_0&=2\int xG(x)k(x)dx\
\kappa_2&=\int x^2k(x)dx
\end{align*}
系数由总体分布函数$F(x)$和核确定$k(x)$。

因此，可以容易的$\hat{F}$到积分均方误差IMSE
\begin{align*}
IMSE(\hat{F})&=\int E[\hat{F}(x)-F(x)]^2dx\
&=C_0n^{-1}-C_1hn{-1}+C_2h^4+o(h4+hn^{-1})
\end{align*}
其中$C_i=\int c_i(x)dx$是和$x$无关的常数。
首项最小化可以的到优化的核宽度选择
$$h_{opt}={\left[\frac{C_1}{4C_2}\right]}^{1/3}{n}^{-1/3}$$
这比密度估计($n^{-1/5}$)收敛速度要快。

渐进正态特性，根据Liapunov中心极限定理，分布上
$$\sqrt{n}[\hat{F}-F]\sim \mathcal{N}\left(0,F(x)[1-F(x))]\right)$$
误差满足正态分布。

1.2.2 窗宽选择

交叉检验法：累计分布函数估计$\hat{F}(x)$的交叉检验函数定义如下
$$C{V}_F(h)=\frac{1}{n}\sum _{i=1}^n\int {\left[1(X_i\le x)-{\hat{F}}_{-i}(x)\right]}^{2}dx$$
其中$1$是示性函数，${\hat{F}}_{-i}(x)=\frac{1}{n-1}{\sum }_{j\ne i}G\left(\frac{x-X_j}{h}\right)$为去一核估计量。

可以证明交叉检验函数期望的首项和$IMSE(\hat{F})$的首项相同。因此用交叉检验和用IMSE得到的效果相同。

1.3 多变量(Multivariable)联合分布密度估计

1.3.1 联合分布的核估计

当我们考察的对象从标量随机变量扩充为$q$维随机向量时，我们需要的估计的密度分布函数就也称为了联合密度分布。我们将问题形式化如下，假定有$n$个$q$维随机向量{Xn}\{X_n\}{Xn​}且i.i.d服从联合密度函数$f(x_1,x_2,\dots ,x_q)$，记$X_{is}$为$X_i$的第$s$个分量。即

s	=		1	2	…	q
$X_1$	=	(	$X_{11}$,	$X_{12}$,	…,	$X_{1q}$	)
$X_2$	=	(	$X_{21}$,	$X_{22}$,	…,	$X_{2q}$	)
…	=	(	…,	…,	…,	…	)
$X_n$	=	(	$X_{n1}$,	$X_{n2}$,	…,	$X_{nq}$	)

联合分布的核函数通过单变量核函数的乘积构造，这样的构造的联合密度核函数是假设$q$个核相互独立时的联合分布函数，但$X$的分量之间并不需要限制是独立的。也就是说，$X$分量之间有依赖时也可以通过这样的核估计出来。我们用下面的方法来估计联合概率密度$f(x)$
$$\hat{f}=\frac{1}{n{h}_1\cdots h_q}\sum _{i=1}^{n}K\left(\frac{X_i-x}{h}\right)$$
其中，核函数
$$K\left(\frac{X_i-x}{h}\right)=\prod _{i=1}^{q}k\left(\frac{X_i-x}{h_i}\right)$$
而$k(x)$则是单变量核函数。

均方误差的计算类似于单变量的其概况，可以得到

定理 1.5：设三阶梯度存在的$q$维联合概率密度分布函数$f(x)\equiv f(x_1,x_2,\dots ,x_q)$有一组i.i.d.的$n$个观测值{Xn∈Rq}\{X_n\in\mathbb{R}^q\}{Xn​∈Rq}。核函数$K(x)$为单变量核函数之积。且当$n\to \infty$时，有格子体积宏观无穷小${\max }_i{h}_i\to 0$，微观无穷大$n{h}_1{h}_2\cdots h_q\to \infty$。则对于$x\in \text{supp}(X)$
\begin{align*}
MSE[\hat{f}(x)]&=bias^2[\hat{f}(x)]+var[\hat{f}(x)] \
&=\left{\frac{\kappa_2}{2}\sum_{s=1}^qh_s2\frac{\partial^2 f(x)}{\partial x_s^{2}+O\left(\sum_{s=1}}qh_s^3 \right)\right}^2\
&+\left{\frac{1}{nh_1h_2\cdots h_q}\left[\kappa^{qf(x)+O\left(\sum_{s=1}}qh_s^2\right)\right]\right}\
&=O\left(\left(\sum_{s=1}^qh_s2\right)^2+(nh_1h_2\cdots h_q)^{-1}\right)\
&=O(L^4+(nV){-1})
\end{align*}
其中$\kappa =\int k^2(x)dx$，${\kappa }_2=\int x^{2}k(x)dx$由单变量核函数性质决定。$L$为核宽度超立方体的对角线长度，而$V$为超立方体的体积。

渐进正态性讨论
如果$n\to \infty$，格子宏观无穷小${\max }_i{h}_i\to 0$，微观无穷大$nV\to \infty$时，并且$nV{\sum }_{s=1}^q{h}_s^6\to 0$，密度估计量具有渐进正态性。
$$\hat{f}(x)-f(x)-bias[\hat{f}(x)]\to \mathcal{N}(0,\frac{{\kappa }^{q}f(x)}{nV})$$
即其无偏误误差服从均值为0的正态分布。

1.3.2 窗框选择

插入法
优化的核宽度选择应当平衡偏误和方差，也就是说，对于所有的$s$应当有
$$h_s^4=O\left((n{h}_1{h}_2\cdots h_q{)}^{-1}\right)$$
因此，优化的$h_s$应满足
$$h_s=c_s{n}^{-1/(q+4)}$$
在应用中，需要对常数$c_s$进行选择，经验法则山，一般选取$c_s=1.06$。但由于总体的分布函数可能各向异性，所以这样一概而论的常数缺乏灵活性。

对于插入法，一般通过$\hat{f}(x)$的偏误和方法首项进行估计，其中包含了总体分布$f(x)$和二阶偏导数，这在高维情况中是复杂的。在实际中插入法没有广泛使用，也不推荐使用。

交叉检验法
自然地将一维交叉检验函数扩充到高维的情况，定义交叉检验目标函数为
\begin{align*}
CV_f(h_1,\ldots,h_q)&=\frac{1}{n^2}\sum_{i=1}n\sum_{j=1}^n\overline{K}h(X_i, X_j)\
&\quad-\frac{2}{n(n-1)}\sum{i=1}^n\sum_{j\ne i}^n K_h(X_i,X_j)
\end{align*}
其中
\begin{align*}
K_h(X_i,X_j)=\prod_{s=1}^q\frac{1}{h_s}k\left(\frac{X_{is}-X_{js}}{h_s}\right)\
\overline{K}h(X_i,X_j)=\prod{s=1}^q\frac{1}{h_s}\bar{k}\left(\frac{X_{is}-X_{js}}{h_s}\right)
\end{align*}
是单变量版本的乘积形式。可以通过数值方法来寻求目标函数的最小化。

从理论分析上交叉检验目标函数$C{V}_f(h_1,\dots ,h_q)$的首项通过下式给出
$$C{V}_{f_0}(h_1,h_2,\dots ,h_q)=\int {\left[\sum _{s=1}^q{B}_s(x)h_s^2\right]}^{2}dx+\frac{{\kappa }^q}{n{h}_1{h}_2\cdots h_q}$$
其中$B_s(x)=\frac{{\kappa }_2}{2}\frac{{\partial }^{2}f(x)}{\partial x_s^2}$， $\kappa =\int k^2(x)dx$，${\kappa }_2=\int x^{2}k(x)dx$。
为了分离出样本数$n$的影响，我们定义$a_s=h_s{n}^{1/(q+4)}$，代换$h_s$得到
$$C{V}_{f_0}(h_1,h_2,\dots ,h_q)=n^{-1/(q+4)}\chi (a_1,a_2,\dots ,a_q)$$
其中$\chi (a_1,a_2,\dots ,a_q)$适合$n$无关的常数，定义为
$$\chi (h_1,h_2,\dots ,h_q)=\int {\left[\sum _{s=1}^q{B}_s(x)a_s^2\right]}^{2}dx+\frac{{\kappa }^q}{a_1{a}_2\cdots a_q}$$

因此可以看到，最大化首项的$h_s$应满足$h_s=O(n^{-1/(q+4)})$。同时可以证明$C{V}_{f_0}$的首项也是$E[C{V}_f]$的首项，也就说说，最优化$h_s$也使得积分均方误差的首项最小化。

最概然交叉检验和单变量情况通过最大化熵来给出最优化窗宽，虽然执行简单，单依然会有厚尾分布时出现缺陷的情况，会出现过度平滑。

1.4 高阶核函数

定义 1.1： 一个$\nu$阶核函数($\nu \ge 2$)应满足如下条件

1. 归一化, $\int k(x)dx=1$
2. 低阶矩为0, $\int x^{l}k(x)dx=0$，$l=1,\cdots ,\nu -1$
3. $\nu$阶矩有限, $\int x^{\nu }k(x)dx={\kappa }_{\nu }\ne 0<\infty$

则称核函数$k(\cdot )$为$\nu$阶核函数。

通常使用的核都属于二阶核函数$\nu =2$。与二阶核类似，对于总体分布函数$f(x)$是$\nu$阶可微，所有的维度使用相同阶核函数时，可以证明
\begin{align*}
bias[\hat{f}(x)]&=O\left(\sum_{s=1}^qh_s\nu\right)\
var[\hat{f}(x)]&=O((nh_1h_2\cdots h_q)^{-1})
\end{align*}
利用这个结果，可以得到均方差和估计的误差

定理 1.6: 对于一个$\nu$阶核函数，$nu\ge 2$，其误差由下式给出
\begin{align*}
MSE[\hat{f}(x)]&=O\left(\sum_{s=1}^qh_s{2\nu}+(nh_1h_2\cdots h_q)^{-1}\right)\
\hat{f}(x)-f(x)&=O_p\left(\sum_{s=1}^qh_s{\nu}+(nh_1h_2\cdots h_q)^{-1/2}\right)
\end{align*}
利用一个高阶和可以同时较少偏误和方法。

值得注意的是，对于$\nu >2$，不存在非负核函数。也就意味着，我们有可能得到负的密度估计。对于有限样本来说，一个非负的二阶核函数经常比高阶核函数得到更稳定的结果。因此，高阶核函数经常被用于理论目次，而不太在实践中运用。

高阶核函数可以通过低阶核函数与多项式乘积的形式进行构造，通过矩约束求解多项式系数。

1.5 展望

1. 放开窗口宽度常数限制，使用变长窗口宽度。
2. 采用变换分布，消除偏度的影响。

参考资料

[1] Q. Li & J. S. Racine, Nonparametric Econometrics Theory and Practice, Peking University Press, 2007
[2] T. Hastie, R. Tibshirani & J. Friedman, The Elements of Statistical Learning, Second Edition, Springer, 2009
[3] B. Silverman, Density Estimation for Statistics and Data Analysis, Springer, 1986