Найти локальный минимум функции python

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

Можете ли вы предложить функцию модуля из numpy/scipy, которая может найти локальные максимумы/минимумы в массиве 1D numpy? Очевидно, что самый простой подход – это взглянуть на ближайших соседей, но я хотел бы иметь принятое решение, которое является частью дистрибутива numpy.

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

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

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

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 для нахождения максимумов и минимумов соответственно.

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

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 уменьшает исходный номер индекса.

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

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

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

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

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

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

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

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

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

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

Почему бы не использовать встроенную функцию 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] 

Начиная с версии 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() 

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 = print turning_points pyplot.plot(a) pyplot.show() 

Хотя этот вопрос действительно старый. Я считаю, что в 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 

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

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

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]  

 # numpy всегда пригодится import numpy as np # Для построения графиков import matplotlib.pyplot as plt # Научный питон from scipy import optimize # Символьный питон from sympy import * 

a, b, h, k, x, v1 = symbols('a, b, h, k, x, v1') 

f = (sqrt(x**2+a**2) + k*sqrt(b**2+(h-x)**2))/v1 

>> df (k*(-h + x)/sqrt(b**2 + (h - x)**2) + x/sqrt(a**2 + x**2))/v1 

params = a:10, b:10, h:10.0, k:5, v1: 5> 

>> df_par x/(5*sqrt(x**2 + 100)) + (x - 10.0)/sqrt((-x + 10.0)**2 + 100) 

df_num = lambda xnum: df_par.subs( x: xnum> ) 

sol = optimize.root(df_num, 5.0)