有边界区间上的核密度估计

一、一个例子

核密度估计应该是大家常用的一种非参数密度估计方法,从某种程度上来说它的性质比直方图更好,可以替代直方图来展示数据的密度分布。但是相信大家会经常遇到一个问题,那就是有些数据是严格大于或等于零的,在这种情况下,零附近的密度估计往往会出现不理想的情况。下面以一个指数分布的模拟数据为例(样本量为1000),R程序代码为:

set.seed(123);
x=rexp(1000,1);
plot(density(x,kernel="epanechnikov"),ylim=c(0,1.2));
lines(seq(0,8,by=0.02),dexp(seq(0,8,by=0.02),1),col="blue");
abline(v=0,col="red");

可以看出,理论上应该单调递减的密度函数在0附近有明显的“陡坡”,而且不应该有密度的小于零的区域也有着正的估计值。当样本量增大时,这种现象也不会得到明显好转,下图是将样本量改为10000时的情形。

set.seed(123);
x=rexp(10000,1);
plot(density(x,kernel="epanechnikov"),ylim=c(0,1.2));
lines(seq(0,8,by=0.02),dexp(seq(0,8,by=0.02),1),col="blue");
abline(v=0,col="red");

我们也许从平时看的书中了解到,当样本量趋于无穷时,核密度估计值将是收敛到真实的密度函数的,但我们可能不会特意去研究那些结论成立的条件。以上这两个简单的例子似乎给了我们一个直观的感觉,那就是当真实密度函数的支撑集(函数f(x)的支撑集指的是使得f(x)≠0的x的集合)有边界时,核密度估计值可能会出现一些不理想的情况。下面就简单地给出一些理论的结果。

二、理论分析

在一些必要的条件下(真实的密度函数f二阶导绝对连续,三阶导平方可积),核密度估计值$\hat{f}(x)$的偏差有表达式$Bias[\hat{f}(x)]=\frac{h^2\sigma_k^2f”(x)}{2}+O(h^4)$,其中h是带宽,$\sigma_k^2=\int u^2k(u)du$,k(u)是支集为[-1,1]的核函数(即在[-1,1]上有值,其余的地方取零)。可以看出这个偏差是随着带宽h的减小以$h^2$的速度趋于零的。

而假设密度函数以0为边界,那么上述表达式将不再成立,而是代之以
$E[\hat{f}_k(x)]=a_0(x)f(x)-ha_1(x)f'(x)+O(h^2)$
其中$a_i(x)=\int_{-1}^{x/h}u^ik(u)du$。可以看出,当$x \ge h$时,$a_0(x)=1$,$a_1(x)=0$,此时的偏差跟之前的那个表达式没有区别;但当$0 \le x<h$时,$a_0(x)$和$a_1(x)$都是非零的,于是偏差总是存在。

也许你会提议说,将估计值除以$a_0(x)$,偏差就可以减小了吧?的确,这样是一种改进的办法,但也要注意到,此时h的一次项不会消除,也就是说原来$h^2$的衰减速度放慢到了h,从效率上说相对于理想的情况是大打了折扣。

这时候一个巧妙的办法是,用另外一个核函数l(x)对f也做一次估计,那么就有
$E[\hat{f}_l(x)]=b_0(x)f(x)-hb_1(x)f'(x)+O(h^2)$
其中的$b_0$和$b_1$意义类似,只不过是针对l(x)的。

对以上两个式子进行线性组合,则会有
$b_1(x)*E[\hat{f}_k(x)]-a_1(x)*E[\hat{f}_l(x)]=[b_1(x)a_0(x)-a_1(x)b_0(x)]f(x)+O(h^2)$
然后把f(x)的系数移到等式左边,O(h)项的偏差就神奇地消失了。

通过观察核密度估计的表达式,我们可以将上面这个过程等价地认为是对f(x)用了一个新的核函数进行估计,这个新的核函数是
$p(x)=\frac{b_1(x)k(x)-a_1(x)l(x)}{b_1(x)a_0(x)-a_1(x)b_0(x)}$

特别地,如果将l(x)取为x*k(x),那么p(x)将有一个简单的形式
$p(x)=\frac{(a_2(x)-a_1(x)x)k(x)}{a_0(x)a_2(x)-a_1^2(x)}$

当$x \ge h$时,这个新的核函数p(x)就是k(x),而当$x \ge h$时(也就是在边界),它会对最初的核函数进行调整。当$x<0$时,不管算出来的估计值是多少,都只需将密度的估计值取为0即可。

三、程序实现

下面这段程序是对本文的第一幅图进行“整容”,代码及效果图如下:

k=function(x) 3/4*(1-x^2)*(abs(x)<=1);
a0=function(u,h)
{
	lb=-1;
	ub=pmin(u/h,1);
	0.75*(ub-lb)-0.25*(ub^3-lb^3);
}
a1=function(u,h)
{
	lb=-1;
	ub=pmin(u/h,1);
	3/8*(ub^2-lb^2)-3/16*(ub^4-lb^4);
}
a2=function(u,h)
{
	lb=-1;
	ub=pmin(u/h,1);
	0.25*(ub^3-lb^3)-0.15*(ub^5-lb^5);
}
kernel.new=function(x,u,h)
{
	k(x)*(a2(u,h)-a1(u,h)*x)/(a0(u,h)*a2(u,h)-a1(u,h)^2);
}
den.est=function(u,ui,h)
{
	sapply(u,function(u) ifelse(u<0,0,mean(kernel.new((u-ui)/h,u,h))/h));
}
set.seed(123);
dat=rexp(1000,1);
x=seq(0,8,by=0.02);
y=den.est(x,dat,2*bw.nrd0(dat));
plot(x,y,type="l",ylim=c(0,1.2),main="Corrected Kernel");
lines(x,dexp(x,1),col="red");

从中可以看出,边界的偏差问题已经得到了很好的改进。

如果真实的密度函数的支集不是[0,+∞],而是某一个闭区间[m,n],那么偏差修正的过程与上面类似,只不过是要将$a_i(x)$定义为$a_i(x)=\int_{(x-n)/h}^{(x-m)/h}u^ik(u)du$。在编程序的时候,也只需把积分的上下限进行相应的调整即可。

四、参考文献

Jeffrey S. Simonoff, 1998. Smoothing Methods in Statistics. Springer-Verlag

相关链接:http://pages.stern.nyu.edu/~jsimonof/SmoothMeth/

有边界区间上的核密度估计》有22个想法

  1. 这个办法看起来挺理论的,学习了!
    不过有个更简单的办法,直接对数据LOG后再在LOG后的数据上建密度函数,再转化过来,不就好了吗。–sliverman 84年那本书上是这么说的

    1. 多谢提醒,我用对数变换的方法试了一下,公式是
      [latex]\hat{f}(x)=\frac{1}{nxh}\sum(\frac{\ln x-\ln x_i}{h})[/latex]
      程序如下:
      k=function(x) 3/4*(1-x^2)*(abs(x)< =1); den.est=function(u,ui,h) { sapply(u,function(u) ifelse(u<0,0, mean(k((log(u)-log(ui))/h))/(u*h))); } set.seed(123); dat=rexp(1000,1); x=seq(0.01,8,by=0.05); y=den.est(x,dat,2*bw.nrd0(log(dat))); plot(x,y,type="l",ylim=c(0,1.2),main="Log Transformed"); lines(x,dexp(x,1),col="red");

      结果发现在很靠近0的时候出现了一个诡异的下降(可能是因为估计式的分母中有x,在0点无定义)。

      1. 直觉感觉:在对数变换之后,原始数据趋向于-INF的分布过于稀疏,nh不趋于正无穷。如果要使对数变换后的数据在-inf的位置上nh趋于正无穷,那么原数据在[0,1]区间上就不应该是均匀分布,在趋于0的位置上的密度应该呈指数增加。但是这个在实际中几乎不能做到,因此对数变换不能消除在0点位置上的一阶偏差。所以原文的方法是更好的处理方法。

      2. 我稍微折腾了一下,基本上就是这个原因,推导在这里。针对上面那个例子,把带宽调大些可以减轻那种现象,相当于是减小方差;但h一大,偏差又大了。如果从MSE方面考虑,还是需要做出一些权衡的。正文中的方法有个缺点是不能保证密度函数的积分是1,不知道是不是有人继续在这方面加以改进。

      3. 恩,我应该更正一下,在0点附近的密度应该成倒数速率增加,因为d(logx)=1/xdx,这也解释了你方差的第一项。
        对于积分为1的问题,可以考虑标准化,但不知能否证明标准化的分母的MSE是否为分子MSE的高阶无穷小。

      1. 这种变换的方法就是把有边界的区间变换到R上吧,比如如果X是[-1,+∞]上的,就令Y=g(X)=log(X+1),然后[latex]f_X(x)=f_Y(g(x))g'(x)[/latex]。

  2. 理论虽然如此,但实际操作中a(x),b(x)都怎么选取?能说说思路吗?
    还有p(x)的公式(倒数第二个)是怎么由倒数第三个导出的?

    1. a_i(x)不是选出来的,它是对核函数进行积分的结果,积分的上下限与区间的边界有关,参见文中的定义式。只要核函数定了,区间定了,a_i(x)就定了。
      把p(x)当作是一个核函数,写出密度估计的表达式,变换后会发现与倒数第三个式子相等。简单地来说,对不同核函数的密度估计进行线性组合,相当于对核函数进行线性组合然后估计密度,前提是保证核函数和密度函数积分为1之类的。

  3. 您好,我们在用SAS编程求核密度的时候也遇到这个问题,在0附近的核密度都大于1,请问怎么解决这个问题呢?

    1. 为什么总是有人认为密度函数值一定要小于1呢?只有累积分布函数(CDF)取值是在[0, 1]上,密度函数没有小于1这种范围。

      1. 密度不是概率,不要求在[0,1]之间。你可以算一下正态分布N(0,0.000001)在0点的密度值。

      2. 是的。谢谢三楼。对于离散的分布列(在连续成为密度函数)的值是在[0,1](由于此时的分布列即是概率值);连续的密度函数值不一定在[0,1]之间,因为它只是函数值不是概率(三楼);例如三楼给的dnorm(0,0,000001)=398942.3;均匀分布[0,0.01]密度函数值全是100!

  4. 文章提供的方法很有实际意义。现实中的很多分布原本就有区间限制。
    我现在就碰到一个这样的问题。根据住院费用的一定赔付数据求出住院赔付的核密度估计。这里的费用明显是非负的。进而需要求解赔付的期望和方差。这里需要用到积分。
    但我遇到的问题是,因为要求根据病人男女不同,医院等级不同等,需要进一步求出每一个细分赔付的核密度估计,并求期望和方差。 根据文章中给出的核密度函数 den.est=function(u,ui,h) 。只需要求 den.est(u,ui,h)*u 在0-Inf上的积分 就可求解期望。 但是这样作用在三个数上的自定义函数,怎样运用tapply 和by 之类的函数实现快捷地对dat在其他条件限制下的核密度函数?一般运用tapply都只是求出了一个值,但是这里需要求出一个函数。如果每次人工将dat按条件挑出来,费时且当限制条件种类很多(如医院等级很多),容易出错。 不知道有没有什么好的方法?
    边学R边用,对R了解不是很深。

  5. I like this boundary bias reduction trick, and it reminds me the somewhat remotely related generalized jackknife bias correction.

发表评论

邮箱地址不会被公开。 必填项已用*标注