r - 核仁密度估算峰

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

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

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

最佳答案

这是处理模式的两个功能。 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))