这就是我如何检测局部最大值/最小值、拐点和鞍座的方法。
首先定义以下函数
import numpy as np
def n_derivative(arr, degree=1):
"""Compute the n-th derivative."""
result = arr.copy()
for i in range(degree):
result = np.gradient(result)
return result
def sign_change(arr):
"""Detect sign changes."""
sign = np.sign(arr)
result = ((np.roll(sign, 1) - sign) != 0).astype(bool)
result[0] = False
return result
def zeroes(arr, threshold=1e-8):
"""Find zeroes of an array."""
return sign_change(arr) | (abs(arr) < threshold)
我们现在可以利用
derivative test
临界点的一阶导数等于零。
def critical_points(arr):
return zeroes(n_derivative(arr, 1))
如果临界点具有非零的二阶导数,则该点为最大值或最小值:
def maxima_minima(arr):
return zeroes(n_derivative(arr, 1)) & ~zeroes(n_derivative(arr, 2))
def maxima(arr):
return zeroes(n_derivative(arr, 1)) & (n_derivative(arr, 2) < 0)
def minima(arr):
return zeroes(n_derivative(arr, 1)) & (n_derivative(arr, 2) > 0)
如果二阶导数等于零,但三阶导数非零,则该点为拐点:
def inflections(arr):
return zeroes(n_derivative(arr, 2)) & ~zeroes(n_derivative(arr, 3))
如果临界点的二阶导数等于零,但三阶导数非零,则这是一个鞍:
def inflections(arr):
return zeroes(n_derivative(arr, 1)) & zeroes(n_derivative(arr, 2)) & ~zeroes(n_derivative(arr, 3))
注意,该方法在数值上不稳定,一方面在某些任意阈值定义上检测到零,另一方面不同的采样可能导致函数/数组不可微。
因此,根据这个定义,您期望的实际上不是鞍点。
为了更好地逼近连续函数,可以在很大程度上过采样(根据
K
在代码)函数中,例如:
import scipy as sp
import scipy.interpolate
data = [
1.04814804, 0.90445908, 0.62026396, 0.60566623, 0.32295758, 0.26658469, 0.19059289,
0.10281547, 0.08582772, 0.05091265, 0.03391474, 0.03844931, 0.03315003, 0.02838656,
0.03420759, 0.03567401, 0.038203, 0.03530763, 0.04394316, 0.03876966, 0.04156067,
0.03937291, 0.03966426, 0.04438747, 0.03690863, 0.0363976, 0.03171374, 0.03644719,
0.02989291, 0.03166156, 0.0323875, 0.03406287, 0.03691943, 0.02829374, 0.0368121,
0.02971704, 0.03427005, 0.02873735, 0.02843848, 0.02101889, 0.02114978, 0.02128403,
0.0185619, 0.01749904, 0.01441699, 0.02118773, 0.02091855, 0.02431763, 0.02472427,
0.03186318, 0.03205664, 0.03135686, 0.02838413, 0.03206674, 0.02638371, 0.02048122,
0.01502128, 0.0162665, 0.01331485, 0.01569286, 0.00901017, 0.01343558, 0.00908635,
0.00990869, 0.01041151, 0.01063606, 0.00822482, 0.01312368, 0.0115005, 0.00620334,
0.0084177, 0.01058152, 0.01198732, 0.01451455, 0.01605602, 0.01823713, 0.01685975,
0.03161889, 0.0216687, 0.03052391, 0.02220871, 0.02420951, 0.01651778, 0.02066987,
0.01999613, 0.02532265, 0.02589186, 0.02748692, 0.02191687, 0.02612152, 0.02309497,
0.02744753, 0.02619196, 0.02281516, 0.0254296, 0.02732746, 0.02567608, 0.0199178,
0.01831929, 0.01776025]
samples = np.arange(len(data))
f = sp.interpolate.interp1d(samples, data, 'cubic')
K = 10
N = len(data) * K
x = np.linspace(min(samples), max(samples), N)
y = f(x)
然后,所有这些定义都可以通过以下方式进行视觉测试:
import matplotlib.pyplot as plt
plt.figure()
plt.plot(samples, data, label='data')
plt.plot(x, y, label='f')
plt.plot(x, n_derivative(y, 1), label='d1f')
plt.plot(x, n_derivative(y, 2), label='d2f')
plt.plot(x, n_derivative(y, 3), label='d3f')
plt.legend()
for w in np.where(inflections(y))[0]:
plt.axvline(x=x[w])
plt.show()
但即使在这种情况下,这一点也不是马鞍。
