使用python修复拐点估计答案

【问题标题】：Fixing inflexion point estimate using python使用python修复拐点估计
【发布时间】：2016-03-08 02:58:36
【问题描述】：

我正在尝试使用 python 找到曲线上的拐点。曲线的数据在这里：https://www.dropbox.com/s/rig8frgewde8i5n/fitted.txt?dl=0。请注意，曲线已与原始数据拟合。原始数据可在此处获得：https://www.dropbox.com/s/1lskykdi1ia1lu7/ww.txt?dl=0

import numpy as np
# Read in array from text file
arr = np.loadtxt(path_to_file)

inflexion_point_1 = np.diff(arr).argmin()
inflexion_point_2 = np.diff(arr).argmax()

这两个拐点在附图中显示为红线。但是，它们的位置似乎并不正确。第一个拐点应靠近黑色箭头指示的区域。我该如何解决这个问题？

另外，这里是一个微分图：

plt.axvline(np.gradient(arr[:365]).argmax())

如您所见，代码的行为与编码相同，即它找到了数组的 np.diff 的 argmax。但是，我想找到一个接近第 110 天左右的位置，即大约到 argmax 的一半。

--编辑--

另外，这里是另一个显示原始数据和拟合曲线的图（使用二次函数）。

【问题讨论】：

我会为您的数据拟合一个函数并获得该函数的导数。由于您的步骤，将其直接应用于您的数据似乎会导致问题。
您使用什么基础模型进行拟合？这对我来说显然是过度拟合。你能修改你的图，让我们在图中看到拟合所基于的实际数据点吗？！
好的，所以它实际上并不合适，只是平滑！？我仍然会尝试将实际函数拟合到蓝线，然后取它的导数。
@BobBaxley 对于平滑函数，拐点也可以表征为一阶导数的局部最大值或局部最小值。如果数据良好且平滑，np.diff(arr).argmax() 和np.diff(arr).argmin() 将是估计拐点的合理方法。但是，这种情况下的数据是有噪声的，所以简单的有限差分不能很好地工作。
@WarrenWeckesser 好点子和好收获。

标签： python numpy scipy

【解决方案1】：

有什么理由不在渐变上直接使用单变量样条？

from scipy.interpolate import UnivariateSpline

#raw data
data = np.genfromtxt('ww.txt')

plt.plot(np.gradient(data), '+')

spl = UnivariateSpline(np.arange(len(data)), np.gradient(data), k=5)
spl.set_smoothing_factor(1000)
plt.plot(spl(np.arange(len(data))), label='Smooth Fct 1e3')
spl.set_smoothing_factor(10000)
plt.plot(spl(np.arange(len(data))), label='Smooth Fct 1e4')
plt.legend(loc='lower left')

max_idx = np.argmax(spl(np.arange(len(data))))
plt.vlines(max_idx, -5, 9, linewidth=5, alpha=0.3)

我们也可以求解最大值：

In [122]:

import scipy.optimize as so
F = lambda x: -spl(x)
so.fmin(F, 102)
Optimization terminated successfully.
         Current function value: -3.339112
         Iterations: 20
         Function evaluations: 40
Out[122]:
array([ 124.91303558])

【讨论】：