4
$\begingroup$

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:

enter image description here

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?

$\endgroup$
4
  • $\begingroup$ You can find the n'th largest number in R by: sort(x, TRUE)[n] - where x is vector of numbers. $\endgroup$
    – Repmat
    Aug 10, 2017 at 6:38
  • 2
    $\begingroup$ @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. $\endgroup$
    – John
    Aug 10, 2017 at 6:43
  • 2
    $\begingroup$ You're looking for local maxima. The answer to this question could help. $\endgroup$ Aug 10, 2017 at 6:43
  • 2
    $\begingroup$ 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! $\endgroup$ Aug 10, 2017 at 6:46

3 Answers 3

7
$\begingroup$

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

enter image description here

$\endgroup$
3
$\begingroup$

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)

resulting plot

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.

$\endgroup$
0
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.