Como encontrar picos / vales locais em uma série de dados?

Aqui está o meu experimento:

Estou usando a findPeaksfunção no pacote quantmod :

Quero detectar picos "locais" dentro de uma tolerância 5, ou seja, os primeiros locais após a série cronológica caírem dos picos locais em 5:

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

A saída é

[1] 3 22 41

Parece errado, pois espero mais "picos locais" do que 3 ...

Alguma ideia?

A fonte desse código é obtida digitando seu nome no prompt R. A saída é

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
}

O teste x[pks - 1] - x[pks] > threshcompara cada valor de pico com o valor imediatamente seguinte na série (não com o próximo nível da série). Ele usa uma estimativa (bruta) do tamanho da inclinação da função imediatamente após o pico e seleciona apenas os picos em que a inclinação excede threshde tamanho. No seu caso, apenas os três primeiros picos são suficientemente nítidos para passar no teste. Você detectará todos os picos usando o padrão:

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

Concordo com a resposta do whuber, mas só queria acrescentar que a parte "+2" do código, que tenta alterar o índice para corresponder ao pico recém-encontrado, na verdade 'ultrapassa' e deve ser "+1". por exemplo, no exemplo em questão, obtemos:

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

quando destacamos esses picos encontrados em um gráfico (vermelho negrito):

vemos que eles estão consistentemente a 1 ponto do pico real.

consequentemente

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

deve ser pks[x[pks] - x[pks + 1] > thresh]oupks[x[pks] - x[pks - 1] > thresh]

ATUALIZAÇÃO GRANDE

seguindo minha própria busca para encontrar uma função adequada de localização de pico, escrevi isso:

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
}

um 'pico' é definido como um máximo local, com mpontos em ambos os lados menores que ele. portanto, quanto maior o parâmetro m, mais rigoroso é o procedimento de pico de financiamento. assim:

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

a função também pode ser usada para encontrar mínimos locais de qualquer vetor seqüencial xvia find_peaks(-x).

Nota: agora coloquei a função no gitHub se alguém precisar: https://github.com/stas-g/findPeaks

Eek: atualização menor. Eu tive que mudar duas linhas de código, os limites (adicionar -1 e +1) para alcançar a equivalência com a função de Stas_G (estava encontrando muitos "picos extras" em conjuntos de dados reais). Peço desculpas a qualquer pessoa que se desvie muito pouco do meu post original.

Eu tenho usado o algoritmo de encontrar picos do Stas_g há algum tempo. Foi benéfico para mim em um dos meus projetos posteriores devido à sua simplicidade. No entanto, eu precisava usá-lo milhões de vezes para um cálculo, então eu o reescrevi no Rcpp (consulte o pacote Rcpp). É aproximadamente 6x mais rápido que a versão R em testes simples. Se alguém estiver interessado, adicionei o código abaixo. Espero que ajude alguém, Saúde!

Algumas pequenas advertências. Esta função retorna os índices de pico na ordem inversa do código R. Requer uma função interna de sinal C ++, que eu incluí. Não foi completamente otimizado, mas não são esperados mais ganhos de desempenho.

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

Primeiramente: o algoritmo também chama falsamente uma queda à direita de um platô plano, pois sign(diff(x, na.pad = FALSE)) será 0 e -1, de modo que seu diff também será -1. Uma correção simples é garantir que a diferença de sinal que antecede a entrada negativa não seja zero, mas 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: o algoritmo fornece resultados muito locais, por exemplo, um 'up' seguido de um 'down' em qualquer execução de três termos consecutivos na sequência. Se alguém está interessado no máximo local de uma função contínua barulhenta, então - provavelmente existem outras coisas melhores por aí, mas esta é a minha solução barata e imediata

1. identifique os picos primeiro usando a média de execução de 3 pontos consecutivos para
   suavizar os dados levemente. Empregue também o controle acima mencionado contra queda plana e depois queda.

2. filtre esses candidatos comparando, para uma versão suavizada, a média dentro de uma janela centralizada em cada pico com a média dos termos locais externos.

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

É verdade que a função também identifica o fim dos platôs, mas acho que há outra solução mais fácil: como o primeiro diff de um pico real resultará em '1' e em '-1', o segundo diff seria '-2', e podemos verificar diretamente

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

usando Numpy

ser = np.random.randint(-40, 40, 100) # 100 points
peak = np.where(np.diff(ser) < 0)[0]

ou

double_difference = np.diff(np.sign(np.diff(ser)))
peak = np.where(double_difference == -2)[0]

usando pandas

ser = pd.Series(np.random.randint(2, 5, 100))
peak_df = ser[(ser.shift(1) < ser) & (ser.shift(-1) < ser)]
peak = peak_df.index