R语言 核密度估计的峰值

mwg9r5ms  于 2023-09-27  发布在  其他
关注(0)|答案(4)|浏览(114)

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

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

但是当计算d$y时,精确的函数是已知的。我怎样才能找到众数的精确值?

1yjd4xko

1yjd4xko1#

这里有两个处理模式的函数。dmode函数找到具有最高峰值的模式(主导模式),并且nmode标识模式的数量。

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))
rsl1atfo

rsl1atfo2#

如果我理解了你的问题,我想你只是想对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
5q4ezhmt

5q4ezhmt3#

我认为你需要两个步骤来归档你需要的东西:
1)找到KDE峰值的x轴值
2)得到峰值的desnity值
所以(如果你不介意使用一个包)使用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']]

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

soat7uwm

soat7uwm4#

使用'multimode'包的一行程序:

multimode::locmodes(x)$locations[1]
Estimated location
Mode: 1.238839 

Estimated value of the density
Mode: 0.1605735 

Critical bandwidth: 2.155968

要在KDE的图上显示计算模式的位置,请添加display = TRUE参数。

multimode::locmodes(x, display = TRUE)

可以使用mod0参数返回多个模式(及其反模式)的位置;模取locmodes(x)$locations的奇数位置,反模取偶数位置。

相关问题