Wie finde ich lokale Gipfel / Täler in einer Reihe von Daten?

Hier ist mein Experiment:

Ich benutze die findPeaksFunktion im quantmod- Paket:

Ich möchte "lokale" Peaks innerhalb einer Toleranz 5 erkennen, dh die ersten Stellen nach der Zeitreihe fallen um 5 von den lokalen Peaks ab:

aa=100:1
bb=sin(aa/3)
cc=aa*bb
plot(cc, type="l")
p=findPeaks(cc, 5)
points(p, cc[p])
p

Die Ausgabe ist

[1] 3 22 41

Es scheint falsch, da ich mehr "lokale Spitzen" als 3 erwarte ...

Irgendwelche Gedanken?

Die Quelle dieses Codes wird durch Eingabe seines Namens an der R-Eingabeaufforderung ermittelt. Die Ausgabe ist

function (x, thresh = 0) 
{
    pks <- which(diff(sign(diff(x, na.pad = FALSE)), na.pad = FALSE) < 0) + 2
    if (!missing(thresh)) {
        pks[x[pks - 1] - x[pks] > thresh]
    }
    else pks
}

Der Test x[pks - 1] - x[pks] > threshvergleicht jeden Spitzenwert mit dem Wert, der ihm unmittelbar in der Reihe folgt (nicht mit dem nächsten Tiefpunkt in der Reihe). Es verwendet eine (grobe) Schätzung der Größe der Steigung der Funktion unmittelbar nach dem Peak und wählt nur die Peaks aus, bei denen diese Steigung die threshGröße überschreitet . In Ihrem Fall sind nur die ersten drei Peaks scharf genug, um den Test zu bestehen. Sie erkennen alle Peaks anhand der Standardeinstellung:

> findPeaks(cc)
[1]  3 22 41 59 78 96

Ich stimme der Antwort von whuber zu, wollte aber nur hinzufügen, dass der "+2" -Teil des Codes, der versucht, den Index zu verschieben, um dem neu gefundenen Peak zu entsprechen, tatsächlich "überschießt" und "+1" sein sollte. Im vorliegenden Beispiel erhalten wir zum Beispiel:

> findPeaks(cc)
[1]  3 22 41 59 78 96

Wenn wir diese gefundenen Peaks in einem Diagramm hervorheben (fett rot):

wir sehen, dass sie konstant 1 Punkt vom tatsächlichen Peak entfernt sind.

Konsequenz

pks[x[pks - 1] - x[pks] > thresh]

sollte pks[x[pks] - x[pks + 1] > thresh]oder seinpks[x[pks] - x[pks - 1] > thresh]

GROSSES UPDATE

Nach meiner eigenen Suche nach einer adäquaten Peak-Finding-Funktion habe ich folgendes geschrieben:

find_peaks <- function (x, m = 3){
    shape <- diff(sign(diff(x, na.pad = FALSE)))
    pks <- sapply(which(shape < 0), FUN = function(i){
       z <- i - m + 1
       z <- ifelse(z > 0, z, 1)
       w <- i + m + 1
       w <- ifelse(w < length(x), w, length(x))
       if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
    })
     pks <- unlist(pks)
     pks
}

Ein „Peak“ ist definiert als ein lokales Maximum, bei dem die mPunkte auf beiden Seiten kleiner sind als es. mJe größer der Parameter , desto strenger ist daher das Spitzenfinanzierungsverfahren. so:

find_peaks(cc, m = 1)
[1]  2 21 40 58 77 95

Die Funktion kann auch verwendet werden, um lokale Minima eines sequentiellen Vektors xüber zu finden find_peaks(-x).

Hinweis: Ich habe jetzt die Funktion auf gitHub gesetzt, falls jemand sie benötigt: https://github.com/stas-g/findPeaks

Eek: Kleines Update. Ich musste zwei Codezeilen, die Grenzen, ändern (ein -1 und ein +1 hinzufügen), um die Gleichwertigkeit mit der Funktion von Stas_G zu erreichen (es wurden einige zu viele zusätzliche Peaks in realen Datensätzen gefunden). Entschuldigungen für irgendjemanden, der von meinem ursprünglichen Beitrag nur geringfügig abgekommen ist.

Ich verwende den Stas_g-Algorithmus zum Auffinden von Peaks seit einiger Zeit. Für eines meiner späteren Projekte war es mir wegen seiner Einfachheit von Vorteil. Ich musste es jedoch millionenfach für eine Berechnung verwenden, damit ich es in Rcpp umschrieb (siehe Rcpp-Paket). In einfachen Tests ist es ungefähr 6x schneller als die R-Version. Wenn jemand interessiert ist, habe ich den Code unten hinzugefügt. Hoffentlich helfe ich jemandem, Prost!

Einige kleine Vorbehalte. Diese Funktion gibt Spitzenindizes in umgekehrter Reihenfolge des R-Codes zurück. Es erfordert eine interne C ++ - Signierfunktion, die ich integriert habe. Es wurde nicht vollständig optimiert, aber weitere Leistungssteigerungen werden nicht erwartet.

//This function returns the sign of a given real valued double.
// [[Rcpp::export]]
double signDblCPP (double x){
  double ret = 0;
  if(x > 0){ret = 1;}
  if(x < 0){ret = -1;}
  return(ret);
}

//Tested to be 6x faster(37 us vs 207 us). This operation is done from 200x per layer
//Original R function by Stas_G
// [[Rcpp::export]]
NumericVector findPeaksCPP( NumericVector vY, int m = 3) {
  int sze = vY.size();
  int i = 0;//generic iterator
  int q = 0;//second generic iterator

  int lb = 0;//left bound
  int rb = 0;//right bound

  bool isGreatest = true;//flag to state whether current index is greatest known value

  NumericVector ret(1);
  int pksFound = 0;

  for(i = 0; i < (sze-2); ++i){
    //Find all regions with negative laplacian between neighbors
    //following expression is identical to diff(sign(diff(xV, na.pad = FALSE)))
    if(signDblCPP( vY(i + 2)  - vY( i + 1 ) ) - signDblCPP( vY( i + 1 )  - vY( i ) ) < 0){
      //Now assess all regions with negative laplacian between neighbors...
      lb = i - m - 1;// define left bound of vector
      if(lb < 0){lb = 0;}//ensure our neighbor comparison is bounded by vector length
      rb = i + m + 1;// define right bound of vector
      if(rb >= (sze-2)){rb = (sze-3);}//ensure our neighbor comparison is bounded by vector length
      //Scan through loop and ensure that the neighbors are smaller in magnitude
      for(q = lb; q < rb; ++q){
        if(vY(q) > vY(i+1)){ isGreatest = false; }
      }

      //We have found a peak by our criterion
      if(isGreatest){
        if(pksFound > 0){//Check vector size.
         ret.insert( 0, double(i + 2) );
       }else{
         ret(0) = double(i + 2);
        }
        pksFound = pksFound + 1;
      }else{ // we did not find a peak, reset location is peak max flag.
        isGreatest = true;
      }//End if found peak
    }//End if laplace condition
  }//End loop
  return(ret);
}//End Fn

Erstens: Der Algorithmus ruft fälschlicherweise einen Drop rechts von einem flachen Plateau auf, da sign(diff(x, na.pad = FALSE)) er 0 und dann -1 ist, sodass sein Diff ebenfalls -1 ist. Eine einfache Lösung besteht darin, sicherzustellen, dass das Vorzeichen vor dem negativen Eintrag nicht null, sondern positiv ist:

    n <- length(x)
    dx.1 <- sign(diff(x, na.pad = FALSE))
    pks <- which(diff(dx.1, na.pad = FALSE) < 0 & dx.1[-(n-1)] > 0) + 1

Zweitens: Der Algorithmus liefert sehr lokale Ergebnisse, z. B. ein "Auf", gefolgt von einem "Ab" in jedem Lauf von drei aufeinanderfolgenden Begriffen in der Sequenz. Wenn man sich stattdessen für lokale Maxima einer gestörten kontinuierlichen Funktion interessiert, dann - gibt es wahrscheinlich andere bessere Dinge da draußen, aber dies ist meine billige und sofortige Lösung

1. Identifizieren Sie die Peaks zuerst mit dem laufenden Durchschnitt von 3 aufeinanderfolgenden Punkten,
   um die Daten ein wenig zu glätten. Verwenden Sie auch die oben erwähnte Kontrolle gegen Flach- und Abfall.

2. Filtern Sie diese Kandidaten, indem Sie für eine Version mit Lössglättung den Durchschnitt innerhalb eines Fensters, das an jedem Peak zentriert ist, mit dem Durchschnitt der lokalen Begriffe außerhalb des Fensters vergleichen.

   "myfindPeaks" <- 
   function (x, thresh=0.05, span=0.25, lspan=0.05, noisey=TRUE)
   {
     n <- length(x)
     y <- x
     mu.y.loc <- y
     if(noisey)
     {
       mu.y.loc <- (x[1:(n-2)] + x[2:(n-1)] + x[3:n])/3
       mu.y.loc <- c(mu.y.loc[1], mu.y.loc, mu.y.loc[n-2])
     }
     y.loess <- loess(x~I(1:n), span=span)
     y <- y.loess[[2]]
     sig.y <- var(y.loess$resid, na.rm=TRUE)^0.5
     DX.1 <- sign(diff(mu.y.loc, na.pad = FALSE))
     pks <- which(diff(DX.1, na.pad = FALSE) < 0 & DX.1[-(n-1)] > 0) + 1
     out <- pks
     if(noisey)
     {
       n.w <- floor(lspan*n/2)
       out <- NULL
       for(pk in pks)
       {
         inner <- (pk-n.w):(pk+n.w)
         outer <- c((pk-2*n.w):(pk-n.w),(pk+2*n.w):(pk+n.w))
         mu.y.outer <- mean(y[outer])
         if(!is.na(mu.y.outer)) 
           if (mean(y[inner])-mu.y.outer > thresh*sig.y) out <- c(out, pk)
       }
     }
     out
   }

Es ist wahr, dass die Funktion auch das Ende von Plateaus identifiziert, aber ich denke, dass es eine andere einfachere Lösung gibt: Da der erste Diff eines realen Peaks '1' ergibt, dann '-1', wäre der zweite Diff '-2', und wir können direkt überprüfen

    pks <- which(diff(sign(diff(x, na.pad = FALSE)), na.pad = FALSE) < 1) + 1