Нахождение локальных максимумов и минимумов - PullRequest
59 голосов
/ 27 июля 2011

Я ищу вычислительно эффективный способ найти локальные максимумы / минимумы для большого списка чисел в R. Надеюсь, без for петли ...

Например, если у меня есть файл данных типа 1 2 3 2 1 1 2 1, я хочу, чтобы функция возвращала 3 и 7, которые являются позициями локальных максимумов.

Ответы [ 12 ]

52 голосов
/ 27 июля 2011

diff(diff(x)) (или diff(x,differences=2): благодаря @ZheyuanLi) по существу вычисляет дискретный аналог второй производной, поэтому должен быть отрицательным при локальных максимумах.+1 ниже учитывает тот факт, что результат diff короче входного вектора.

edit : добавлена ​​коррекция @ Tommy для случаев, когда delta-x не равен1 ...

tt <- c(1,2,3,2,1, 1, 2, 1)
which(diff(sign(diff(tt)))==-2)+1

Мое предложение выше (http://statweb.stanford.edu/~tibs/PPC/Rdist/) предназначено для случая, когда данные являются более шумными.

36 голосов
/ 27 июля 2011

@ Решение Бена довольно сладкое. Это не обрабатывает следующие случаи, хотя:

# all these return numeric(0):
x <- c(1,2,9,9,2,1,1,5,5,1) # duplicated points at maxima 
which(diff(sign(diff(x)))==-2)+1 
x <- c(2,2,9,9,2,1,1,5,5,1) # duplicated points at start
which(diff(sign(diff(x)))==-2)+1 
x <- c(3,2,9,9,2,1,1,5,5,1) # start is maxima
which(diff(sign(diff(x)))==-2)+1

Вот более надежная (и более медленная, уродливая) версия:

localMaxima <- function(x) {
  # Use -Inf instead if x is numeric (non-integer)
  y <- diff(c(-.Machine$integer.max, x)) > 0L
  rle(y)$lengths
  y <- cumsum(rle(y)$lengths)
  y <- y[seq.int(1L, length(y), 2L)]
  if (x[[1]] == x[[2]]) {
    y <- y[-1]
  }
  y
}

x <- c(1,2,9,9,2,1,1,5,5,1)
localMaxima(x) # 3, 8
x <- c(2,2,9,9,2,1,1,5,5,1)
localMaxima(x) # 3, 8
x <- c(3,2,9,9,2,1,1,5,5,1)
localMaxima(x) # 1, 3, 8
20 голосов
/ 27 июля 2011

Используйте функцию zoo library rollapply:

x <- c(1, 2, 3, 2, 1, 1, 2, 1)
library(zoo)
 xz <- as.zoo(x)
 rollapply(xz, 3, function(x) which.min(x)==2)
#    2     3     4     5     6     7 
#FALSE FALSE FALSE  TRUE FALSE FALSE 
 rollapply(xz, 3, function(x) which.max(x)==2)
#    2     3     4     5     6     7 
#FALSE  TRUE FALSE FALSE FALSE  TRUE 

Затем извлеките индекс, используя «coredata» для тех значений, где «which.max» является «центральным значением», сигнализирующим о локальном максимуме.Очевидно, вы могли бы сделать то же самое для локальных минимумов, используя which.min вместо which.max.

 rxz <- rollapply(xz, 3, function(x) which.max(x)==2)
 index(rxz)[coredata(rxz)]
#[1] 3 7

Я предполагаю, что вам не нужны начальные или конечные значения, но если вы это сделаете, вы можете дополнитьконцы ваших векторов перед обработкой, скорее как теломеры на хромосомах.

(Я отмечаю пакет ppc («Пиковые контрасты вероятности» для проведения масс-спектрометрического анализа, просто потому, что я не знал о его доступности до чтения@ Комментарий BenBolker выше, и я думаю, что добавление этих нескольких слов увеличит шансы, что кто-то с интересом массовой спецификации увидит это в поиске.)

12 голосов
/ 21 июня 2013

Есть несколько хороших решений, но это зависит от того, что вам нужно.

Просто diff(tt) возвращает различия.

Вы хотите определить, когда переходите от увеличения значений к уменьшению значений. Один из способов сделать это - @Ben:

 diff(sign(diff(tt)))==-2

Проблема здесь в том, что при этом будут обнаруживаться только те изменения, которые немедленно переходят от строго возрастающих к строго убывающим.

Небольшое изменение позволит повторить значения на пике (возвращая TRUE для последнего вхождения пика):

 diff(diff(x)>=0)<0

Тогда вам просто нужно правильно заполнить переднюю и заднюю части, если вы хотите обнаружить максимумы в начале или конце

Вот все, что заключено в функции (включая поиск долин):

 which.peaks <- function(x,partial=TRUE,decreasing=FALSE){
     if (decreasing){
         if (partial){
             which(diff(c(FALSE,diff(x)>0,TRUE))>0)
         }else {
             which(diff(diff(x)>0)>0)+1
         }
     }else {
         if (partial){
             which(diff(c(TRUE,diff(x)>=0,FALSE))<0)
         }else {
             which(diff(diff(x)>=0)<0)+1
         }
     }
 }
11 голосов
/ 28 марта 2017

Я сделал удар сегодня. Я знаю, что вы сказали, надеюсь, без циклов for, но я придерживался с помощью функции apply. Несколько компактный и быстрый и позволяет задавать пороговое значение, чтобы вы могли превысить 1.

Функция:

inflect <- function(x, threshold = 1){
  up   <- sapply(1:threshold, function(n) c(x[-(seq(n))], rep(NA, n)))
  down <-  sapply(-1:-threshold, function(n) c(rep(NA,abs(n)), x[-seq(length(x), length(x) - abs(n) + 1)]))
  a    <- cbind(x,up,down)
  list(minima = which(apply(a, 1, min) == a[,1]), maxima = which(apply(a, 1, max) == a[,1]))
}

Для визуализации / воспроизведения с порогами вы можете запустить следующий код:

# Pick a desired threshold # to plot up to
n <- 2
# Generate Data
randomwalk <- 100 + cumsum(rnorm(50, 0.2, 1)) # climbs upwards most of the time
bottoms <- lapply(1:n, function(x) inflect(randomwalk, threshold = x)$minima)
tops <- lapply(1:n, function(x) inflect(randomwalk, threshold = x)$maxima)
# Color functions
cf.1 <- grDevices::colorRampPalette(c("pink","red"))
cf.2 <- grDevices::colorRampPalette(c("cyan","blue"))
plot(randomwalk, type = 'l', main = "Minima & Maxima\nVariable Thresholds")
for(i in 1:n){
  points(bottoms[[i]], randomwalk[bottoms[[i]]], pch = 16, col = cf.1(n)[i], cex = i/1.5)
}
for(i in 1:n){
  points(tops[[i]], randomwalk[tops[[i]]], pch = 16, col = cf.2(n)[i], cex = i/1.5)
}
legend("topleft", legend = c("Minima",1:n,"Maxima",1:n), 
       pch = rep(c(NA, rep(16,n)), 2), col = c(1, cf.1(n),1, cf.2(n)), 
       pt.cex =  c(rep(c(1, c(1:n) / 1.5), 2)), cex = .75, ncol = 2)

enter image description here

3 голосов
/ 25 июля 2017

Ответ @ 42- это здорово, но у меня был случай, когда я не хотел использовать zoo.Это легко реализовать с помощью dplyr, используя lag и lead:

library(dplyr)
test = data_frame(x = sample(1:10, 20, replace = TRUE))
mutate(test, local.minima = if_else(lag(x) > x & lead(x) > x, TRUE, FALSE)

Как и в случае решения rollapply, вы можете управлять размером окна и краями с помощью lag / lead аргументы n и default соответственно.

2 голосов
/ 18 марта 2016

У меня были некоторые проблемы с получением местоположений для работы в предыдущих решениях, и я нашел способ получить минимумы и максимумы напрямую. Приведенный ниже код сделает это и отобразит его, отметив минимумы зеленым цветом и максимумы красным. В отличие от функции which.max() это вытянет все индексы минимумов / максимумов из фрейма данных. Нулевое значение добавляется в первой функции diff() для учета отсутствующей уменьшенной длины результата, возникающего при каждом использовании этой функции. Вставка этого в самый внутренний вызов функции diff() избавляет от необходимости добавлять смещение вне логического выражения. Это не имеет большого значения, но я чувствую, что это более чистый способ сделать это.

# create example data called stockData
stockData = data.frame(x = 1:30, y=rnorm(30,7))

# get the location of the minima/maxima. note the added zero offsets  
# the location to get the correct indices
min_indexes = which(diff(  sign(diff( c(0,stockData$y)))) == 2)
max_indexes = which(diff(  sign(diff( c(0,stockData$y)))) == -2)

# get the actual values where the minima/maxima are located
min_locs = stockData[min_indexes,]
max_locs = stockData[max_indexes,]

# plot the data and mark minima with red and maxima with green
plot(stockData$y, type="l")
points( min_locs, col="red", pch=19, cex=1  )
points( max_locs, col="green", pch=19, cex=1  )
2 голосов
/ 25 сентября 2014

Вот решение для минимумов :

@ решение Бена

x <- c(1,2,3,2,1,2,1)
which(diff(sign(diff(x)))==+2)+1 # 5

Пожалуйста, рассмотрите случаи на посту Томми!

@ Решение Томми:

localMinima <- function(x) {
  # Use -Inf instead if x is numeric (non-integer)
  y <- diff(c(.Machine$integer.max, x)) > 0L
  rle(y)$lengths
  y <- cumsum(rle(y)$lengths)
  y <- y[seq.int(1L, length(y), 2L)]
  if (x[[1]] == x[[2]]) {
    y <- y[-1]
  }
  y
}

x <- c(1,2,9,9,2,1,1,5,5,1)
localMinima(x) # 1, 7, 10
x <- c(2,2,9,9,2,1,1,5,5,1)
localMinima(x) # 7, 10
x <- c(3,2,9,9,2,1,1,5,5,1)
localMinima(x) # 2, 7, 10

Обратите внимание: ни localMaxima, ни localMinima не могут обрабатывать дублированные максимумы / минимумы при запуске!

0 голосов
/ 26 марта 2019

Нахождение локальных максимумов и минимумов для не такой простой последовательности, например 1 0 1 1 2 0 1 1 0 1 1 1 0 1 Я бы дал их позиции в (1), 5, 7.5, 11 и (14) для максимумов и 2, 6, 9, 13 для минимумов.

#Position                1 1 1 1 1
#      1 2 3 4 5 6 7 8 9 0 1 2 3 4
x <- c(1,0,1,1,2,0,1,1,0,1,1,1,0,1) #Frequency
#      p v     p v  p  v   p   v p  p..Peak, v..Valey

peakPosition <- function(x, inclBorders=TRUE) {
  if(inclBorders) {y <- c(min(x), x, min(x))
  } else {y <- c(x[1], x)}
  y <- data.frame(x=sign(diff(y)), i=1:(length(y)-1))
  y <- y[y$x!=0,]
  idx <- diff(y$x)<0
  (y$i[c(idx,F)] + y$i[c(F,idx)] - 1)/2
}

#Find Peaks
peakPosition(x)
#1.0  5.0  7.5 11.0 14.0

#Find Valeys
peakPosition(-x)
#2  6  9 13

peakPosition(c(1,2,3,2,1,1,2,1)) #3 7
0 голосов
/ 25 марта 2019

Поздно на вечеринку, но это может быть интересно для других. В настоящее время вы можете использовать (внутреннюю) функцию find_peaks из пакета ggpmisc. Вы можете параметризовать его, используя threshold, span и strict аргументы. Поскольку пакет ggpmisc предназначен для использования с ggplot2, вы можете напрямую построить минимумы и максимумы , используя функции stat_peaks и stat_valleys:

set.seed(1)
x <- 1:10
y <- runif(10)
# Maxima
x[ggpmisc:::find_peaks(y)]
[1] 4 7
y[ggpmisc:::find_peaks(y)]
[1] 0.9082078 0.9446753
# Minima
x[ggpmisc:::find_peaks(-y)]
[1] 5
y[ggpmisc:::find_peaks(-y)]
[1] 0.2016819    
# Plot
ggplot(data = data.frame(x, y), aes(x = x, y = y)) + geom_line() + stat_peaks(col = "red") + stat_valleys(col = "green")

enter image description here

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...