2013-04-27 24 views
8

Muszę znaleźć jak najdokładniej pik oceny gęstości jądra (wartość modalna ciągłej zmiennej losowej). mogę znaleźć wartość przybliżoną:Pik oszacowania gęstości ziarnistości

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

ale kiedy obliczania d$y precyzyjna funkcja jest znana. Jak mogę zlokalizować dokładną wartość trybu?

Odpowiedz

5

Jeśli dobrze rozumiem twoje pytanie, myślę, że po prostu pragniesz lepszej dyskretyzacji x i y. Aby to zrobić, możesz zmienić wartość n w funkcji density (domyślnie jest to n=512).

Na przykład porównać

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

z:

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

Dzięki! ;) wydaje się działać dość dokładnie – 16per9

10

Oto dwie funkcje do czynienia z trybów. Funkcja dmode znajduje tryb z najwyższym pikiem (tryb dominacji), a n.mode określa liczbę trybów.

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

myślę, że potrzebne są dwa kroki w celu archiwizacji, czego potrzebujesz:

1) Znajdź wartość oś x piku KDE

2) Get wartość desnity piku

Tak (jeśli nie przeszkadza przy użyciu pakietu) rozwiązania przy użyciu pakietu hdrcde będzie wyglądać następująco:

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: Można również użyć

hdrResult[['falpha']] 

jeśli jest wystarczająco dokładna dla Ciebie!