Recherche de maxima et minima locaux

Question

Je cherche un moyen efficace de trouver les maxima/minima locaux pour une grande liste de nombres dans R. J'espère qu'il n'y aura pas de for boucles...

Par exemple, si j'ai un fichier de données comme 1 2 3 2 1 1 2 1 Je veux que la fonction renvoie 3 et 7, qui sont les positions des maxima locaux.

Answer 1

diff(diff(x)) (ou diff(x,differences=2) : merci à @ZheyuanLi) calcule essentiellement l'analogue discret de la dérivée seconde, et devrait donc être négatif aux maxima locaux. Les +1 ci-dessous prend en compte le fait que le résultat de l'opération diff est plus court que le vecteur d'entrée.

éditer : ajout de la correction de @Tommy pour les cas où delta-x n'est pas 1...

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

Ma suggestion ci-dessus ( http://statweb.stanford.edu/~tibs/PPC/Rdist/ ) est destiné au cas où les données sont plus bruyantes.

Answer 2

La solution de @Ben est très intéressante. Elle ne gère cependant pas les cas suivants :

# 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

Voici une version plus robuste (et plus lente, plus laide) :

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

Answer 3

Utilisez la fonction rollapply de la bibliothèque zoo :

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 

Il faut ensuite extraire l'index à l'aide des "coredata" pour les valeurs pour lesquelles "which.max" est une "valeur centrale" signalant un maximum local. Vous pouvez évidemment faire la même chose pour les minima locaux en utilisant which.min au lieu de which.max .

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

Je suppose que vous ne voulez pas les valeurs de départ ou d'arrivée, mais si c'est le cas, vous pourriez remplir les extrémités de vos vecteurs avant le traitement, un peu comme le font les télomères sur les chromosomes.

(Je mentionne le paquet ppc ("Peak Probability Contrasts" pour faire des analyses de spectrométrie de masse, simplement parce que je n'étais pas au courant de sa disponibilité avant de lire le commentaire de @BenBolker ci-dessus, et je pense que l'ajout de ces quelques mots augmentera les chances que quelqu'un qui s'intéresse à la spectrométrie de masse voie ceci lors d'une recherche).

Answer 4

J'ai fait un essai aujourd'hui. Je sais que vous avez dit qu'il fallait espérer sans boucle for, mais je me suis contenté d'utiliser la fonction apply. Elle est assez compacte et rapide et permet de spécifier un seuil, de sorte que l'on peut aller au-delà de 1.

La fonction :

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

Pour le visualiser/jouer avec les seuils, vous pouvez exécuter le code suivant :

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

Answer 5

Il existe de bonnes solutions, mais tout dépend de ce dont vous avez besoin.

Juste diff(tt) renvoie les différences.

Vous souhaitez détecter le moment où vous passez de valeurs croissantes à des valeurs décroissantes. Une façon de le faire est fournie par @Ben :

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

Le problème est que cette méthode ne permet de détecter que les changements qui passent immédiatement d'une augmentation stricte à une diminution stricte.

Une légère modification permettra d'obtenir des valeurs répétées au niveau du pic (retour à l'état initial). TRUE pour la dernière occurrence de la valeur maximale) :

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

Il suffit ensuite de remplir correctement le recto et le verso si l'on veut détecter des maxima au début ou à la fin de la période de référence.

Voici tout ce qui se trouve dans une fonction (y compris la recherche des vallées) :

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