Znajdowanie lokalnych maksimów / minimów za pomocą Numpy w tablicy numpy 1D

Czy możesz zasugerować funkcję modułu z numpy / scipy, która może znaleźć lokalne maksima / minima w 1D tablicy numpy? Oczywiście najprostszym podejściem jest przyjrzenie się najbliższym sąsiadom, ale chciałbym mieć akceptowane rozwiązanie, które jest częścią numpy distro.

Jeśli szukasz wszystkich wpisów w tablicy 1d amniejszych niż ich sąsiedzi, możesz spróbować

numpy.r_[True, a[1:] < a[:-1]] & numpy.r_[a[:-1] < a[1:], True]

Możesz również wygładzić swoją tablicę przed tym krokiem, używając numpy.convolve().

Myślę, że nie ma do tego dedykowanej funkcji.

W SciPy> = 0,11

import numpy as np
from scipy.signal import argrelextrema

x = np.random.random(12)

# for local maxima
argrelextrema(x, np.greater)

# for local minima
argrelextrema(x, np.less)

Produkuje

>>> x
array([ 0.56660112,  0.76309473,  0.69597908,  0.38260156,  0.24346445,
    0.56021785,  0.24109326,  0.41884061,  0.35461957,  0.54398472,
    0.59572658,  0.92377974])
>>> argrelextrema(x, np.greater)
(array([1, 5, 7]),)
>>> argrelextrema(x, np.less)
(array([4, 6, 8]),)

Zauważ, że są to indeksy x, które są lokalnymi maks./min. Aby uzyskać wartości, spróbuj:

>>> x[argrelextrema(x, np.greater)[0]]

scipy.signalzapewnia również argrelmaxi argrelmindo znajdowania odpowiednio maksimów i minimów.

W przypadku krzywych z niezbyt dużym szumem polecam następujący mały fragment kodu:

from numpy import *

# example data with some peaks:
x = linspace(0,4,1e3)
data = .2*sin(10*x)+ exp(-abs(2-x)**2)

# that's the line, you need:
a = diff(sign(diff(data))).nonzero()[0] + 1 # local min+max
b = (diff(sign(diff(data))) > 0).nonzero()[0] + 1 # local min
c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max


# graphical output...
from pylab import *
plot(x,data)
plot(x[b], data[b], "o", label="min")
plot(x[c], data[c], "o", label="max")
legend()
show()

+1Jest ważny, ponieważ diffzmniejsza oryginalny numer indeksu.

Inne podejście (więcej słów, mniej kodu), które może pomóc:

Lokalizacje lokalnych maksimów i minimów są również lokalizacjami przejść przez zero pierwszej pochodnej. Na ogół znacznie łatwiej jest znaleźć przejście przez zero niż bezpośrednio znaleźć lokalne maksima i minima.

Niestety, pierwsza pochodna ma tendencję do „wzmacniania” szumu, więc gdy w oryginalnych danych obecny jest znaczny szum, pierwszą pochodną najlepiej jest zastosować dopiero po zastosowaniu pewnego stopnia wygładzenia oryginalnych danych.

Ponieważ wygładzanie jest, w najprostszym sensie, filtrem dolnoprzepustowym, wygładzanie jest często najlepsze (cóż, najłatwiej) przy użyciu jądra splotu, a "kształtowanie" tego jądra może zapewnić zaskakującą ilość możliwości zachowania / ulepszenia funkcji . Proces znajdowania optymalnego jądra można zautomatyzować za pomocą różnych środków, ale najlepszym może być zwykła brutalna siła (dużo szybka do znalezienia małych jąder). Dobre jądro będzie (zgodnie z przeznaczeniem) znacznie zniekształcić oryginalne dane, ale NIE wpłynie to na lokalizację interesujących nas szczytów / dolin.

Na szczęście dość często odpowiednie jądro można utworzyć za pomocą prostego pliku SWAG („wykształcone przypuszczenie”). Szerokość jądra wygładzającego powinna być nieco szersza niż najszerszy oczekiwany „interesujący” szczyt w oryginalnych danych, a jego kształt będzie przypominał ten pik (falka o pojedynczej skali). Dla jąder zachowujących średnią (czym powinien być każdy dobry filtr wygładzający), suma elementów jądra powinna być dokładnie równa 1,00, a jądro powinno być symetryczne względem swojego środka (co oznacza, że ​​będzie miało nieparzystą liczbę elementów.

Biorąc pod uwagę optymalne jądro wygładzające (lub niewielką liczbę jąder zoptymalizowanych pod kątem różnej zawartości danych), stopień wygładzenia staje się współczynnikiem skalującym („zysk”) jądra splotu.

Określenie „prawidłowego” (optymalnego) stopnia wygładzenia (wzmocnienia jądra splotu) można nawet zautomatyzować: Porównaj odchylenie standardowe danych pierwszej pochodnej z odchyleniem standardowym wygładzonych danych. Sposób, w jaki stosunek dwóch odchyleń standardowych zmienia się wraz ze zmianami stopnia wygładzania krzywki, aby przewidzieć efektywne wartości wygładzania. Wystarczy kilka ręcznych uruchomień danych (które są naprawdę reprezentatywne).

Wszystkie wcześniejsze rozwiązania opublikowane powyżej obliczają pierwszą pochodną, ​​ale nie traktują jej jako miary statystycznej, ani też powyższe rozwiązania nie próbują wykonywać wygładzania zachowującego / wzmacniającego funkcję (aby pomóc subtelnym skokom „przeskoczyć” nad szumem).

Na koniec zła wiadomość: znajdowanie „prawdziwych” szczytów staje się prawdziwym bólem, gdy szum ma również cechy, które wyglądają jak prawdziwe szczyty (nakładające się pasmo). Kolejnym bardziej złożonym rozwiązaniem jest generalnie użycie dłuższego jądra splotu („szersza apertura jądra”), które bierze pod uwagę zależność między sąsiednimi „rzeczywistymi” pikami (takimi jak minimalne lub maksymalne współczynniki występowania pików) lub użycie wielu splot przebiega przy użyciu jąder o różnych szerokościach (ale tylko wtedy, gdy jest szybszy: podstawową prawdą matematyczną jest to, że splot liniowy wykonywany w sekwencji zawsze może być razem splatany w pojedynczy splot). Często jednak o wiele łatwiej jest najpierw znaleźć sekwencję przydatnych jąder (o różnych szerokościach) i połączyć je razem, niż bezpośrednio znaleźć ostateczne jądro w jednym kroku.

Mamy nadzieję, że dostarczy to wystarczających informacji, aby umożliwić Google (i być może dobry tekst dotyczący statystyk) wypełnienie luk. Naprawdę chciałbym mieć czas, aby podać działający przykład lub link do niego. Jeśli ktoś natknie się na taki w Internecie, opublikuj go tutaj!

Począwszy od wersji 1.1 SciPy, możesz również użyć find_peaks . Poniżej znajdują się dwa przykłady zaczerpnięte z samej dokumentacji.

Używając tego heightargumentu, można wybrać wszystkie maksima powyżej określonego progu (w tym przykładzie wszystkie nieujemne maksima; może to być bardzo przydatne, jeśli mamy do czynienia z zaszumioną linią bazową; jeśli chcesz znaleźć minima, po prostu pomnóż wprowadzone dane autor -1:):

import matplotlib.pyplot as plt
from scipy.misc import electrocardiogram
from scipy.signal import find_peaks
import numpy as np

x = electrocardiogram()[2000:4000]
peaks, _ = find_peaks(x, height=0)
plt.plot(x)
plt.plot(peaks, x[peaks], "x")
plt.plot(np.zeros_like(x), "--", color="gray")
plt.show()

Kolejnym niezwykle pomocnym argumentem jest distanceokreślenie minimalnej odległości między dwoma szczytami:

peaks, _ = find_peaks(x, distance=150)
# difference between peaks is >= 150
print(np.diff(peaks))
# prints [186 180 177 171 177 169 167 164 158 162 172]

plt.plot(x)
plt.plot(peaks, x[peaks], "x")
plt.show()

Dlaczego nie skorzystać z wbudowanej funkcji Scipy signal.find_peaks_cwt do wykonania tego zadania?

from scipy import signal
import numpy as np

#generate junk data (numpy 1D arr)
xs = np.arange(0, np.pi, 0.05)
data = np.sin(xs)

# maxima : use builtin function to find (max) peaks
max_peakind = signal.find_peaks_cwt(data, np.arange(1,10))

# inverse  (in order to find minima)
inv_data = 1/data
# minima : use builtin function fo find (min) peaks (use inversed data)
min_peakind = signal.find_peaks_cwt(inv_data, np.arange(1,10))

#show results
print "maxima",  data[max_peakind]
print "minima",  data[min_peakind]

wyniki:

maxima [ 0.9995736]
minima [ 0.09146464]

pozdrowienia

Aktualizacja: nie byłem zadowolony z gradientu, więc uznałem, że jest bardziej niezawodny w użyciu numpy.diff. Daj mi znać, jeśli zrobisz to, co chcesz.

Jeśli chodzi o szum, matematycznym problemem jest zlokalizowanie maksimów / minimów, jeśli chcemy spojrzeć na szum, możemy użyć czegoś takiego jak splot, o którym wspomniano wcześniej.

import numpy as np
from matplotlib import pyplot

a=np.array([10.3,2,0.9,4,5,6,7,34,2,5,25,3,-26,-20,-29],dtype=np.float)

gradients=np.diff(a)
print gradients


maxima_num=0
minima_num=0
max_locations=[]
min_locations=[]
count=0
for i in gradients[:-1]:
        count+=1

    if ((cmp(i,0)>0) & (cmp(gradients[count],0)<0) & (i != gradients[count])):
        maxima_num+=1
        max_locations.append(count)     

    if ((cmp(i,0)<0) & (cmp(gradients[count],0)>0) & (i != gradients[count])):
        minima_num+=1
        min_locations.append(count)


turning_points = {'maxima_number':maxima_num,'minima_number':minima_num,'maxima_locations':max_locations,'minima_locations':min_locations}  

print turning_points

pyplot.plot(a)
pyplot.show()

Chociaż to pytanie jest naprawdę stare. Uważam, że w numpy jest znacznie prostsze podejście (jeden liniowiec).

import numpy as np

list = [1,3,9,5,2,5,6,9,7]

np.diff(np.sign(np.diff(list))) #the one liner

#output
array([ 0, -2,  0,  2,  0,  0, -2])

Aby znaleźć lokalne maksimum lub minimum, zasadniczo chcemy znaleźć, kiedy różnica między wartościami na liście (3-1, 9-3 ...) zmienia się z dodatniej na ujemną (maks.) Lub z ujemnej na dodatnią (min). Dlatego najpierw znajdujemy różnicę. Następnie znajdujemy znak, a następnie znajdujemy zmiany w znaku, ponownie biorąc różnicę. (Coś w rodzaju pierwszej i drugiej pochodnej w rachunku różniczkowym, tylko że mamy dane dyskretne i nie mamy funkcji ciągłej.)

Dane wyjściowe w moim przykładzie nie zawierają ekstrema (pierwszej i ostatniej wartości na liście). Podobnie jak w przypadku rachunku różniczkowego, jeśli druga pochodna jest ujemna, masz max, a jeśli jest dodatnia, masz min.

Tak więc mamy następujące zestawienie:

[1,  3,  9,  5,  2,  5,  6,  9,  7]
    [0, -2,  0,  2,  0,  0, -2]
        Max     Min         Max

Żadne z tych rozwiązań nie zadziałało, ponieważ chciałem również znaleźć piki w centrum powtarzających się wartości. na przykład w

ar = np.array([0,1,2,2,2,1,3,3,3,2,5,0])

odpowiedź powinna brzmieć

array([ 3,  7, 10], dtype=int64)

Zrobiłem to za pomocą pętli. Wiem, że nie jest super czysty, ale spełnia swoje zadanie.

def findLocalMaxima(ar):
# find local maxima of array, including centers of repeating elements    
maxInd = np.zeros_like(ar)
peakVar = -np.inf
i = -1
while i < len(ar)-1:
#for i in range(len(ar)):
    i += 1
    if peakVar < ar[i]:
        peakVar = ar[i]
        for j in range(i,len(ar)):
            if peakVar < ar[j]:
                break
            elif peakVar == ar[j]:
                continue
            elif peakVar > ar[j]:
                peakInd = i + np.floor(abs(i-j)/2)
                maxInd[peakInd.astype(int)] = 1
                i = j
                break
    peakVar = ar[i]
maxInd = np.where(maxInd)[0]
return maxInd 

import numpy as np
x=np.array([6,3,5,2,1,4,9,7,8])
y=np.array([2,1,3,5,3,9,8,10,7])
sortId=np.argsort(x)
x=x[sortId]
y=y[sortId]
minm = np.array([])
maxm = np.array([])
i = 0
while i < length-1:
    if i < length - 1:
        while i < length-1 and y[i+1] >= y[i]:
            i+=1

        if i != 0 and i < length-1:
            maxm = np.append(maxm,i)

        i+=1

    if i < length - 1:
        while i < length-1 and y[i+1] <= y[i]:
            i+=1

        if i < length-1:
            minm = np.append(minm,i)
        i+=1


print minm
print maxm

minmi maxmzawierają odpowiednio wskaźniki minimów i maksimów. W przypadku ogromnego zestawu danych da wiele maksym / minimów, więc w takim przypadku najpierw wygładź krzywą, a następnie zastosuj ten algorytm.

Inne rozwiązanie wykorzystujące zasadniczo operatora dylatacji:

import numpy as np
from scipy.ndimage import rank_filter

def find_local_maxima(x):
   x_dilate = rank_filter(x, -1, size=3)
   return x_dilate == x

a dla minimów:

def find_local_minima(x):
   x_erode = rank_filter(x, -0, size=3)
   return x_erode == x

Również z scipy.ndimagewas może zastąpić rank_filter(x, -1, size=3)z grey_dilationi rank_filter(x, 0, size=3)z grey_erosion. Nie będzie to wymagało sortowania lokalnego, więc jest nieco szybsze.

Inny:

def local_maxima_mask(vec):
    """
    Get a mask of all points in vec which are local maxima
    :param vec: A real-valued vector
    :return: A boolean mask of the same size where True elements correspond to maxima. 
    """
    mask = np.zeros(vec.shape, dtype=np.bool)
    greater_than_the_last = np.diff(vec)>0  # N-1
    mask[1:] = greater_than_the_last
    mask[:-1] &= ~greater_than_the_last
    return mask