改进算法以找到一维 numpy 数组的局部最大值答案

【问题标题】：Improving algorithm to find local maxima of 1D numpy array改进算法以找到一维 numpy 数组的局部最大值
【发布时间】：2018-02-05 20:41:36
【问题描述】：

我试图找到函数f(x) = (sin(x)/x)^2 的局部最大值。对于近似解决方案，我初始化了两个变量x 和y，并首先绘制了一个图形以进行可视化表示。

x = np.linspace(-20.0, 20.0, num=int((20+20)/0.01)) 
y = np.power(np.sin(x)/x, 2)
plt.plot(x, y, 'r.', markersize= 1)
plt.show()

这显示graph。

然后我尝试创建一个算法来找到千里马：

def returnMaxima(num, x, y):
    """
    number, np.array, np.array -> list
    num: number of maxima needed | x: x 1D array | y: y 1D array
    returns [[x1,y1], [x2,y2]...] in descending order of y
    """
    allMaximaPoints = [] # stores all Maxima points
    reqMaximaPoints = [] # stores num Maxima points
    for i in range(y.size): 
        # for first y value
        if i == 0: 
            if y[i] > y[i+1]:
                allMaximaPoints += [[x[i], y[i]], ]
        # for last y value
        elif i == y.size - 1:
            if y[i] > y[i-1]:
                allMaximaPoints += [[x[i], y[i]], ]
        # for rest y values
        else: 
            if y[i] > y[i-1] and y[i] > y[i+1]:
                allMaximaPoints += [[x[i], y[i]], ]
    # extract largest maximas from allMaximaPoints
    while num > 0: 
        reqMaximaPoints += [max(allMaximaPoints, key=lambda item:item[1]),]
        del allMaximaPoints[allMaximaPoints.index(max(allMaximaPoints, key=lambda item:item[1]))]
        num -= 1
    return reqMaximaPoints

当我尝试returnMaxima(2, x, y) 时，我得到[[-4.4961240310077528, 0.04719010162459622], [4.4961240310077528, 0.04719010162459622]]。

这是不正确的，因为它跳过了 x = 0 处的局部最大值。我怀疑这是因为在x=0 处与y[i] 的最大值相邻的y[i-1] 和y[i+1] 值大约等于y[i] 导致代码

else:
    if y[i] > y[i-1] and y[i] > y[i+1]:
        allMaximaPoints += [[x[i], y[i]], ]

不考虑这一点。这是因为当我将x = np.linspace(-20.0, 20.0, num=int((20+20)/0.01)) 更改为x = np.linspace(-20.0, 20.0, num=int((20+20)/0.1)) 即x 中的较大步长时，正确找到了x=0 处的局部最大值。但是，即使我将上述代码中的> 符号更改为>=，x=0 处的最大值仍然不计算在内。

为什么会这样？我应该如何改进我的代码以获得正确的结果？谢谢！

【问题讨论】：

我为 maxima 找到的最好的库是 peakutils: pypi.python.org/pypi/PeakUtils 。它允许您选择灵敏度和阈值以使其对噪声具有鲁棒性。

标签： python arrays algorithm numpy

【解决方案1】：

使用scipy.signal.find_peaks_cwt 之类的东西可能会更好。比如：

indices = scipy.signal.find_peaks_cwt(y, [1, 2, 3, 4], noise_perc=50)
plt.plot(x, y)
plt.plot(x[indices], y[indices], 'r.')
plt.show()

结果：

【讨论】：

【解决方案2】：

简答：x 永远不会是 0（您可能也不希望它也是）。

我们可以通过运行来测试这个

In [53]: 0 in x
Out[53]: False

这实际上是幸运的，因为 x=0 会导致运行时警告并由于除以 0 而产生“nan”（非数字）值。问题是函数未定义为 0。

否则，您对它为什么不起作用的原因基本上是正确的。如果您将> 换成>=，您会看到输出不同：

In [55]: returnMaxima(2, x, y)
Out[55]:
[[-0.0050012503125778096, 0.99999166252624228],
[0.0050012503125778096, 0.99999166252624228]]

这实际上是“正确”的输出，因为这些是最大值，并且 x 的值最接近 0。

【讨论】：