我们可以举一个极端的例子:在R中输入:
plot(density(rep(0, 1000)))
可以看到它得到了正态分布的曲线,但实际上呢?从数据上判断,它更有可能是一个退化的单点分布。
但是这并不意味着核密度估计是不可取的,至少他可以解决许多模拟中存在的异方差问题。比如说我们要估计一下下面的一组数据:
set.seed(10)
dat<-c(rgamma(300,shape=2,scale=2),rgamma(100,shape=10,scale=2))
可以看出它是由300个服从gamma(2,2)与100个gamma(10,2)的随机数构成的,他用参数统计的办法是没有办法得到一个好的估计的。那么我们尝试使用核密度估计:
plot(density(dat),ylim=c(0,0.2))
将利用正态核密度与标准密度函数作对比
dfn<-function(x,a,alpha1,alpha2,theta){
a*dgamma(x,shape=alpha1,scale=theta)+(1-a)*dgamma(x,shape=alpha2,scale=theta)}
pfn<-function(x,a,alpha1,alpha2,theta){
a*pgamma(x,shape=alpha1,scale=theta)+(1-a)*pgamma(x,shape=alpha2,scale=theta)}
curve(dfn(x,0.75,2,10,2),add=T,col="red")
如果换用gamma分布的核效果无疑会更好,但是遗憾的是r中并没有提供那么多的核供我们挑选(其实我们知道核的选择远没有窗宽的选择来得重要),所以也无需介怀。
R中提供的核:kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight","cosine", "optcosine")。
我们先来看看窗宽的选择对核密度估计的影响:
dfn1<-function(x){
0.5*dnorm(x,3,1)+0.5*dnorm(x,-3,1)}
par(mfrow=c(2,2))
curve(dfn1(x),from=-6,to=6)
data<-c(rnorm(200,3,1),rnorm(200,-3,1))
plot(density(data,bw=8))
plot(density(data,bw=0.8))
plot(density(data,bw=0.08))
窗宽究竟该如何选择呢?
我们这里不加证明的给出最佳窗宽选择公式:
(这个基于积分均方误差最小的角度得到的)
这里介绍两个可操作的窗宽估计办法:(这两种方法都比较容易导致过分光滑)
1、 Silverman大拇指法则
这里使用R(phi’’)/sigma^5估计R(f’’),phi代表标准正态密度函数,得到h的表达式:
h=(4/(3n))^(*1/5)*sigma
2、 极大光滑原则
h=3*(R(K)/(35n))^(1/5)*sigma
当然也有比较麻烦的窗宽估计办法,比如缺一交叉验证,插入法等,可以参阅《computational statistics》一书
我们用上面的两种办法得到的窗宽是多少,他的核密度估计效果好吗?
我们还是以上面的混合正态数据为例来看看效果。
使用大拇指法则,将数据n=400,sigma=3.030658,带入公式,h=0.9685291
使用极大光滑原则,假设K为正态核,R(K)=1/(sqrt(2*pi)),h=1.121023
可以看出他们都比我们认为的h=0.8要大一些,作图如下:
plot(density(data,bw=0.9685))
plot(density(data,bw=1.1210))