python - 使用 curve_fit 拟合梯形

Question

我正在尝试使用 scipy.optimize 中的 curve_fit 库将梯形拟合到一组时间序列。我用来生成梯形的函数如下:

def trapezoid(x, a, b, c, tau1, tau2):
    y = np.zeros(len(x))
    c = -np.abs(c)
    a = np.abs(a)
    y[:int(tau1)] =  a*x[:int(tau1)] + b
    y[int(tau1):int(tau2)] =  a*tau1 + b 
    y[int(tau2):] = c*(x[int(tau2):]-tau2) + (a*tau1 + b)
    return y

其中 a 和 c 是斜率，tau1 和 tau2 标记平坦阶段的开始和结束。

为了适合我只使用:

popt, pcov = curve_fit(trapezoid, xdata, ydata, method = 'lm')

对于大多数情况，它工作得很好，如下所示:

但是，我也遇到了一些无法拟合数据的情况，但看起来应该没问题:

这些情况的问题在于，它设置的 tau2(平坦阶段的结束)小于 tau1(平坦阶段的开始)。

有人可以提出解决这个问题的方法吗？是通过施加约束还是以其他方式？

拟合不起作用的示例数组:

数组([1.2, 1.21, 1.2, 1.19, 1.21, 1.22, 2.47, 2.53, 2.49, 2.39, 2.28, 2.16、2.07、1.99、1.91、1.83、1.74、1.65、1.57、1.5、1.45、1.41、 1.38、1.35、1.33、1.29、1.24、1.19、1.14、1.11、1.07、1.04、1.、 0.95、0.91、0.87、0.84、0.8、0.77、0.74、0.72、0.7、0.68、0.66、 0.63、0.61、0.59、0.57、0.55、0.52、0.5、0.48、0.45、0.43、0.41、 0.39, 0.38, 0.37, 0.37, 0.36, 0.35, 0.34, 0.34, 0.33])

其产量:tau1:8.45，tau2:5.99

Answer 1

您可能会发现 lmfit ( http://lmfit.github.io/lmfit-py/ ) 对于此问题很有用。 Lmfit 提供了一个稍微更高级别的曲线拟合接口(interface)，仍然基于 scipy 优化器，但具有一些更好的抽象和功能。

特别是对于您的问题，lmfit 参数是 Python 对象，可以有边界、可以固定或根据其他变量编写为简单的代数约束。这可以支持强加 tau2 > tau1。 这个想法本质上是设置 tau2=tau1+taudiff 并在 taudiff 上设置下限 0。虽然您可以重写函数以在代码中执行此操作，但使用 lmfit 您不必这样做，而是可以将该逻辑放入参数中。

将脚本转换为使用 lmfit 会得到如下所示的结果:

from lmfit import Model

# use your same model function
def trapezoid(x, a, b, c, tau1, tau2):
    y = np.zeros(len(x))
    c = -np.abs(c)
    a = np.abs(a)
    y[:int(tau1)] =  a*x[:int(tau1)] + b
    y[int(tau1):int(tau2)] =  a*tau1 + b 
    y[int(tau2):] = c*(x[int(tau2):]-tau2) + (a*tau1 + b)
    return y

# turn model function into lmfit Model
tmod = Model(trapezoid)

# create Parameters for this model: they will be *named* according
# to the signature of the model function, and be used as keys in
# an ordered-directory-derived object.  Here you can also give
# initial values
params = tmod.make_params(a=1, b=2, c=0.5, tau1=5, tau2=-1)

# now you can set bounds or constraints.  

# 1st, add a new variable "taudiff"
params.add('taudiff', value=0.1, min=0, vary=True)

# constraint tau2 to be taudiff+tau1 -- this is no longer a "free variable:
params['tau2'].expr = "taudiff + tau1"

# now do fit to data:
result = tmod.fit(ydata, params, x=xdata)

# print report of fit
print(result.fit_report())

# get best fit params:
for parname, param in result.params:
    print(parname, param.value, param.stderr, param.expr)

# get best fit array for plotting
pylab.plot(xdata, ydata)
pylab.plot(xdata, result.best_fit)

希望有帮助。