Поиск локальных максимумов/минимумов с помощью Numpy в массиве 1D numpy

Ответ 1

Если вы ищете все записи в массиве 1d a меньше, чем их соседи, вы можете попробовать

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

Вы также можете smooth ваш массив до этого шага с помощью numpy.convolve().

Я не думаю, что для этого есть специальная функция.

Ответ 2

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

Производит

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

Обратите внимание, что это индексы x, которые являются локальными max/min. Чтобы получить значения, попробуйте:

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

scipy.signal также обеспечивает argrelmax и argrelmin для нахождения максимумов и минимумов соответственно.

Ответ 3

Для кривых с не слишком большим шумом я рекомендую следующий небольшой фрагмент кода:

from numpy import *

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

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

+1 важен, поскольку diff уменьшает исходный номер индекса.

Ответ 4

Другой подход (больше слов, меньше кода), который может помочь:

Расположение локальных максимумов и минимумов также является местом пересечения нуля первой производной. Как правило, гораздо легче найти пересечения нуля, чем непосредственно находить локальные максимумы и минимумы.

К сожалению, первая производная имеет тенденцию "усиливать" шум, поэтому, когда в исходных данных присутствует значительный шум, первая производная лучше всего использовать только после того, как исходные данные имеют некоторую степень сглаживания.

Поскольку сглаживание в простейшем смысле является фильтром нижних частот, сглаживание часто лучше всего (ну, что наиболее легко) выполняется с использованием ядра свертки, и "формирование" этого ядра может обеспечить удивительное количество сохраняющих функцию/повышение возможностей. Процесс нахождения оптимального ядра может быть автоматизирован с использованием различных средств, но лучшим может быть простая грубая сила (достаточно много для поиска небольших ядер). Хорошее ядро ​​будет (как и предполагалось) массивно искажать исходные данные, но не будет влиять на местоположение интересующих вас пиков/долин.

К счастью, довольно часто подходящее ядро ​​может быть создано с помощью простого SWAG ( "образованное предположение" ). Ширина сглаживающего ядра должна быть немного шире, чем самый ожидаемый "интересный" пик в исходных данных, и его форма будет напоминать этот пик (одномасштабный вейвлет). Для ядер, сохраняющих среднее значение (каким должен быть любой хороший сглаживающий фильтр) сумма элементов ядра должна быть точно равна 1.00, а ядро ​​должно быть симметричным относительно своего центра (что означает, что у него будет нечетное число элементов.

Учитывая оптимальное сглаживающее ядро ​​(или небольшое количество ядер, оптимизированных для различного содержимого данных), степень сглаживания становится фактором масштабирования для ( "усиления" ) ядра свертки.

Определение "правильной" (оптимальной) степени сглаживания (усиление ядра свертки) может быть даже автоматизировано: сравните стандартное отклонение данных первой производной со стандартным отклонением сглаженных данных. Как отношение двух стандартных отклонений изменяется с изменением степени сглаживания кулачка, чтобы предсказать эффективные значения сглаживания. Необходимо выполнить несколько ручных операций с данными (которые действительно являются репрезентативными).

Все предыдущие решения, вышеперечисленные выше, вычисляют первую производную, но они не рассматривают ее как статистическую меру и не выполняют вышеуказанные решения для выполнения функции сохранения/улучшения сглаживания (чтобы помочь тонким пикам "прыгать выше" шума).

Наконец, плохая новость: поиск "реальных" пиков становится королевской болью, когда шум также имеет функции, которые выглядят как реальные пики (перекрывающиеся полосы пропускания). Следующее более комплексное решение, как правило, должно использовать более длинное сверточное ядро ​​( "более широкую апертуру ядра" ), которое учитывает взаимосвязь между смежными "реальными" пиками (такими как минимальные или максимальные скорости для возникновения пика) или использовать множественные свертки передаются с использованием ядер разной ширины (но только в том случае, если они быстрее: фундаментальная математическая истина состоит в том, что линейные свертки, выполняемые в последовательности, всегда могут быть свернуты вместе в одну свертку). Но часто гораздо легче сначала найти последовательность полезных ядер (различной ширины) и свернуть их вместе, чем напрямую находить последнее ядро ​​за один шаг.

Надеемся, это даст достаточно информации, чтобы позволить Google (и, возможно, хороший текст статистики) заполнить пробелы. Мне очень жаль, что у меня не было времени предоставить обработанный пример или ссылку на него. Если кто-то наткнулся на один онлайн, отправьте его здесь!

Ответ 5

Почему бы не использовать встроенную функцию Scipy signal.find_peaks_cwt для выполнения этой работы?

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]

Результаты:

maxima [ 0.9995736]
minima [ 0.09146464]

С уважением

Ответ 6

Начиная с версии SciPy 1.1, вы также можете использовать find_peaks. Ниже приведены два примера, взятых из самой документации.

Используя аргумент height, можно выбрать все максимумы выше определенного порога (в этом примере все неотрицательные максимумы; это может быть очень полезно, если приходится иметь дело с шумной базовой линией; если вы хотите найти минимумы, просто умножьте ввод -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()

Еще один чрезвычайно полезный аргумент - это distance, которое определяет минимальное расстояние между двумя пиками:

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

Ответ 7

Update: Я был недоволен градиентом, поэтому я нашел более надежным использование numpy.diff. Пожалуйста, дайте мне знать, если он делает то, что вы хотите.

Что касается проблемы шума, математическая проблема заключается в том, чтобы найти максимумы/минимумы, если мы хотим посмотреть на шум, мы можем использовать нечто вроде convolve, о котором говорилось ранее.

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

Ответ 8

Хотя этот вопрос действительно старый. Я считаю, что в numpy гораздо более простой подход (один лайнер).

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

Чтобы найти локальный макс или мин, мы, по сути, хотим найти, когда разница между значениями в списке (3-1, 9-3...) изменяется от положительного к отрицательному (макс) или отрицательному к положительному (мин). Поэтому сначала мы находим разницу. Затем мы находим знак, а затем находим изменения знака, снова принимая разность. (Вроде как первая и вторая производные в исчислении, только у нас есть дискретные данные и не имеют непрерывной функции.)

Результат в моем примере не содержит экстремумов (первое и последнее значения в списке). Кроме того, как и исчисление, если вторая производная отрицательна, у вас есть max, и если она положительная, у вас есть мин.

Таким образом, мы имеем следующий матч:

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

Ответ 9

Ни одно из этих решений не работало для меня, так как я хотел найти пики в центре повторяющихся значений. например, в

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

ответ должен быть

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

Я сделал это, используя цикл. Я знаю, что он не супер чистый, но он выполняет свою работу.

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 

Ответ 10

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

minm и maxm содержат индексы минимумов и максимумов соответственно. Для огромного набора данных он даст много максим/минимумов, поэтому в этом случае сначала сгладьте кривую, а затем примените этот алгоритм.