Finden lokaler Maxima / Minima mit Numpy in einem 1D-Numpy-Array

Können Sie eine Modulfunktion von numpy / scipy vorschlagen, die lokale Maxima / Minima in einem 1D-numpy-Array findet? Natürlich ist der einfachste Ansatz, einen Blick auf die nächsten Nachbarn zu werfen, aber ich hätte gerne eine akzeptierte Lösung, die Teil der numpy Distribution ist.

Wenn Sie nach allen Einträgen im 1d-Array suchen, die akleiner als ihre Nachbarn sind, können Sie es versuchen

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

Sie können Ihr Array auch vor diesem Schritt mit glättennumpy.convolve() .

Ich glaube nicht, dass es dafür eine spezielle Funktion gibt.

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

Produziert

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

Beachten Sie, dass dies die Indizes von x sind, die lokal max / min sind. Um die Werte zu erhalten, versuchen Sie:

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

scipy.signalbietet auch argrelmaxund argrelminzum Finden von Maxima bzw. Minima.

Für Kurven mit nicht zu viel Rauschen empfehle ich das folgende kleine Code-Snippet:

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

Das +1ist wichtig, weil diffdie ursprüngliche Indexnummer reduziert wird.

Ein anderer Ansatz (mehr Wörter, weniger Code), der helfen kann:

Die Orte der lokalen Maxima und Minima sind auch die Orte der Nulldurchgänge der ersten Ableitung. Es ist im Allgemeinen viel einfacher, Nulldurchgänge zu finden, als lokale Maxima und Minima direkt zu finden.

Leider neigt die erste Ableitung dazu, das Rauschen zu "verstärken". Wenn also in den Originaldaten signifikantes Rauschen vorhanden ist, wird die erste Ableitung am besten erst verwendet, nachdem auf die Originaldaten ein gewisser Grad an Glättung angewendet wurde.

Da das Glätten im einfachsten Sinne ein Tiefpassfilter ist, wird das Glätten häufig am besten (am einfachsten) unter Verwendung eines Faltungskerns durchgeführt, und das "Formen" dieses Kernels kann eine überraschende Menge an Funktionen zum Erhalten / Verbessern von Merkmalen bereitstellen . Der Prozess des Findens eines optimalen Kernels kann mit einer Vielzahl von Mitteln automatisiert werden, aber das Beste kann einfache Brute Force sein (viel schnell, um kleine Kernel zu finden). Ein guter Kernel wird (wie beabsichtigt) die Originaldaten massiv verzerren, aber NICHT die Position der interessierenden Spitzen / Täler beeinflussen.

Glücklicherweise kann ziemlich oft ein geeigneter Kernel über eine einfache SWAG ("fundierte Vermutung") erstellt werden. Die Breite des Glättungskerns sollte etwas breiter sein als der breiteste erwartete "interessante" Peak in den Originaldaten, und seine Form ähnelt diesem Peak (einem einskalierten Wavelet). Für mittelschonende Kernel (was ein guter Glättungsfilter sein sollte) sollte die Summe der Kernelelemente genau gleich 1,00 sein, und der Kernel sollte symmetrisch zu seiner Mitte sein (was bedeutet, dass er eine ungerade Anzahl von Elementen hat.

Bei einem optimalen Glättungskern (oder einer kleinen Anzahl von Kerneln, die für unterschiedliche Dateninhalte optimiert sind) wird der Grad der Glättung zu einem Skalierungsfaktor für (den "Gewinn") des Faltungskerns.

Die Bestimmung des "richtigen" (optimalen) Glättungsgrades (Faltungskernverstärkung) kann sogar automatisiert werden: Vergleichen Sie die Standardabweichung der Daten der ersten Ableitung mit der Standardabweichung der geglätteten Daten. Wie sich das Verhältnis der beiden Standardabweichungen mit Änderungen des Glättungsgrads ändert, kann verwendet werden, um effektive Glättungswerte vorherzusagen. Ein paar manuelle Datenläufe (die wirklich repräsentativ sind) sollten alles sein, was benötigt wird.

Alle oben aufgeführten früheren Lösungen berechnen die erste Ableitung, behandeln sie jedoch weder als statistisches Maß, noch versuchen die oben genannten Lösungen, die Glättung von Merkmalen zu erhalten / zu verbessern (um subtilen Spitzen zu helfen, über das Rauschen zu "springen").

Schließlich die schlechte Nachricht: Das Finden von "echten" Peaks wird zu einem königlichen Schmerz, wenn das Rauschen auch Merkmale aufweist, die wie echte Peaks aussehen (überlappende Bandbreite). Die nächste komplexere Lösung besteht im Allgemeinen darin, einen längeren Faltungskern (eine "breitere Kernelöffnung") zu verwenden, der die Beziehung zwischen benachbarten "realen" Peaks (wie minimale oder maximale Raten für das Auftreten von Peaks) berücksichtigt, oder mehrere zu verwenden Faltungsdurchläufe werden mit Kerneln unterschiedlicher Breite durchgeführt (aber nur, wenn es schneller ist: Es ist eine grundlegende mathematische Wahrheit, dass nacheinander durchgeführte lineare Faltungen immer zusammen zu einer einzigen Faltung gefaltet werden können). Es ist jedoch oft viel einfacher, zuerst eine Folge nützlicher Kernel (unterschiedlicher Breite) zu finden und sie zusammenzufalten, als den endgültigen Kernel in einem einzigen Schritt direkt zu finden.

Hoffentlich bietet dies genügend Informationen, damit Google (und möglicherweise ein guter Statistiktext) die Lücken füllen kann. Ich wünschte wirklich, ich hätte die Zeit, ein funktionierendes Beispiel oder einen Link zu einem zu liefern. Wenn jemand online auf eines stößt, poste es bitte hier!

Ab SciPy Version 1.1 können Sie auch find_peaks verwenden . Nachfolgend finden Sie zwei Beispiele aus der Dokumentation.

Mit dem heightArgument können alle Maxima über einem bestimmten Schwellenwert ausgewählt werden (in diesem Beispiel alle nicht negativen Maxima; dies kann sehr nützlich sein, wenn man sich mit einer verrauschten Grundlinie befassen muss; wenn Sie Minima finden möchten, multiplizieren Sie einfach Ihre Eingabe von -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()

Ein weiteres äußerst hilfreiches Argument ist distance, das den Mindestabstand zwischen zwei Spitzen definiert:

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

Warum nicht die in Scipy integrierte Funktion signal.find_peaks_cwt verwenden , um die Arbeit zu erledigen?

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]

Ergebnisse:

maxima [ 0.9995736]
minima [ 0.09146464]

Grüße

Update: Ich war mit dem Farbverlauf nicht zufrieden und fand es daher zuverlässiger numpy.diff. Bitte lassen Sie mich wissen, ob es tut, was Sie wollen.

In Bezug auf das Problem des Rauschens besteht das mathematische Problem darin, Maxima / Minima zu lokalisieren, wenn wir das Rauschen betrachten wollen, können wir so etwas wie Faltung verwenden, das zuvor erwähnt wurde.

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

Während diese Frage wirklich alt ist. Ich glaube, es gibt einen viel einfacheren Ansatz bei Numpy (einem Einzeiler).

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

Um ein lokales Maximum oder Min zu finden, möchten wir im Wesentlichen herausfinden, wann sich die Differenz zwischen den Werten in der Liste (3-1, 9-3 ...) von positiv zu negativ (max) oder negativ zu positiv (min) ändert. Deshalb finden wir zuerst den Unterschied. Dann finden wir das Zeichen, und dann finden wir die Zeichenänderungen, indem wir den Unterschied erneut nehmen. (Ähnlich wie bei einer ersten und zweiten Ableitung im Kalkül haben nur wir diskrete Daten und keine stetige Funktion.)

Die Ausgabe in meinem Beispiel enthält nicht die Extrema (den ersten und den letzten Wert in der Liste). Ebenso wie beim Kalkül haben Sie, wenn die zweite Ableitung negativ ist, max und wenn sie positiv ist, haben Sie eine min.

Somit haben wir folgendes Matchup:

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

Keine dieser Lösungen funktionierte für mich, da ich auch Spitzenwerte im Zentrum sich wiederholender Werte finden wollte. zum Beispiel in

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

Die Antwort sollte sein

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

Ich habe das mit einer Schleife gemacht. Ich weiß, dass es nicht super sauber ist, aber es erledigt den Job.

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

minmund maxmenthalten Indizes von Minima bzw. Maxima. Bei einem großen Datensatz werden viele Maxima / Minimas angegeben. In diesem Fall wird die Kurve zuerst geglättet und dann dieser Algorithmus angewendet.

Eine andere Lösung, die im Wesentlichen einen erweiterten Operator verwendet:

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

und für die Minima:

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

Auch von können scipy.ndimageSie rank_filter(x, -1, size=3)mit grey_dilationund rank_filter(x, 0, size=3)mit ersetzen grey_erosion. Dies erfordert keine lokale Sortierung und ist daher etwas schneller.

Noch einer:

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