核密度估计

4 篇文章 0 订阅
2 篇文章 1 订阅

核密度估计

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

直方图密度估计

基本概念

要说核密度估计就不得不提到直方图估计,直方图通常用来描述数据的频率,可以使我们对数据的分布有一个基础的大概的了解.直方图估计与直方图的区别在于,直方图估计需要对频率进行归一化,使其满足概率密度函数的定义.
考虑一元情况下的直方图密度估计,假定有数据 x 1 , x 2 , ⋯   , x n ∈ [ a , b ) x_1,x_2,\cdots,x_n\in[a,b) x1,x2,,xn[a,b).现在对区间 [ a , b ) [a,b) [a,b)做如下的划分:
a = a 0 < a 1 < a 2 < ⋯ < a k = b a=a_0<a_1<a_2<\cdots<a_k=b a=a0<a1<a2<<ak=b
将任意的小区间定义为 I i = [ a i − 1 , a i ) I_i=[a_{i-1},a_i) Ii=[ai1,ai), i = 1 , 2 , ⋯   , k i=1,2,\cdots,k i=1,2,,k.显然有:
⋃ i = 1 k I i = [ a , b ) 且 I i ∩ I j = ∅ , i ≠ j \bigcup_{i=1}^kI_i=[a,b)且I_i\cap I_j=\varnothing,i\neq j i=1kIi=[a,b)IiIj=,i=j
定义 n i = n_i= ni=# { x i ∈ I i } \left\{ x_i\in I_i \right\} {xiIi}为落在 I i I_i Ii区间中数据的个数.这样就能定义直方图密度估计:
p ^ ( x ) = { n i n ( a i − a i − 1 ) , 当 x ∈ I i ; 0 , 当 x ∉ [ a , b ) , \hat p(x)= \begin{cases} \frac{n_i}{n(a_i-a_{i-1})}, & \text{当} x\in I_i ;\\ 0, & \text{当} x\notin[a,b), \\ \end{cases} p^(x)={n(aiai1)ni,0,xIi;x/[a,b),
在实际操作中,我们经常取相同的区间,即 I i ( i = 1 , 2 , ⋯   , k ) I_i(i=1,2,\cdots,k) Ii(i=1,2,,k)的宽度均为 h h h,这时有:
p ^ ( x ) = { n i n h , 当 x ∈ I i ; 0 , 当 x ∉ [ a , b ) , \hat p(x)= \begin{cases} \frac{n_i}{nh}, & \text{当} x\in I_i ;\\ 0, & \text{当} x\notin[a,b), \\ \end{cases} p^(x)={nhni,0,xIi;x/[a,b),
上式中, h h h既是归一化参数,又表示了每一组的组距,称为窗宽或者带宽.同时我们可以验证 p ^ ( x ) \hat{p}(x) p^(x)是概率密度函数:
∫ a b p ^ ( x ) d x = ∑ i = 1 k ∫ I i n i / ( n h ) d x = ∑ i = 1 k n i / ( n h ) d x = ∑ i = 1 k n i / n = 1 \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}^kn_i/(nh)\mathrm{d}x=\sum_{i=1}^kn_i/n=1 abp^(x)dx=i=1kIini/(nh)dx=i=1kni/(nh)dx=i=1kni/n=1
由于位于同一小区间内的所有点的直方图密度估计都相等,因此直方图对应的分布函数 F ^ h ( x ) \hat{F}_h(x) F^h(x)显然是单调递增的阶梯函数,当 h h h小到每个区间仅能容下一个数据时,直方图的分布函数就是经验分布函数了.举个例子:
在这里插入图片描述
这是4月2日全国各城市累计确诊新冠肺炎人数的直方图密度估计.

核密度估计

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

定义

给定一组数据 x 1 , x 2 , ⋯   , x n x_1,x_2,\cdots,x_n x1,x2,,xn是取自于一个未知的连续分布 p ( x ) p(x) p(x),任意点 x x x处的一种核密度估计定义为:
p ^ ( x ) = 1 n h ∑ i = 1 n K ( x − x i h ) \hat{p}(x)=\frac{1}{nh}\sum_{i=1}^nK \left ( \frac{x-x_i}{h} \right ) p^(x)=nh1i=1nK(hxxi)
其中 K ( ⋅ ) K(\cdot) K()称为核函数(kernel function). p ^ ( x ) \hat{p}(x) p^(x)作为密度函数就需要保证非负和积分值为1,实际上只需要令核函数满足密度函数的定义即可:
K ( x ) ≥ 0 , ∫ K ( x ) d x = 1 K(x)\ge0, \int K(x)\mathrm{d}x=1 K(x)0,K(x)dx=1
这时:
∫ p ^ ( x ) d x = ∫ 1 n ∑ i = 1 n 1 h K ( x − x i h ) d x = 1 n ∑ i = 1 n ∫ 1 h K ( x − x i h ) d x = 1 n ∑ i = 1 n ∫ K ( u ) d u = 1 n ⋅ n = 1 ( 其 中 u = x − x i h ) . \begin{aligned} \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{aligned} p^(x)dx=n1i=1nh1K(hxxi)dx=n1i=1nh1K(hxxi)dx=n1i=1nK(u)du=n1n=1(u=hxxi).
上式验证了 p ^ ( x ) \hat{p}(x) p^(x)符合概率密度函数的定义,核密度估计中一个重要的部分就是核函数,常用的核函数有:

核函数名称核函数 K ( u ) K(u) K(u)
Parzen窗(均匀分布) 1 2 I ( ∣ u ∣ ≤ 1 ) \frac{1}{2}I(\lvert u\rvert \leq1) 21I(u1)
三角(Triangle) ( 1 − ∣ u ∣ ) I ( ∣ u ∣ ≤ 1 ) (1-\lvert u\rvert)I(\lvert u\rvert \leq1) (1u)I(u1)
Epanechikov 3 4 ( 1 − u 2 ) I ( ∣ u ∣ ≤ 1 ) \frac{3}{4}(1-u^2)I(\lvert u\rvert \leq1) 43(1u2)I(u1)
四次(Quartic) 15 16 ( 1 − u 2 ) I ( ∣ u ∣ ≤ 1 ) \frac{15}{16}(1-u^2)I(\lvert u\rvert \leq1) 1615(1u2)I(u1)
三权(Triweight) 35 32 ( 1 − u 2 ) 3 I ( ∣ u ∣ ≤ 1 ) \frac{35}{32}(1-u^2)^3I(\lvert u\rvert \leq1) 3235(1u2)3I(u1)
高斯(Gauss) 1 2 π exp ⁡ ( − 1 2 u 2 ) \frac{1}{\sqrt{2\pi}}\exp(-\frac{1}{2}u^2) 2π 1exp(21u2)
余弦(Cosinus) π 4 c o s ( π 2 u ) I ( ∣ u ∣ ≤ 1 ) \frac{\pi}{4}cos(\frac{\pi}{2}u)I(\lvert u\rvert \leq1) 4πcos(2πu)I(u1)
指数(Exponent) exp ⁡ ( ∣ u ∣ ) \exp(\lvert u \rvert) exp(u)

在这里插入图片描述
可以看出,相同窗宽时不同的函数对数据密度曲线的估计情况差别不大,但如果调整窗宽后:
在这里插入图片描述
同一核函数时,不断扩大核函数窗宽,数据的密度估计曲线越发接近于正态曲线,可见窗宽的选择相比较核函数的选择更为重要。那么如何选择合适的窗宽?理论上是从密度估计与真实密度之间的误差开始的。

窗宽选择

确定窗宽 h h 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=nh1K22+4h4μ2(K)2p22+o(nh1)+o(h4)
其中 μ 2 ( K ) = ∫ t 2 K ( t ) d t \mu_2(K)=\int t^2K(t)\mathrm{d}t μ2(K)=t2K(t)dt, t = ( x − x i ) / h t=(x-x_i)/h t=(xxi)/h.当 h → ∞ , n h → ∞ h\to\infty,nh\to\infty h,nh时,可以得到渐进的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)}=nh1K22+4h4μ2(K)2p22
因此最优窗宽为:
h o p t = ( ∥ K ∥ 2 2 ∥ p ′ ′ ∥ 2 2 μ 2 ( K ) 2 n ) 1 / 5 ∽ n 1 / 5 h_{opt}=\left( \frac{\lVert K\rVert_2^2}{\lVert p''\rVert_2^2\mu_2(K)^2n} \right)^{1/5} \backsim n^{1/5} hopt=(p22μ2(K)2nK22)1/5n1/5
通过上式可以知道,对于不同的核函数,我们可以得到相应的最优窗宽,例如高斯核函数时,可以得到 μ 2 = 1 , ∫ K ( u ) 2 d u = ∫ 1 / 2 π exp ⁡ ( − u 2 ) d u = π − 1 / 2 \mu_2=1,\int K(u)^2\mathrm{d}u=\int1/2\pi\exp(-u^2)\mathrm{d}u=\pi^{-1/2} μ2=1,K(u)2du=1/2πexp(u2)du=π1/2,那么最优窗宽就是 h o p t = 1.06 σ n − 1 / 5 h_{opt}=1.06\sigma n^{-1/5} hopt=1.06σn1/5.

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Infinity343

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值