用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（它的开始）。

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

不适合的示例数组：

array([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

只需将 t1,t2 设置为最小值和最大值即可

def trapezoid(x, a, b, c, tau1, tau2):
    y = np.zeros(len(x))
    c = -np.abs(c)
    a = np.abs(a)
    (tau1,tau2) = (min(tau1,tau2),max(tau1,tau2))

    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)



x_data = np.arange(len(A))
popt, pcov = curve_fit(trapezoid, x_data, A, method = 'lm')
print popt
fit = trapezoid(x_data,*popt)

导致：

Answer 2

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

特别是对于您的问题，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)

希望对您有所帮助。