R How to find the secondary peak of a distribution For example, if I generate density with R code
D = c(rnorm(100,1,1), rnorm(100,5,1))

Then the following density will follow:

We can find the primary peak location by 
density(D)$x[which(density(D)$y == max(density(D)$y))]

But how to find the secondary peak?
 A: You can use mixture models to capture the biomodality
library(flexmix)
set.seed(42)

D <- c(rnorm(100,1,1), rnorm(100,5,1))
kde <- density(D)
m1 <- FLXMRglm(family = "gaussian")
m2 <- FLXMRglm(family = "gaussian")
fit <- flexmix(D ~ 1, data = as.data.frame(D), k = 2, model = list(m1, m2))
c1 <- parameters(fit, component=1)[[1]]
c2 <- parameters(fit, component=2)[[1]]




> c1
                  Comp.1
coef.(Intercept) 1.022880
sigma            1.031319


> c2
                  Comp.2
coef.(Intercept) 4.9042434
sigma            0.9081448

plot(kde)
abline(v=1, col='blue')
abline(v=c1[[1]], lty=2, col='blue')
abline(v=5, col='red')
abline(v=c2[[1]], lty=2, col='red')


A: If you can assume that you have a mixture of normal distributions, simply use a mixture model:
set.seed(42)
D = c(rnorm(100,1,1), rnorm(100,5,1))

library(mixtools)
mD <- normalmixEM(D)
mD$mu
#[1] 1.079553 4.918794
summary(mD)

plot(mD, which=2)
lines(density(D, "SJ"), lwd = 2)


If you really need the exact peak locations of the combined density function, you have all necessary values available (mixing proportion, means and standard deviations) for calculating the maxima. I don't have time to figure out the maths right now, but it shouldn't be too hard.
A: I realize that I am late to the game on this one, but I solved this problem using the sm and features packages...
library(sm)
library(features)

D = c(rnorm(100,1,1), rnorm(100,5,1))

#find points on density curve
kde <- sm::sm.density(D)
Dcurve <- data.frame(x = kde$eval.points, y = kde$estimate)

#find critical points (slope = 0) of curve
f <- features::features(x = Dcurve$x, y = Dcurve$y)

#filter for points with negative curvature (maxima)
f$cpts[f$curvature < 0]

[1] 1.06 4.83

