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

Question

你能推荐一种从numpy/scipy模块中找到一维numpy数组中局部最大值/最小值的功能吗？显然，最简单的方法是查看最近邻居，但我希望有一个被numpy批准的解决方案。

- Navi

---

可能是 如何高效地在NumPy中找到平滑多维数组的局部最小值？ 的重复问题。 - Sven Marnach

不，那是在2D空间（我说的是1D）。它涉及到自定义函数。我有自己简单的实现方式，但我想知道是否有更好的方式，可以利用Numpy/Scipy模块。 - Navi

也许您可以更新问题，包括（1）您有一个一维数组和（2）您正在寻找什么类型的局部最小值。只是一个比两个相邻条目更小的条目吗？ - Sven Marnach

如果你所说的数据带有噪音，可以看一下scipy.signal.find_peaks_cwt。 - lakshayg

Answer 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 函数，分别用于寻找极大值和极小值。

- danodonovan

---

12的意义是什么？ - marshmallow

np.random.random(12) 生成12个随机值，它们用于演示函数 argrelextrema。请注意，我会尽力使翻译更加通俗易懂，但不会改变原来的意思。 - sebix

如果输入是test02=np.array([10,4,4,4,5,6,7,6])，那么它无法工作。它无法识别连续的值作为局部最小值。 - Leos313

谢谢你，@Cleb。我想指出其他问题：数组的极端点怎么办？第一个元素也是局部最大值，最后一个元素也是局部最小值。而且，它也没有返回找到多少个连续的值。不过，在这个问题的代码中，我提出了一个解决方案在这里。谢谢！ - Leos313

谢谢，这是我迄今为止找到的最佳解决方案之一。 - Noufal E

Answer 2

如果你想要找到一维数组a中所有小于它们相邻元素的条目，可以尝试使用以下代码：

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

在这一步之前，您还可以使用numpy.convolve()平滑您的数组。 链接

我认为没有专门用于此的函数。

- Sven Marnach

---

@Navi：问题在于“局部最小值”的概念因用例而异，因此很难为此目的提供“标准”函数。平滑有助于考虑更多不仅是最近邻居的因素。使用不同的整数（例如3）代替1可能会很奇怪，因为它只考虑两个方向上的第三个元素，而不是直接相邻的元素。 - Sven Marnach

@Sven Marnach：你提供的方法会延迟信号。这里有第二种方法，使用了scipy.signal中的filtfilt函数。具体可以参考这个链接：http://scipy-cookbook.readthedocs.org/items/FiltFilt.html - bobrobbob

仅出于某种目的，将“<”替换为“>”将会给你局部最大值而不是最小值。 - DarkCygnus

@SvenMarnach，我使用了您上面提供的解决方案来解决我在这里发布的问题https://stackoverflow.com/questions/57403659/python-find-multiple-maximum-and-minimums-in-a-list?noredirect=1#comment101288806_57403659，但是我得到了输出`[False False]`。这里可能有什么问题？ - Msquare

@Msquare，问题可能是你的“a”不是一个numpy数组，而是一个Python列表。尝试使用“a = np.array(a)”。 - toliveira

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

- Cleb

---

谢谢你的回答。我想知道，用(-1)乘以输入来找到最小值是否是推荐的方法。在我的脑海里，一直有这样的信念，这不可能是正确的方法。有什么想法吗？ - Martin

@OkLetsdothis：我认为这是相当标准的。在优化问题中，经常使用那个“技巧”；当你试图最大化一个目标函数时，可以将其乘以-1，然后使用最小化方法来解决问题。 - Cleb

Answer 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会减少原始索引号。

- R. C.

---

嵌套的numpy函数使用得很好！但请注意，这会错过数组两端的最大值 :) - danodonovan

如果存在重复的值，这个也会表现得很奇怪。比如说，如果你取数组 [1, 2, 2, 3, 3, 3, 2, 2, 1]，局部最大值显然在中间的三个之间。但是如果你运行提供的函数，你会得到最大值在索引2、6处，最小值在索引1、3、5、7处，对我来说这没有太多意义。 - Korem

为了避免使用“np.diff()”，而应该使用“np.gradient()”来代替，这样可以避免出现“+1”的情况。 - ankostis

我知道这个帖子已经很多年了，但值得一提的是，如果你的曲线太嘈杂，你可以先尝试低通滤波来进行平滑处理。对于我来说，大多数局部最大/最小值的使用都是为了在某个局部区域内找到全局最大/最小值（例如，大的峰和谷，而不是数据中的每一个变化）。 - marcman

请注意，仅仅在上面的代码中删除“+1”并将“np.diff”替换为“np.gradient()”会产生每个极小值/极大值的索引以及它们的最低/最高邻居。 - Reed Espinosa

Answer 5

另一种方法(更多的文字，更少的代码)可能会有所帮助:

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

不幸的是，一阶导数往往会"放大"噪声，因此当原始数据中存在显著的噪声时，最好先对原始数据进行一定程度的平滑处理，然后再使用一阶导数。

由于平滑在最简单的意义上是一个低通滤波器，所以平滑通常最好(也最容易)使用卷积核来完成，并且"塑形"该卷积核可以提供惊人的特征保留/增强能力。找到最佳核的过程可以通过各种手段自动化完成，但最好的方法可能是简单的暴力搜索(对于找到小型核而言足够快)。好的核将(按预期)大量扭曲原始数据，但它不会影响感兴趣的峰/谷的位置。

幸运的是，很多时候可以通过一个简单的SWAG(“educated guess”)创建一个合适的核。平滑核的宽度应该比原始数据中最宽的"有趣"峰略宽一些，其形状将类似于该峰(一个单尺度小波)。对于保持平均值的核(任何好的平滑滤波器都应该是这样)，核元素的总和应该恰好等于1.00，并且该核应该关于其中心对称(意味着它将具有奇数个元素)。

给定一个最佳平滑核(或为不同数据内容优化的少量核)，平滑程度成为(卷积核的"增益")的缩放因子。

确定正确(最佳)的平滑程度(卷积核增益)甚至可以自动化:比较一阶导数数据的标准差与平滑数据的标准差。如何随着平滑程度的变化，两个标准差的比值变化可用于预测有效的平滑值。几次手动数据运行(真正代表性的)应该就足够了。

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

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

希望这提供了足够的信息让Google（和可能是良好的统计文本）填补空白。我真的希望有时间提供一个实例或链接。如果有人在网上找到一个，请在此处发布！

- BobC

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

为了找到局部最大值或最小值，我们主要想找到列表中数值的差异从正变负（最大值）或从负变正（最小值）时的位置。因此，首先我们找到差异，然后找到符号，接着通过再次求差异来找到符号变化的位置。（有点像微积分中的一阶导数和二阶导数，只不过我们有离散数据，没有连续的函数。）

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

因此我们有以下的匹配：

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

- Dave

---

我认为这个（好的！）回答和 R.C. 在 2012 年的回答是一样的？如果我正确理解他的解决方案，他提供了三种单行解决方案，取决于调用者想要最小值、最大值还是两者都想要。 - Brandon Rhodes

它忽略了具有重复元素的情况。例如[1,2,3,1,2,2,2,1,4,5]。如何解决？ - Элёржон Кимсанов

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

敬礼

- A. STEFANI

---

为了避免除法带来的精度损失，为什么不直接将最大值乘以-1，从而实现从极大值到极小值的转换呢？ - Livius

我试图将“1/data”更改为“data * -1”，但是它引发了一个错误，您能分享一下如何实现您的方法吗？ - A. STEFANI

也许是因为我们不想要求最终用户另外安装scipy。 - Damian Yerrick

Answer 8

更新： 我对渐变效果不满意，所以我发现使用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()

- Mike Vella

---

你知道这个梯度是如何计算的吗？如果你有嘈杂的数据，可能梯度会变化很大，但这并不意味着一定存在最大值或最小值。 - Navi

是的，我知道，但是嘈杂的数据是一个不同的问题。为此，我想使用卷积。 - Mike Vella

我在做一个项目时需要类似的东西，所以使用了上面提到的numpy.diff方法。我认为提一下对我的数据来说上述代码错过了一些极大值和极小值可能会有所帮助。通过将两个if语句中的中间项分别改为<=和>=，我能够捕捉到所有的点。 - user723888

Answer 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 

- Misha Smirnov

---

不错的函数！很容易理解。你可能还想看看scipy.signal.find_peaks。它会返回重复组中间的最大值。它还有更多的参数来忽略噪声，可以处理N维数据。https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html - John Henckel

使用了这个算法： - Элёржон Кимсанов

Answer 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包含极小值和极大值的索引。对于一个巨大的数据集，它会给出很多最大值/最小值，因此在这种情况下先平滑曲线，然后再应用该算法。

- prtkp

---

这看起来很有趣。没有库，它是如何工作的？ - john k

沿着曲线从起点开始遍历，观察你是否持续向上或向下移动，一旦从上升变为下降，就意味着你达到了最大值；如果你从下降转为上升，则达到了最小值。 - prtkp

使用了以下方法：`list = [1,2,3,4,5,1,2,3,4,5,0,1,2,2,2,2,2, 1,2,3,1]maxlist = []res = np.diff(list)for i, el in enumerate(res):if el < 0: maxlist.append(list[i-1])` - Элёржон Кимсанов