I need to find as precisely as possible the peak of the kernel density estimation (modal value of the continuous random variable). I can find the approximate value:
x<-rlnorm(100)
d<-density(x)
plot(d)
i<-which.max(d$y)
d$y[i]
d$x[i]
But when calculating d$y
precise function is known. How can I locate the exact value of the mode?
Here are two functions for dealing with modes. The dmode function finds the mode with the highest peak (dominate mode) and n.modes identify the number of 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))