¿Cómo encontrar picos/valles locales en una serie de datos?

Question

Este es mi experimento:

Estoy utilizando el findPeaks en la función quantmod paquete:

Quiero detectar picos "locales" dentro de una tolerancia 5, es decir, las primeras ubicaciones después de que la serie temporal se aleje de los picos locales en 5:

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

El resultado es

[1] 3 22 41

Me parece mal, ya que espero más "picos locales" que 3...

¿Alguna idea?

Answer 1

Estoy de acuerdo con la respuesta de whuber, pero sólo quería añadir que el "+2" parte del código, que intenta cambiar el índice para que coincida con el pico recién encontrado en realidad 'overshoots' y debe ser "+1". por ejemplo, en el ejemplo que nos ocupa obtenemos:

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

cuando resaltamos estos picos encontrados en un gráfico (negrita roja):

vemos que están constantemente a 1 punto del pico real.

en consecuencia

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

debe ser pks[x[pks] - x[pks + 1] > thresh] o pks[x[pks] - x[pks - 1] > thresh]

GRAN ACTUALIZACIÓN

siguiendo mi propia búsqueda para encontrar una función adecuada para encontrar picos escribí esto:

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
}

un "pico" se define como un máximo local con m puntos a cada lado de él que sean menores que él. por lo tanto, cuanto mayor sea el parámetro m más estricto es el procedimiento de financiación de picos. así:

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

la función también puede utilizarse para encontrar mínimos locales de cualquier vector secuencial x vía find_peaks(-x) .

Nota: ahora he puesto la función en gitHub por si alguien la necesita: https://github.com/stas-g/findPeaks

Answer 2

La fuente de este código se obtiene escribiendo su nombre en el indicador R. La salida es

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
}

La prueba x[pks - 1] - x[pks] > thresh compara cada valor de pico con el valor inmediatamente después de la serie (no a la siguiente depresión de la serie). Utiliza una estimación (aproximada) del tamaño del pendiente de la función inmediatamente después del pico y selecciona sólo aquellos picos en los que dicha pendiente supera thresh en tamaño. En tu caso, sólo los tres primeros picos son suficientemente nítidos para pasar la prueba. Detectarás todos los picos utilizando el valor predeterminado:

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

Answer 3

Eek: Actualización menor. Tuve que cambiar dos líneas de código, los límites, (añadir un -1 y +1) para llegar a la equivalencia con la función de Stas_G (que estaba encontrando algunos demasiados 'picos extra' en conjuntos de datos reales). Disculpas a todos los que se hayan extraviado por mi mensaje original.

Llevo bastante tiempo utilizando el algoritmo find peaks de Stas_g. Me resultó beneficioso para uno de mis últimos proyectos debido a su simplicidad. Sin embargo, necesitaba usarlo millones de veces para un cálculo así que lo reescribí en Rcpp(Ver paquete Rcpp). Es aproximadamente 6 veces más rápido que la versión R en pruebas sencillas. Si alguien está interesado he añadido el código a continuación. Espero ayudar a alguien, ¡Salud!

Algunas advertencias menores. Esta función devuelve los índices de los picos en orden inverso al código R. Requiere una función interna de C++ Sign, que he incluido. No se ha optimizado completamente, pero no se espera que aumente el rendimiento.

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

Answer 4

En primer lugar: El algoritmo también llama falsamente una caída a la derecha de una meseta plana ya que sign(diff(x, na.pad = FALSE)) será 0 y luego -1 para que su diff también será -1. Una solución sencilla consiste en asegurarse de que el signo de la diferencia que precede a la entrada negativa no sea cero, sino positivo. entrada negativa no sea cero sino positiva:

    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

Segundo: El algoritmo da muy resultados locales, por ejemplo, un "arriba" seguido de un "abajo" en cualquier serie de tres términos consecutivos en la secuencia. Si uno está interesado en los máximos locales de una función continua con ruido, entonces -- probablemente hay otras cosas mejores por ahí, pero esta es mi solución barata e inmediata

1. identificar primero los picos utilizando la media de 3 puntos consecutivos para
   suavizar ligeramente los datos. Emplea también el control antes mencionado contra la caída de los datos planos.

2. filtrar estos candidatos comparando, para una versión suavizada de loess, la media dentro de una ventana centrada en cada pico con la media de los términos locales fuera de ella.

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

Answer 5

Es cierto que la función también identifica el final de las mesetas, pero creo que hay otra solución más fácil: Dado que la primera diferencia de un pico real resultará en '1' y luego en '-1', la segunda diferencia sería '-2', y podemos comprobarlo directamente. segundo diff sería '-2', y podemos comprobar directamente

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