使用 Numpy 在 1D numpy 数组中查找局部最大值/最小值

问题描述：

您能否推荐一个 numpy/scipy 中的模块函数，用于在 1D numpy 数组中找到局部最大值/最小值？显然，最简单的方法是查看最近的邻居，但我希望有一个公认的解决方案，它是 numpy 发行版的一部分。

解决方案 1：

在 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 的局部最大值/最小值索引。要获取这些值，请尝试：

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

scipy.signal还提供argrelmax和argrelmin用于分别寻找最大值和最小值。

解决方案 2：

如果你正在寻找 1d 数组中所有a小于其邻居的条目，你可以尝试

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

您还可以在此步骤之前使用 来平滑numpy.convolve()您的阵列。

我认为没有专门的功能可以实现这一点。

解决方案 3：

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

解决方案 4：

对于噪声不太大的曲线，我推荐以下小代码片段：

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

很+1重要，因为它diff减少了原始索引号。

解决方案 5：

另一种方法（更多的文字，更少的代码）可能会有所帮助：

局部极大值和极小值的位置也是一阶导数的零交叉点的位置。通常，找到零交叉点比直接找到局部极大值和极小值要容易得多。

不幸的是，一阶导数往往会“放大”噪声，因此当原始数据中存在明显噪声时，最好仅在对原始数据应用一定程度的平滑后才使用一阶导数。

由于平滑在最简单的意义上是一种低通滤波器，因此平滑通常最好（最容易）通过使用卷积核来完成，并且“塑造”该内核可以提供令人惊讶的特征保留/增强能力。可以使用各种方法自动找到最佳内核，但最好的方法可能是简单的蛮力（对于找到小内核来说速度足够快）。好的内核将（如预期的那样）大量扭曲原始数据，但不会影响感兴趣的峰/谷的位置。

幸运的是，通常可以通过简单的 SWAG（“有根据的猜测”）创建合适的核。平滑核的宽度应比原始数据中预期最宽的“有趣”峰值稍宽，其形状将类似于该峰值（单尺度小波）。对于均值保留核（任何良好的平滑滤波器都应该如此），核元素的总和应精确等于 1.00，并且核应关于其中心对称（这意味着它将具有奇数个元素。

给定一个最佳平滑核（或针对不同数据内容优化的少量核），平滑程度成为卷积核（“增益”）的比例因子。

确定“正确”（最佳）平滑度（卷积核增益）甚至可以自动化：将一阶导数数据的标准差与平滑数据的标准差进行比较。两个标准差的比率如何随平滑度的变化而变化可用于预测有效平滑值。只需进行几次手动数据运行（真正具有代表性）就足够了。

上面发布的所有先前的解决方案都计算了一阶导数，但它们并不将其视为统计测量，上述解决方案也没有尝试执行特征保留/增强平滑（以帮助细微峰值“超越”噪声）。

最后，坏消息是：当噪声也具有看起来像真实峰值（重叠带宽）的特征时，找到“真实”峰值会变得非常麻烦。下一个更复杂的解决方案通常是使用更长的卷积核（“更宽的核孔径”），其中考虑相邻“真实”峰值之间的关系（例如峰值发生的最小或最大速率），或者使用具有不同宽度的核进行多次卷积（但前提是它更快：这是一个基本的数学事实，即按顺序执行的线性卷积总是可以卷积成单个卷积）。但首先找到一系列有用的核（宽度不同）并将它们卷积在一起通常比直接在单个步骤中找到最终的核要容易得多。

希望这能提供足够的信息，让 Google（也许还有一篇好的统计文本）填补空白。我真的希望我有时间提供一个可行的示例或一个链接。如果有人在网上看到一个，请在这里发布！

解决方案 6：

我相信 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...）中值之间的差值何时从正变为负（最大值）或从负变为正（最小值）。因此，我们首先找到差值。然后我们找到符号，然后我们通过再次取差值来找到符号的变化。（有点像微积分中的一阶和二阶导数，只是我们有离散数据，没有连续函数。）

我的示例中的输出不包含极值（列表中的第一个和最后一个值）。此外，就像微积分一样，如果二阶导数为负，则有最大值，如果为正，则有最小值。

因此，我们有以下对决：

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

解决方案 7：

为什么不使用 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]

问候

解决方案 8：

更新：
我对渐变不太满意，所以我发现使用它更可靠numpy.diff。

关于噪音问题，数学问题是找到最大值/最小值，如果我们想看噪音，我们可以使用前面提到的卷积之类的方法。

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

解决方案 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分别包含最小值和最大值的指标。对于庞大的数据集，它将给出大量的最大值/最小值，因此在这种情况下，首先平滑曲线，然后应用此算法。

解决方案 11：

另一种解决方案本质上是使用扩张运算符：

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

对于最小值：

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

另外，scipy.ndimage你可以rank_filter(x, -1, size=3)用grey_dilation和rank_filter(x, 0, size=3)替换grey_erosion。这不需要本地排序，因此速度会稍微快一些。

解决方案 12：

另一个：

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

解决方案 13：

还有...另一个答案。

这个不需要额外的包（numpy 除外）。例如，

points = [ 0, 0, 1, 2, 3, 3, 2, 2, 3, 1, 1 ]
minimums   ^  ^              ^  ^     ^  ^

将返回所有局部最小值的列表

result = [ 0, 1, 6, 7, 9, 10 ]

它可以很容易地扩展以寻找最大值。

def find_valleys(points: np.ndarray, edges=True) -> list:
    """
    Find the indices of all points that are local minimums.

    :param np.ndarray points: a 1D array of numeric data
    :param bool edges: allows the first and last indices to be returned, defaults to True
    :return list: a list of integers, indices into the array
    """
    dif = np.diff(points)
    p = -1 if edges else 1
    s = 0
    result = []
    for i,d in enumerate(dif):
        if d < 0: s = i + 1
        if p < 0 and d > 0:   # found a valley
            result.extend(range(s,i + 1))
        if d: p = d
    if p < 0 and edges:
        result.extend(range(s,i + 2))
    return result