2013-04-27 1 views
8

Мне нужно как можно точнее определить пик оценки плотности ядра (модальное значение непрерывной случайной величины). Я могу найти приближенное значение:Пик оценки плотности ядра

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

Но при расчете d$y точная функция известна. Как определить точное значение режима?

ответ

5

Если я понимаю ваш вопрос, я думаю, что вы просто хотите более тонкую дискретизацию x и y. Для этого вы можете изменить значение n в функции density (по умолчанию 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 
+0

Спасибо! ;) он работает довольно точно – 16per9

10

Вот две функции для работы с режимами. Функция 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)) 
0

Я думаю, вам нужно два шага в архив, что вам нужно:

1) Найти значение х-оси пика KDE

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']]) 

Edit: Вы можете также использовать

hdrResult[['falpha']] 

, если он достаточно точен для вас!