核密度估计的峰值

Posted

技术标签:

【中文标题】核密度估计的峰值【英文标题】:Peak of the kernel density estimation 【发布时间】:2013-04-21 18:14:40 【问题描述】:

我需要尽可能精确地找到(连续随机变量的模态值)。我可以找到近似值:

x<-rlnorm(100)
d<-density(x)
plot(d)
i<-which.max(d$y)
d$y[i]
d$x[i]

但是在计算d$y 时,精确的函数是已知的。如何找到模式的确切值?

【问题讨论】:

【参考方案1】:

这里有两个处理模式的函数。 dmode 函数查找具有最高峰值的模式(主导模式),n.modes 标识模式的数量。

    dmode <- function(x) 
      den <- density(x, kernel=c("gaussian"))
        ( den$x[den$y==max(den$y)] )   
      

    n.modes <- function(x)   
       den <- density(x, kernel=c("gaussian"))
       den.s <- smooth.spline(den$x, den$y, all.knots=TRUE, spar=0.8)
         s.0 <- predict(den.s, den.s$x, deriv=0)
         s.1 <- predict(den.s, den.s$x, deriv=1)
       s.derv <- data.frame(s0=s.0$y, s1=s.1$y)
       nmodes <- length(rle(den.sign <- sign(s.derv$s1))$values)/2
       if ((nmodes > 10) == TRUE)  nmodes <- 10 
          if (is.na(nmodes) == TRUE)  nmodes <- 0  
       ( nmodes )
    

# Example
x <- runif(1000,0,100)
  plot(density(x))
    abline(v=dmode(x))

【讨论】:

【参考方案2】:

如果我理解您的问题,我认为您只是想要对 xy 进行更精细的离散化。为此,您可以在density 函数中更改n 的值(默认为n=512)。

例如比较

set.seed(1)
x = rlnorm(100)
d = density(x)
i = which.max(d$y)
d$y[i]; d$x[i]
0.4526; 0.722

与:

d = density(x, n=1e6)
i = which.max(d$y)
d$y[i]; d$x[i]
0.4525; 0.7228

【讨论】:

【参考方案3】:

我认为您需要两个步骤来归档您需要的内容:

1) 求 KDE 峰值的 x 轴值

2) 获取峰的密度值

因此(如果您不介意使用软件包)使用 hdrcde 软件包的解决方案如下所示:

require(hdrcde)

x<-rlnorm(100)
d<-density(x)

# calcualte KDE with help of the hdrcde package
hdrResult<-hdr(den=d,prob=0)

# define the linear interpolation function for the density estimation
dd<-approxfun(d$x,d$y)
# get the density value of the KDE peak
vDens<-dd(hdrResult[['mode']])

编辑:您也可以使用

hdrResult[['falpha']]

如果它对你来说足够精确!

【讨论】:

以上是关于核密度估计的峰值的主要内容,如果未能解决你的问题,请参考以下文章

非参数核密度估计是用啥软件实现

如何估计密度函数并计算其峰值?

如何估计密度函数并计算其峰值?

什么是gis核密度计算

r语言绘制核密度图怎么计算重叠

比较核密度估计图