核密度估计

核密度估计

核密度估计属于非参数估计方法。概率分布是统计推断的核心问题,一旦给出联合概率密度,就能够回答变量子集之间的所有问题。可以说,参数统计的核心内容就是对密度的估计。通常在实际问题中,很多数据的分布无从得知,只能通过假设确定,而核密度估计正是为了适应这一问题而产生的。

直方图密度估计

基本概念

要说核密度估计就不得不提到直方图估计,直方图通常用来描述数据的频率,可以使我们对数据的分布有一个基础的大概的了解.直方图估计与直方图的区别在于,直方图估计需要对频率进行归一化,使其满足概率密度函数的定义.
考虑一元情况下的直方图密度估计,假定有数据$x_1,x_2,\cdots ,x_n\in [a,b)$.现在对区间$[a,b)$做如下的划分:
$$a=a_0<a_1<a_2<\cdots <a_k=b$$
将任意的小区间定义为$I_i=[a_{i-1},a_i)$,$i=1,2,\cdots ,k$.显然有:
$$\bigcup _{i=1}^k{I}_i=[a,b)且I_i\cap I_j=\varnothing ,i\ne j$$
定义$n_i=$# { x i ∈ I i } \left\{ x_i\in I_i \right\} {xi​∈Ii​}为落在$I_i$区间中数据的个数.这样就能定义直方图密度估计:
$$\hat{p}(x)=\{\begin{array}{ll}\frac{n_i}{n(a_i-a_{i-1})},& \text{当}x\in I_i;\\ 0,& \text{当}x\notin [a,b),\end{array}$$
在实际操作中,我们经常取相同的区间,即$I_i(i=1,2,\cdots ,k)$的宽度均为$h$,这时有：
$$\hat{p}(x)=\{\begin{array}{ll}\frac{n_i}{nh},& \text{当}x\in I_i;\\ 0,& \text{当}x\notin [a,b),\end{array}$$
上式中,$h$既是归一化参数,又表示了每一组的组距,称为窗宽或者带宽.同时我们可以验证$\hat{p}(x)$是概率密度函数:
$${\int }_a^b\hat{p}(x)\mathrm{d}x=\sum _{i=1}^k{\int }_{I_i}{n}_i/(nh)\mathrm{d}x=\sum _{i=1}^k{n}_i/(nh)\mathrm{d}x=\sum _{i=1}^k{n}_i/n=1$$
由于位于同一小区间内的所有点的直方图密度估计都相等,因此直方图对应的分布函数${\hat{F}}_h(x)$显然是单调递增的阶梯函数,当$h$小到每个区间仅能容下一个数据时,直方图的分布函数就是经验分布函数了.举个例子:

这是4月2日全国各城市累计确诊新冠肺炎人数的直方图密度估计.

核密度估计

上节提到的直方图密度估计存在的问题是不光滑,核密度估计方法能够很好的解决这个问题.

定义

给定一组数据$x_1,x_2,\cdots ,x_n$是取自于一个未知的连续分布$p(x)$,任意点$x$处的一种核密度估计定义为:
$$\hat{p}(x)=\frac{1}{nh}\sum _{i=1}^{n}K\left(\frac{x-x_i}{h}\right)$$
其中$K(\cdot )$称为核函数(kernel function).$\hat{p}(x)$作为密度函数就需要保证非负和积分值为1,实际上只需要令核函数满足密度函数的定义即可:
$$K(x)\ge 0,\int K(x)\mathrm{d}x=1$$
这时:
$$\begin{array}{rl}\int \hat{p}(x)\mathrm{d}x& =\int \frac{1}{n}\sum _{i=1}^n\frac{1}{h}K\left(\frac{x-x_i}{h}\right)\mathrm{d}x=\frac{1}{n}\sum _{i=1}^n\int \frac{1}{h}K\left(\frac{x-x_i}{h}\right)\mathrm{d}x\\ & =\frac{1}{n}\sum _{i=1}^n\int K(u)\mathrm{d}u=\frac{1}{n}\cdot n=1\left(其中u=\frac{x-x_i}{h}\right).\end{array}$$
上式验证了$\hat{p}(x)$符合概率密度函数的定义,核密度估计中一个重要的部分就是核函数,常用的核函数有:

核函数名称	核函数$K(u)$
Parzen窗(均匀分布)	$\frac{1}{2}I(|u|\le 1)$
三角(Triangle)	$(1-|u|)I(|u|\le 1)$
Epanechikov	$\frac{3}{4}(1-u^2)I(|u|\le 1)$
四次(Quartic)	$\frac{15}{16}(1-u^2)I(|u|\le 1)$
三权(Triweight)	$\frac{35}{32}(1-u^2{)}^{3}I(|u|\le 1)$
高斯(Gauss)	$\frac{1}{\sqrt{2\pi }}\exp (-\frac{1}{2}{u}^2)$
余弦(Cosinus)	$\frac{\pi }{4}cos(\frac{\pi }{2}u)I(|u|\le 1)$
指数(Exponent)	$\exp (|u|)$

可以看出,相同窗宽时不同的函数对数据密度曲线的估计情况差别不大,但如果调整窗宽后:

同一核函数时,不断扩大核函数窗宽,数据的密度估计曲线越发接近于正态曲线,可见窗宽的选择相比较核函数的选择更为重要。那么如何选择合适的窗宽？理论上是从密度估计与真实密度之间的误差开始的。

窗宽选择

确定窗宽$h$使得估计的积分均方误差(mean integral square error,MISE)达到最小:
M I S E { p ^ h ( x ) } = E [ ∫ − ∞ ∞ { p h ^ ( x ) − p ( x ) } 2 d x ] = ∫ − ∞ ∞ M S E [ p h ^ ( x ) ] d x = 1 n h ∥ K ∥ 2 2 + h 4 4 μ 2 ( K ) 2 ∥ p ′ ′ ∥ 2 2 + o ( 1 n h ) + o ( h 4 ) \begin{aligned} \mathrm{MISE}\{\hat{p}_h(x) \} &= E \left[ \int_{-\infty}^\infty \{ \hat{p_h}(x)-p(x)\}^2\mathrm{d}x \right] \\ &= \int_{-\infty}^\infty \mathrm{MSE}[\hat{p_h}(x)] \mathrm{d}x \\ &= \frac{1}{nh}\lVert K \rVert_2^2+\frac{h^4}{4}\mu_2(K)^2 \lVert p''\rVert_2^2+o(\frac{1}{nh}) +o(h^4) \end{aligned} MISE{p^​h​(x)}​=E[∫−∞∞​{ph​^​(x)−p(x)}2dx]=∫−∞∞​MSE[ph​^​(x)]dx=nh1​∥K∥22​+4h4​μ2​(K)2∥p′′∥22​+o(nh1​)+o(h4)​
其中${\mu }_2(K)=\int t^{2}K(t)\mathrm{d}t$,$t=(x-x_i)/h$.当$h\to \infty ,nh\to \infty$时,可以得到渐进的MISE(AMISE)
A M I S E = { p ^ ( x ) } = 1 n h ∥ K ∥ 2 2 + h 4 4 μ 2 ( K ) 2 ∥ p ′ ′ ∥ 2 2 \mathrm{AMISE}=\{ \hat{p}(x)\}=\frac{1}{nh}\lVert K\rVert_2^2+\frac{h^4}{4}\mu_2(K)^2\lVert p''\rVert_2^2 AMISE={p^​(x)}=nh1​∥K∥22​+4h4​μ2​(K)2∥p′′∥22​
因此最优窗宽为:
$$h_{opt}={\left(\frac{\Vert K{\Vert }_2^2}{\Vert p^{\prime \prime }{\Vert }_2^2{\mu }_2(K{)}^{2}n}\right)}^{1/5}\backsim n^{1/5}$$
通过上式可以知道,对于不同的核函数,我们可以得到相应的最优窗宽,例如高斯核函数时,可以得到${\mu }_2=1,\int K(u{)}^2\mathrm{d}u=\int 1/2\pi \exp (-u^2)\mathrm{d}u={\pi }^{-1/2}$,那么最优窗宽就是$h_{opt}=1.06\sigma n^{-1/5}$.