Encontrar máximos y mínimos locales

Estoy buscando una forma computacionalmente eficiente de encontrar máximos / mínimos locales para una gran lista de números en R. Con suerte sin bucles for

Por ejemplo, si tengo un archivo de datos como 1 2 3 2 1 1 2 1 , quiero que la función devuelva 3 y 7, que son las posiciones de los máximos locales.

diff(diff(x)) (o diff(x,differences=2) : gracias a @ZheyuanLi) básicamente calcula el análogo discreto de la segunda derivada, por lo que debe ser negativo a máximos locales. El +1 continuación se ocupa del hecho de que el resultado de diff es más corto que el vector de entrada.

editar : se agregó @ la corrección de Tommy para los casos donde delta-x no es 1 …

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

Mi sugerencia anterior ( http://statweb.stanford.edu/~tibs/PPC/Rdist/ ) está destinada al caso en que los datos son más ruidosos.

La solución de @ Ben es muy dulce. Sin embargo, no maneja los siguientes casos:

 # 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 

Aquí hay una versión más robusta (y más lenta, más fea):

 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 

Use la función de biblioteca del zoológico 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 

A continuación, extraiga el índice utilizando los 'coredatos' para aquellos valores en los que 'which.max' es un "valor central" que indica un máximo local. Obviamente podrías hacer lo mismo para los mínimos locales usando which.min lugar de which.max .

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

Supongo que no quieres los valores iniciales o finales, pero si lo haces, puedes rellenar los extremos de tus vectores antes de procesarlos, como hacen los telómeros en los cromosomas.

(Estoy observando el paquete ppc ("Contrastes de probabilidad máxima" para realizar análisis de espectrometría de masas, simplemente porque no estaba al tanto de su disponibilidad hasta leer el comentario de BenBolker anterior, y creo que agregar estas pocas palabras boostá las posibilidades de que alguien con un el interés de la especificación masiva verá esto en una búsqueda).

Apuñalé esto hoy. Sé que lo dijo con suerte sin bucles pero me quedé con el uso de la función Aplicar. Algo compacto y rápido, y permite la especificación de umbral para que pueda ir más allá de 1.

La función:

 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])) } 

Para visualizarlo / jugar con umbrales puede ejecutar el siguiente código:

 # 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

Se proporcionan algunas buenas soluciones, pero depende de lo que necesite.

Solo diff(tt) devuelve las diferencias.

Desea detectar cuándo pasa de valores crecientes a valores decrecientes. Una forma de hacerlo es provista por @Ben:

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

El problema aquí es que esto solo detectará los cambios que van de estrictamente boost a estrictamente decrecientes.

Un ligero cambio permitirá valores repetidos en el pico (regresando TRUE para la última ocurrencia del valor máximo):

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

Entonces, simplemente necesita rellenar correctamente la parte frontal y posterior si desea detectar los máximos al principio o al final de

Aquí está todo envuelto en una función (incluido el hallazgo de valles):

  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 } } } 

Aquí está la solución para los mínimos :

La solución de @ Ben

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

Por favor, considere los casos en la publicación de Tommy!

@ La solución de Tommy:

 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 

Tenga en cuenta que ni localMaxima ni localMinima pueden manejar máximos / mínimos duplicados al inicio.

Tuve algunos problemas para lograr que las ubicaciones funcionaran en soluciones anteriores y se me ocurrió una forma de obtener los mínimos y los máximos directamente. El siguiente código hará esto y lo trazará, marcando los mínimos en verde y los máximos en rojo. A diferencia de la función which.max() esto which.max() todos los índices de los mínimos / máximos de un dataframe. El valor cero se agrega en la primera función diff() para dar cuenta de la menor longitud perdida del resultado que ocurre cada vez que usa la función. Insertar esto en la llamada de función diff() más interna ahorra tener que agregar un desplazamiento fuera de la expresión lógica. No importa mucho, pero creo que es una forma más limpia de hacerlo.

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

Answer by @ 42- es genial, pero tuve un caso de uso en el que no quería usar el zoo . Es fácil implementar esto con dplyr usando lag y 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) 

Al igual que la solución rollapply , puede controlar el tamaño de la ventana y los casos de borde a través de los argumentos lag / lead n y default , respectivamente.

En el paquete pracma , usa el

 tt <- c(1,2,3,2,1, 1, 2, 1) tt_peaks <- findpeaks(tt, zero = "0", peakpat = NULL, minpeakheight = -Inf, minpeakdistance = 1, threshold = 0, npeaks = 0, sortstr = FALSE) [,1] [,2] [,3] [,4] [1,] 3 3 1 5 [2,] 2 7 6 8 

Eso devuelve una matriz con 4 columnas. La primera columna muestra los valores absolutos de los picos locales. La segunda columna son los índices La tercera y la cuarta columna son el inicio y el final de los picos (con superposición potencial).

Consulte https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/findpeaks para obtener más detalles.

Una advertencia: lo usé en una serie de números enteros, y el pico fue un índice demasiado tarde (para todos los picos) y no sé por qué. Así que tuve que eliminar manualmente "1" de mi vector de índice (no es gran cosa).

Lo publiqué en otro lugar, pero creo que esta es una forma interesante de hacerlo. No estoy seguro de cuál es su eficiencia computacional, pero es una forma muy concisa de resolver el problema.

 vals=rbinom(1000,20,0.5) text=paste0(substr(format(diff(vals),scientific=TRUE),1,1),collapse="") sort(na.omit(c(gregexpr('[ ]-',text)[[1]]+1,ifelse(grepl('^-',text),1,NA), ifelse(grepl('[^-]$',text),length(vals),NA))))