Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
364 views
in Technique[技术] by (71.8m points)

r - Extracting peak values from a density ggplot

Here is a minimal example of my problem.

I would like a vertical line passing through the two peaks in the graph and a label of its corresponding x axis.

    library(dplyr)
    library(ggplot2)

   data.frame(Values=c(1,1,1,2,2,3,3,5,4,5,5,5,5,6,6,6,6,5,5,5,
                        4,4,4,7,8,8,8,8,8,8,8,4,3,2,2,3,5,6,7,4,
                        2,2,3,9,9,9,9,9,9,9,9,9,9,5,1,3,3,1))%>%
    ggplot(aes(x=Values))+
      geom_density()+
      theme_classic()
question from:https://stackoverflow.com/questions/65886842/extracting-peak-values-from-a-density-ggplot

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

Building on Ben Bolker's comment:

This code doesn't strictly tell you all the mathematical local maxima, it just tells you points which are the max within a width-20 window whose endpoints don't equal the max.

values <- c(1,1,1,2,2,3,3,5,4,5,5,5,5,6,6,6,6,5,5,5,
                        4,4,4,7,8,8,8,8,8,8,8,4,3,2,2,3,5,6,7,4,
                        2,2,3,9,9,9,9,9,9,9,9,9,9,5,1,3,3,1)

dens <- density(values)

i1 <- 
  zoo::rollapply(dens$y, 20, function(x){
    if(head(x, 1) < max(x) & tail(x, 1) < max(x))
      which.max(x)
    else 
      NA
  })

max_inds <- unique(i1 + seq_along(i1) - 1)
x_at_max <- dens$x[max_inds[!is.na(max_inds)]]


data.frame(Values=c(1,1,1,2,2,3,3,5,4,5,5,5,5,6,6,6,6,5,5,5,
                        4,4,4,7,8,8,8,8,8,8,8,4,3,2,2,3,5,6,7,4,
                        2,2,3,9,9,9,9,9,9,9,9,9,9,5,1,3,3,1))%>%
    ggplot(aes(x=Values))+
      geom_density()+
      theme_classic() +
    geom_vline(xintercept = x_at_max)

enter image description here

Edit:

Actually, looks like this gives the same result as the approach in the answer linked by Ben Bolker. You should probably use that one. I assume it's faster, and also it will give you all the local maximua even if they're not the max of a width-20 window.

r <- rle(dens$y)
which(rep(x = diff(sign(diff(c(-Inf, r$values, -Inf)))) == -2, 
          times = r$lengths))
# [1] 240 378

max_inds[!is.na(max_inds)]
# [1] 240 378

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...