# 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? • You can find the n'th largest number in R by: sort(x, TRUE)[n] - where x is vector of numbers. Aug 10, 2017 at 6:38 • @Repmat but the second largest number in density values might not be the value of the second peak. It could be a point next to the first peak which is still larger than the second peak. – John Aug 10, 2017 at 6:43 • You're looking for local maxima. The answer to this question could help. Aug 10, 2017 at 6:43 • Alternatively, you could write a function with some optimization technique (think Newton's method etc) to identify the maxima in density(D)$y and then store them away. Then sort the vector, and voila! Aug 10, 2017 at 6:46

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


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.

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