核平滑方法理论-I

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

0. Introduction

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

优点：

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

缺点：

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

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

1.1 单变量(Univariable)密度估计

1.1.1 单变量的核密度估计

定理 1.1: 均匀核估计量
f ^ ( x ) = 1 n h ∑ i = 1 n k ( X i − x h ) \hat{f}(x)=\frac{1}{nh}\sum_{i=1}^{n}k\left(\frac{X_i-x}{h}\right) f^(x)=nh1i=1∑nk(hXi−x)
要 f ^ \hat{f} f^是 f f f的一致估计量，只要核函数 k ( ⋅ ) k(\cdot) k(⋅)满足

归一化, ∫ k ( v ) d v = 1 \int k(v)dv=1 ∫k(v)dv=1

对称性, k ( v ) = k ( − v ) k(v)=k(-v) k(v)=k(−v)

二阶矩有限, ∫ v 2 k ( v ) d v < ∞ \int v^2k(v)dv<\infty ∫v2k(v)dv<∞

并且可证明 f ^ \hat{f} f^有一个渐进正态分布，也就是说 f ^ ( x ) \hat{f}(x) 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 ⋆表示卷积， F ^ ( x ) = ∑ i = 1 n δ ( x − X i ) \hat{F}(x)=\sum_{i=1}^n\delta(x-X_i) F^(x)=∑i=1nδ(x−Xi)，是一堆针刺。这也就是说，通过总体的密度分布 f ^ \hat{f} f^是利用宽度为 h h h的核函数 k h k_h kh平滑了针板函数 F ^ ( x ) \hat F(x) 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展开方便的证明 f ^ ( x ) \hat{f}(x) f^(x)具有均方误差一致收敛速度满足下面定理

定理 1.2：设三阶可微概率密度 f ( x ) f(x) f(x)有一组i.i.d.的 n n n个观测值 { X n } \{X_n\} {Xn}。核函数 k ( ⋅ ) k(\cdot) k(⋅)满足归一性，对称性和二阶矩存在，且当 n → ∞ n\to\infty n→∞时，有 h h h宏观无穷小 h → 0 h\to 0 h→0，微观无穷大 n h → ∞ nh\to\infty nh→∞。则对于 x ∈ supp ( X ) x\in\text{supp}(X) x∈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*}
其中 κ = ∫ k 2 ( x ) d x \kappa=\int k^2(x)dx κ=∫k2(x)dx， κ 2 = ∫ x 2 k ( x ) d x \kappa_2=\int x^2k(x)dx κ2=∫x2k(x)dx由核函数性质决定。并且 ∣ f ( 1 ) ( x ) ∣ < ∞ |f^{(1)}(x)|<\infty ∣f(1)(x)∣<∞， ∫ ∣ x 3 k ( x ) ∣ d x < ∞ \int|x^3k(x)|dx<\infty ∫∣x3k(x)∣dx<∞。

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

更近一步，如果将MSE作为判据，为了使MSE最小( d M S E ( f ^ ( x ) ) d h = 0 \frac{d\,MSE(\hat{f}(x))}{dh}=0 dhdMSE(f^(x))=0)，应该选取的核宽度为
h o p t = c ( x ) n − 1 / 5 h_{opt}=c(x)n^{-1/5} hopt=c(x)n−1/5
其中 c ( x ) = { κ f ( x ) ( κ 2 f ( 2 ) ( x ) ) 2 } 1 / 5 c(x)=\{\frac{\kappa\,f(x)}{\left(\kappa_2\,f^{(2)}(x)\right)^2}\}^{1/5} c(x)={(κ2f(2)(x))2κf(x)}1/5

注意到上面的窗口宽度随着 x x 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 o p t h_{opt} hopt
h o p t = c 0 n − 1 / 5 h_{opt}=c_0n^{-1/5} hopt=c0n−1/5
其中 c 0 = κ 2 − 2 / 5 κ 1 / 5 { ∫ [ f ( 2 ) ( x ) ] 2 d x } − 1 / 5 > 0 c_0=\kappa_2^{-2/5}\kappa^{1/5}\left\{\int [f^{(2)}(x)]^2 dx\right\}^{-1/5}>0 c0=κ2−2/5κ1/5{∫[f(2)(x)]2dx}−1/5>0。

1.1.2 窗宽选择

插入法(plug-in methods)
为了求出在IMSE条件下最有的窗宽，需要确定常数 c 0 c_0 c0中的 ∫ [ f ( 2 ) ( x ) ] 2 d x \int \left[\,f^{(2)}(x)\,\right]^2\,dx ∫[f(2)(x)]2dx。由于 f f f是未知的，所以这个量无法事先知道。如果选择一个 h h h初始的“试验值”(pilot value)，然后将这个值代入 h o p t h_{opt} hopt的计算式求出的优化 h h h，则这种方法称为“插入法”(plug-in methods)。
Silverman(1986)提出假定 f f f是一个以 σ 2 \sigma^2 σ2为方差的正态分布，则其二阶导可确定， ∫ [ f ( 2 ) ( x ) ] 2 d x = 3 8 π σ 5 \int \left[\,f^{(2)}(x)\,\right]^2\,dx=\frac{3}{8\sqrt{\pi}\sigma^5} ∫[f(2)(x)]2dx=8π σ53，代入优化窗宽，可以得到试验窗宽估计
h p i l o t = ( 4 π ) − 1 / 10 [ ( 3 / 8 ) π ] − 1 / 5 σ n − 1 / 5 ≈ 1.06 σ n − 1 / 5 h_{pilot}=(4\pi)^{-1/10}[(3/8)\sqrt{\pi}]^{-1/5}\sigma n^{-1/5}\approx1.06\sigma n^{-1/5} hpilot=(4π)−1/10[(3/8)π ]−1/5σn−1/5≈1.06σn−1/5
用此试验值进一步迭代计算 ∫ [ f ^ ( 2 ) ( x ) ] 2 d x \int \left[\hat f^{(2)}(x)\,\right]^2\,dx ∫[f^(2)(x)]2dx，定出最终的优化结果 h o p t h_{opt} hopt。
Silverman还提出一种更加稳健的分散程度度量，就是用 min ⁡ { σ , q 1 / 4 / 1.34 } \min\{\sigma, q_{1/4}/1.34\} min{σ,q1/4/1.34}来代替 σ \sigma σ，其中 q 1 / 4 q_{1/4} q1/4表示四分位矩。
交错鉴定法
交错鉴定法是一种完全由数据驱动的方法，其核心在于用一部分样本拟合模型来检验另一部分样本的拟合程度。通过不断改变训练集合测试集，来评价模型的好坏。当每次都只留一个样本作为检验对象，其他样本均做训练集时，所得到的估计量称为去一估计量(leave-one-out estimator)。
通过这种方法，我们可以来估计 f ^ \hat{f} f^和 f 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*}
其中第三项和 f ^ \hat{f} f^无关，视为常数
∫ f ( x ) 2 d x = C \int f(x)^2dx=C ∫f(x)2dx=C
第二项采用去一估计量估计，即
∫ f ^ ( x ) f ( x ) d x = E X [ f ^ ( X ) ] = 1 n ∑ i = 1 n f ^ − i ( X i ) + O ( n − 1 / 2 ) \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}) ∫f^(x)f(x)dx=EX[f^(X)]=n1i=1∑nf^−i(Xi)+O(n−1/2)
其中 E x [ ⋅ ] E_x[\cdot] Ex[⋅]是对 x x x求期望，用来区别对观测量 X i X_i Xi求期望。在 X i X_i Xi处的去一估计量 f ^ − i ( X i ) \hat{f}_{-i}(X_i) f^−i(Xi)定义为
f ^ − i ( X i ) = 1 ( n − 1 ) h ∑ j ≠ i n k ( X i − X j h ) \hat{f}_{-i}(X_i)=\frac{1}{(n-1)h}\sum_{j\neq i}^n k\left(\frac{X_i-X_j}{h}\right) f^−i(Xi)=(n−1)h1j=i∑nk(hXi−Xj)
表示用除了 X i X_i Xi这个观测量外的其他观测量来估计 X i X_i Xi处的密度函数。
第一项直接代入 f ^ ( x ) \hat{f}(x) 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*}
其中 k ˉ ( t ) = ∫ k ( x ) k ( t − x ) d x \bar{k}(t)=\int k(x)k(t-x)\,dx kˉ(t)=∫k(x)k(t−x)dx是 k ( ⋅ ) k(\cdot) k(⋅)的重卷积核(two-fold convolution)，一般是两个独立同分布的随机变量之和的分布。可证明， k ˉ ( ⋅ ) \bar{k}(\cdot) kˉ(⋅)也是偶函数。

定理 1.3 总体分布函数为 f ( x ) f(x) f(x)，通过去一核估计交叉检验得到的估计量 f ^ \hat{f} f^的积分平方误差 C V CV CV为
C V f ( h ) = 1 n 2 h ∑ i = 1 n ∑ j = 1 n k ˉ ( X i − X j h ) − 2 n ( n − 1 ) h ∑ i = 1 n ∑ j ≠ i n k ( X i − X j h ) + C CV_f(h)=\frac{1}{n^2h}\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\neq i}^nk\left(\frac{X_i-X_j}{h}\right)+C CVf(h)=n2h1i=1∑nj=1∑nkˉ(hXi−Xj)−n(n−1)h2i=1∑nj=i∑nk(hXi−Xj)+C
其中 k ˉ ( t ) = ∫ k ( x ) k ( t − x ) d x \bar{k}(t)=\int k(x)k(t-x)\,dx kˉ(t)=∫k(x)k(t−x)dx是 k ( ⋅ ) k(\cdot) k(⋅)的重卷积核。

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

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

1.2 单变量累计分布函数

1.2.1 累计分布函数的核估计

为了得到平滑的CDF估计量，我们从核函数出发，将密度分布函数估计进行积分
F ^ ( x ) = ∫ − ∞ x f ^ ( x ) d x = 1 n ∑ i = 1 n G ( x − X i h ) \hat{F}(x)=\int_{-\infty}^x\hat{f}(x)dx=\frac{1}{n}\sum_{i=1}^nG\left(\frac{x-X_i}{h}\right) F^(x)=∫−∞xf^(x)dx=n1i=1∑nG(hx−Xi)
其中 G ( x ) = ∫ − ∞ x k ( x ) d x G(x)=\int_{-\infty}^xk(x)dx G(x)=∫−∞xk(x)dx是核的累计分布函数。其均方误差有下面定理给出

定理 1.4：总体的累计分布函数 F ( x ) F(x) F(x)二阶连续可微，且二阶倒数Holder连续， k ( x ) k(x) k(x)为对称的核函数， G ( x ) = ∫ ∞ x k ( x ) G(x)=\int_\infty^xk(x) G(x)=∫∞xk(x)为核积分函数。则当 n → ∞ n\to\infty n→∞时，
\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 ) F(x) F(x)和核确定 k ( x ) k(x) k(x)。

因此，可以容易的 F ^ \hat{F} 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 = ∫ c i ( x ) d x C_i=\int c_i(x)dx Ci=∫ci(x)dx是和 x x x无关的常数。
首项最小化可以的到优化的核宽度选择
h o p t = [ C 1 4 C 2 ] 1 / 3 n − 1 / 3 h_{opt}=\left[\frac{C_1}{4C_2}\right]^{1/3}n^{-1/3} hopt=[4C2C1]1/3n−1/3
这比密度估计( n − 1 / 5 n^{-1/5} n−1/5)收敛速度要快。

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

1.2.2 窗宽选择

交叉检验法：累计分布函数估计 F ^ ( x ) \hat{F}(x) F^(x)的交叉检验函数定义如下
C V F ( h ) = 1 n ∑ i = 1 n ∫ [ 1 ( X i ≤ x ) − F ^ − i ( x ) ] 2 d x CV_F(h)=\frac{1}{n}\sum_{i=1}^n\int \left[\mathbf{1}(X_i\le x)-\hat{F}_{-i}(x)\right]^2dx CVF(h)=n1i=1∑n∫[1(Xi≤x)−F^−i(x)]2dx
其中 1 \mathbf{1} 1是示性函数， F ^ − i ( x ) = 1 n − 1 ∑ j ≠ i G ( x − X j h ) \hat{F}_{-i}(x)=\frac{1}{n-1}\sum_{j\ne i}G\left(\frac{x-X_j}{h}\right) F^−i(x)=n−11∑j=iG(hx−Xj)为去一核估计量。

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

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

1.3.1 联合分布的核估计

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

s	=		1	2	…	q
X 1 X_1 X1	=	(	X 11 X_{11} X11,	X 12 X_{12} X12,	…,	X 1 q X_{1q} X1q	)
X 2 X_2 X2	=	(	X 21 X_{21} X21,	X 22 X_{22} X22,	…,	X 2 q X_{2q} X2q	)
…	=	(	…,	…,	…,	…	)
X n X_n Xn	=	(	X n 1 X_{n1} Xn1,	X n 2 X_{n2} Xn2,	…,	X n q X_{nq} Xnq	)

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

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

定理 1.5：设三阶梯度存在的 q q q维联合概率密度分布函数 f ( x ) ≡ f ( x 1 , x 2 , … , x q ) f(x)\equiv f(x_1, x_2, \ldots, x_q) f(x)≡f(x1,x2,…,xq)有一组i.i.d.的 n n n个观测值 { X n ∈ R q } \{X_n\in\mathbb{R}^q\} {Xn∈Rq}。核函数 K ( x ) K(x) K(x)为单变量核函数之积。且当 n → ∞ n\to\infty n→∞时，有格子体积宏观无穷小 max ⁡ i h i → 0 \max_{i}h_i\to 0 maxihi→0，微观无穷大 n h 1 h 2 ⋯ h q → ∞ nh_1h_2\cdots h_q\to\infty nh1h2⋯hq→∞。则对于 x ∈ supp ( X ) x\in\text{supp}(X) x∈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*}
其中 κ = ∫ k 2 ( x ) d x \kappa=\int k^2(x)dx κ=∫k2(x)dx， κ 2 = ∫ x 2 k ( x ) d x \kappa_2=\int x^2k(x)dx κ2=∫x2k(x)dx由单变量核函数性质决定。 L L L为核宽度超立方体的对角线长度，而 V V V为超立方体的体积。

渐进正态性讨论
如果 n → ∞ n\to\infty n→∞，格子宏观无穷小 max ⁡ i h i → 0 \max_{i}h_i\to 0 maxihi→0，微观无穷大 n V → ∞ nV\to\infty nV→∞时，并且 n V ∑ s = 1 q h s 6 → 0 nV\sum_{s=1}^qh_s^6\to 0 nV∑s=1qhs6→0，密度估计量具有渐进正态性。
f ^ ( x ) − f ( x ) − b i a s [ f ^ ( x ) ] → N ( 0 , κ q f ( x ) n V ) \hat{f}(x)-f(x)-bias[\hat{f}(x)]\rightarrow\mathcal N(0, \frac{\kappa^qf(x)}{nV}) f^(x)−f(x)−bias[f^(x)]→N(0,nVκqf(x))
即其无偏误误差服从均值为0的正态分布。

1.3.2 窗框选择

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

对于插入法，一般通过 f ^ ( x ) \hat{f}(x) f^(x)的偏误和方法首项进行估计，其中包含了总体分布 f ( x ) 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 , … , h q ) CV_f(h_1,\ldots,h_q) CVf(h1,…,hq)的首项通过下式给出
C V f 0 ( h 1 , h 2 , … , h q ) = ∫ [ ∑ s = 1 q B s ( x ) h s 2 ] 2 d x + κ q n h 1 h 2 ⋯ h q CV_{f_0}(h_1, h_2, \ldots, h_q)=\int\left[\sum_{s=1}^qB_s(x)\,h_s^2\right]^2dx+\frac{\kappa^q}{nh_1h_2\cdots h_q} CVf0(h1,h2,…,hq)=∫[s=1∑qBs(x)hs2]2dx+nh1h2⋯hqκq
其中 B s ( x ) = κ 2 2 ∂ 2 f ( x ) ∂ x s 2 B_s(x)=\frac{\kappa_2}{2}\frac{\partial^2f(x)}{\partial x_s^2} Bs(x)=2κ2∂xs2∂2f(x)， κ = ∫ k 2 ( x ) d x \kappa=\int k^2(x)dx κ=∫k2(x)dx， κ 2 = ∫ x 2 k ( x ) d x \kappa_2=\int x^2k(x)dx κ2=∫x2k(x)dx。
为了分离出样本数 n n n的影响，我们定义 a s = h s n 1 / ( q + 4 ) a_s=h_s\,n^{1/(q+4)} as=hsn1/(q+4)，代换 h s h_s hs得到
C V f 0 ( h 1 , h 2 , … , h q ) = n − 1 / ( q + 4 ) χ ( a 1 , a 2 , … , a q ) CV_{f_0}(h_1, h_2, \ldots, h_q)=n^{-1/(q+4)}\chi(a_1, a_2,\ldots, a_q) CVf0(h1,h2,…,hq)=n−1/(q+4)χ(a1,a2,…,aq)
其中 χ ( a 1 , a 2 , … , a q ) \chi(a_1, a_2,\ldots, a_q) χ(a1,a2,…,aq)适合 n n n无关的常数，定义为
χ ( h 1 , h 2 , … , h q ) = ∫ [ ∑ s = 1 q B s ( x ) a s 2 ] 2 d x + κ q a 1 a 2 ⋯ a q \chi(h_1, h_2, \ldots, h_q)=\int\left[\sum_{s=1}^qB_s(x)\,a_s^2\right]^2dx+\frac{\kappa^q}{a_1a_2\cdots a_q} χ(h1,h2,…,hq)=∫[s=1∑qBs(x)as2]2dx+a1a2⋯aqκq

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

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

1.4 高阶核函数

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

归一化, ∫ k ( x ) d x = 1 \int k(x)dx =1 ∫k(x)dx=1

低阶矩为0, ∫ x l k ( x ) d x = 0 \int x^l\,k(x)dx=0 ∫xlk(x)dx=0， l = 1 , ⋯ , ν − 1 l=1,\cdots, \nu-1 l=1,⋯,ν−1

ν \nu ν阶矩有限, ∫ x ν k ( x ) d x = κ ν ≠ 0 < ∞ \int x^\nu k(x)dx=\kappa_\nu\ne 0<\infty ∫xνk(x)dx=κν=0<∞

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

通常使用的核都属于二阶核函数 ν = 2 \nu=2 ν=2。与二阶核类似，对于总体分布函数 f ( x ) f(x) 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 ν阶核函数， n u ≥ 2 nu\ge 2 nu≥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*}
利用一个高阶和可以同时较少偏误和方法。

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

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

1.5 展望

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

参考资料

[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

核密度估计基础-Part I相关推荐

核密度估计python_核密度估计Kernel Density Estimation(KDE)
在介绍核密度评估Kernel Density Estimation(KDE)之前,先介绍下密度估计的问题.由给定样本集合求解随机变量的分布密度函数问题是概率统计学的基本问题之一.解决这一问题的方法包括 ...
核密度估计python_非参数估计：核密度估计KDE
核密度估计Kernel Density Estimation(KDE)概述密度估计的问题由给定样本集合求解随机变量的分布密度函数问题是概率统计学的基本问题之一.解决这一问题的方法包括参数估计和非参 ...
核密度估计Kernel Density Estimation(KDE)-代码详细解释
在介绍核密度评估Kernel Density Estimation(KDE)之前,先介绍下密度估计的问题.由给定样本集合求解随机变量的分布密度函数问题是概率统计学的基本问题之一.解决这一问题的方法包括 ...
非参数估计：核密度估计KDE
http://blog.csdn.net/pipisorry/article/details/53635895 核密度估计Kernel Density Estimation(KDE)概述密度估计的问 ...
核密度估计与自适应带宽的核密度估计
最近看论文,发现一个很不错的概率密度估计方法.在此小记一下. 先来看看准备知识. 密度估计经常在统计学中作为一种使用有限的样本来估计其概率密度函数的方法. 我们在研究随机变量的过程中,随机变量的概率密 ...
核密度聚类（二）核密度估计、自适应核密度的数学原理
1 分布密度函数给定一个样本集,怎么得到该样本集的分布密度函数,解决这一问题有两个方法: 1.1 参数估计方法简单来讲,即假定样本集符合某一概率分布,然后根据样本集拟合该分布中的参数,例如:似然估 ...
【R语言】-核密度估计图绘制
本期介绍了利用R语言ggplot2包绘制核密度估计图. 核密度估计图(Kernel Density Estimation, KDE),是在概率论中用来估计未知的密度函数,属于非参数检验方法之一,由Ro ...
基于密度的聚类（Density-based clustering）-- 核密度估计（kernel density estimation）
In density-based clustering, clusters are defined as areas of higher density than the remainder of t ...
机器学习算法（二十一）：核密度估计 Kernel Density Estimation(KDE)
目录 1 分布密度函数 1.1 参数估计方法 1.2 非参数估计 2 直方图到核密度估计 2.1 核函数 2.2 带宽的选择 2.2.1 自适应或可变带宽的核密度估计 2.3 多维 1 分布密度函数 ...

核密度估计基础-Part I