R语言与非参数统计核密度估计 r语言密度估计曲线

R语言与非参数统计(核密度估计)

分类:R语言数理统计

R语言与非参数统计(核密度估计)

核密度估计是在概率论中用来估计未知的密度函数,属于非参数检验方法之一,由Rosenblatt (1955)和EmanuelParzen(1962)提出,又名Parzen窗(Parzen window)。

假设我们有n个数X1-Xn,我们要计算某一个数X的概率密度有多大。核密度估计的方法是这样的:

其中K为核密度函数,h为设定的窗宽。

核密度估计的原理其实是很简单的。在我们对某一事物的概率分布的情况下。如果某一个数在观察中出现了,我们可以认为这个数的概率密度很大,和这个数比较近的数的概率密度也会比较大,而那些离这个数远的数的概率密度会比较小。基于这种想法,针对观察中的第一个数,我们都可以f(x-xi)去拟合我们想象中的那个远小近大概率密度。当然其实也可以用其他对称的函数。针对每一个观察中出现的数拟合出多个概率密度分布函数之后,取平均。如果某些数是比较重要,某些数反之,则可以取加权平均。

但是核密度的估计并不是,也不能够找到真正的分布函数。我们可以举一个极端的例子:在R中输入:

[plain]viewplaincopyprint?
  1. plot(density(rep(0,1000)))

可以看到它得到了正态分布的曲线,但实际上呢?从数据上判断,它更有可能是一个退化的单点分布。

但是这并不意味着核密度估计是不可取的,至少他可以解决许多模拟中存在的异方差问题。比如说我们要估计一下下面的一组数据:

[plain]viewplaincopyprint?
  1. set.seed(10)
  2. dat<-c(rgamma(300,shape=2,scale=2),rgamma(100,shape=10,scale=2))


R语言与非参数统计(核密度估计) r语言密度估计曲线
可以看出它是由300个服从gamma(2,2)与100个gamma(10,2)的随机数构成的,他用参数统计的办法是没有办法得到一个好的估计的。那么我们尝试使用核密度估计:

[plain]viewplaincopyprint?
  1. plot(density(dat),ylim=c(0,0.2))

将利用正态核密度与标准密度函数作对比[plain]viewplaincopyprint?
  1. dfn<-function(x,a,alpha1,alpha2,theta){
  2. a*dgamma(x,shape=alpha1,scale=theta)+(1-a)*dgamma(x,shape=alpha2,scale=theta)}
  3. pfn<-function(x,a,alpha1,alpha2,theta){
  4. a*pgamma(x,shape=alpha1,scale=theta)+(1-a)*pgamma(x,shape=alpha2,scale=theta)}
  5. curve(dfn(x,0.75,2,10,2),add=T,col="red")


得到下图:

(红色的曲线为真实密度曲线)

可以看出核密度与真实密度相比,得到大致的估计是不成问题的。至少趋势是得到了的。如果换用gamma分布的核效果无疑会更好,但是遗憾的是r中并没有提供那么多的核供我们挑选(其实我们知道核的选择远没有窗宽的选择来得重要),所以也无需介怀。

R中提供的核:kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight","cosine", "optcosine")。
我们先来看看窗宽的选择对核密度估计的影响:
[plain] view plaincopyprint?
  1. dfn1<-function(x){
  2. 0.5*dnorm(x,3,1)+0.5*dnorm(x,-3,1)}
  3. par(mfrow=c(2,2))
  4. curve(dfn1(x),from=-6,to=6)
  5. data<-c(rnorm(200,3,1),rnorm(200,-3,1))
  6. plot(density(data,bw=8))
  7. plot(density(data,bw=0.8))
  8. plot(density(data,bw=0.08))

得到下图,我们可以清楚的看到带宽为0.8恰好合适,其余的不是拟合不足便是过拟合。

窗宽究竟该如何选择呢?

我们这里不加证明的给出最佳窗宽选择公式:

(这个基于积分均方误差最小的角度得到的)

这里介绍两个可操作的窗宽估计办法:(这两种方法都比较容易导致过分光滑)

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要大一些,作图如下:

[plain]viewplaincopyprint?
  1. plot(density(data,bw=0.9685))
  2. plot(density(data,bw=1.1210))


由我们给出的

以Gauss核为例做核密度估计

用Gauss核做核密度估计的R程序如下(还是使用我们的混合正态密度的例子):

[plain]viewplaincopyprint?
  1. ker.density=function(x,h){
  2. x=sort(x)
  3. n=length(x);s=0;t=0;y=0
  4. for(iin2:n)
  5. s[i]=0
  6. for(iin1:n){
  7. for(jin1:n)
  8. s[i]=s[i]+exp(-((x[i]-x[j])^2)/(2*h*h))
  9. t[i]=s[i]
  10. }
  11. for(iin1:n)
  12. y[i]=t[i]/(n*h*sqrt(2*pi))
  13. z=complex(re=x,im=y)
  14. hist(x,freq=FALSE)
  15. lines(z)
  16. }
  17. ker.density(data,0.8)


作图如下:

最后说一个R的内置函数density()。其实我觉得如果不是为了简要介绍核密度估计的一些常识我们完全可以只学会这个函数

先看看函数的基本用法:

density(x, ...)
## Default S3 method:
density(x, bw = "nrd0", adjust = 1,
 kernel = c("gaussian", "epanechnikov", "rectangular",
 "triangular", "biweight",
 "cosine", "optcosine"),
 weights = NULL, window = kernel, width,
 give.Rkern = FALSE,
 n = 512, from, to, cut = 3, na.rm = FALSE, ...)

对重要参数做出较为详细的说明:

X:我们要进行核密度估计的数据

Bw:窗宽,这里可以由我们自己制定,也可以使用默认的办法nrd0: Bandwidth selectors for Gaussiankernels。我们还可以使用bw.SJ(x,nb = 1000, lower = 0.1 * hmax, upper =hmax,method = c("ste","dpi"), tol = 0.1 * lower),这里的method=”dpi”就是前面提到过的插入法,”ste”代表solve-the-equationplug-in,也是插入法的改进

Kernel:核的选择

Weights:对比较重要的数据采取加权处理

对于上述混合正态数据data,有

> density(data)

Call:

density.default(x= data)

Data: data (400obs.); Bandwidth 'bw' = 0.8229

xy

Min.:-7.5040Min.:0.0000191

1stQu.:-3.50761st Qu.:0.0064919

Median :0.4889 Median:0.0438924

Mean:0.4889Mean:0.0624940

3rdQu.:4.4853 3rdQu.:0.1172919

Max.:8.4817Max.:0.1615015

知道带宽:h=0.8229(采取正态密度核)那么带入密度估计式就可以写出密度估计函数。

最后以faithful数据集为例说明density的用法:

R数据集faithful是oldfaithful火山爆发的数据,其中“eruption”是火山爆发的持续时间,waiting是时间间隔

对数据“eruption”做核密度估计

R程序:

[plain]viewplaincopyprint?
  1. data(faithful)
  2. A<-faithful
  3. x<-A[,"eruptions"]
  4. density(x)
  5. plot(density(x))


知道h= 0.3348

作图:

关于核密度估计R中还有不少函数包提供了大量的支持:

可以研读一下如下几个包,也可以自己编程去实现

ks Kernel smoothing

Kendall

Kendall rankcorrelation and Mann-Kendall trend test

KernSmooth

Functions forkernel smoothing for Wand & Jones (1995)

kappalab

Non-additivemeasure and integral manipulation functions

kerfdr

semi-parametrickernel-based approach to local fdr estimations

kernlab

Kernel MethodsLab

  

爱华网本文地址 » http://www.413yy.cn/a/25101013/184352.html

更多阅读

spss教程:两独立样本的非参数检验

spss教程:两独立样本的非参数检验——简介在对总体分布不了解的情况下,通过对两组独立样本的分析来推断样本来自的两个总体分布是否存在显著性差异。Spss提供多种两独立样本的非参数检验方法。其中包括曼-惠特尼U检验、K-S检验、W-W游

郭德纲与“非诚勿扰”之类比 郭德纲上非诚勿扰

郭德纲与“非诚勿扰”之类比“郭德纲徒弟李鹤彪打伤记者”引发的连锁反应还未停,继而李菁、何云伟前日宣布退出德云社后,暂停所有小剧场演出,进行内部自查。书店接上级通知,所有有关郭德纲、德云社的书籍、音像制品全下架。

惯性离心力和向心力的是与非 健康传言是与非

惯性离心力和向心力的是与非作者在《惯性力是一个真实的力》一文中论述了匀加速直线运动的惯性力问题,也就顺理成章地要研究匀速圆周运动的惯性离心力问题,因为那还有一个应该正本清源还其本来面目的被虚拟的力——惯性离心力。本文

每一块需要8个与非门快闪存储器芯片规划和 与非门芯片

每一块须要8个与非门快闪存储器芯片。计划和实行的全方位服务,由于,产量很少哈瓦那雪茄?和所谓的慢爱好宁静的单独散步龟,干部对策略,因而高利润产品补助亏损产品。而且他们本人的钱买良多人一拥而上,但许多人有什么用。增加陷阱把持扩大

声明:《R语言与非参数统计核密度估计 r语言密度估计曲线》为网友柠檬姓酸分享!如侵犯到您的合法权益请联系我们删除