python - 有没有一种方法可以在scipy.optimize.curve_fit的参数之间添加依赖项

我试图描述海洋模型中由70年代核试验引起的碳同位素扩散的特征。
大气信号是一个强烈的尖峰，它将随洋流传播到更深的地方（深流慢得多）。

我的目标是检测各种深度级别的浓度升高和增长率。

我假设海洋中的碳同位素浓度表现为具有3个分段的分段线性函数：


直到时间（b）的恒定初始值（t_0）
浓度从时间（t_0）到（t_1）随速率m1线性增加。
时间（t_1）之后浓度随速率m2线性下降


我在python中使用以下代码表示函数：

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as sio

def piecewise_linear( t, t0, t1, b, m1, m2 ):
    condlist = [ t < t0,
                (t >= t0 ) & ( t < t1 ),
                t >= t1
               ]
    funclist = [lambda t: b,
                lambda t: b + m1 * ( t - t0 ),
                lambda t: b + m1 * ( t - t0 ) + m2 * ( t - t1 )
               ]
    return np.piecewise( t, condlist, funclist )

对于给定的时间数组t，我希望能够适合此函数的两个“类型”：


完整的3段线代表了上层海洋，信号在该区域快速传播，峰值被完全捕获。
一种特殊情况，在时间序列结束时浓度没有达到峰值（这代表深海中的信号，在该信号中传播信号需要很长时间）


举个例子

t = np.arange( 0, 15, 0.1 )
y_full = piecewise_linear( t, 5, 10, 2, 2, -4 )
y_cut = piecewise_linear( t, 5, 15, 2, 2, -4 )
plt.plot( t, y_full )
plt.plot( t, y_cut )
plt.legend( [ 'surface', 'deep ocean' ] )

对于第一种情况，当我在添加一些随机噪声后尝试将函数拟合到信号时，我得到了很好的结果：

noise = np.random.normal( 0, 1, len( y_full ) ) * 1
y = y_full
yy = y_full + noise
bounds = ( [ 0, 0, 0, 0, -np.inf ], [ np.inf, np.inf, np.inf, np.inf, 0 ] )
fit,_ = sio.curve_fit( piecewise_linear, t, yy, bounds=bounds )
print( fit )
y_fit = piecewise_linear( t, *tuple( fit ) )
plt.plot( t, yy, color='0.5' )
plt.plot( t, y_fit, linewidth=3 )
plt.plot( t, y, linestyle='--', linewidth=3 )

导致

>>[  5.00001407  10.01945313   2.13055863   1.95208167  -3.95199719]

但是，当我尝试评估第二种情况（深海）时，通常会得到如下不良结果：

noise = np.random.normal( 0, 1, len(y_full ) ) * 1#
y = y_cut
yy = y_cut+noise
bounds = ( [ 0, 0, 0, 0, -np.inf], [ np.inf, np.inf, np.inf, np.inf, 0 ] )
fit,_ = sio.curve_fit( piecewise_linear, t, yy, bounds=bounds )
print( fit )
y_fit = piecewise_linear( t, *tuple( fit ) )
plt.plot( t, yy, color='0.5' )
plt.plot( t, y_fit, linewidth=3 )
plt.plot( t, y, linestyle='--', linewidth=3 )
plt.legend( [ 'noisy data', 'fit', 'original' ] )

我懂了

>>[ 1.83838997  0.40000014  1.51810839  2.56982348 -1.0622842 ]

优化确定t_0大于t_1，在这种情况下这是没有意义的。

有没有办法将条件t_0 < t_1建立到曲线拟合中？还是我必须测试给出哪种类型的曲线，然后将其拟合到两个不同的函数（3段或2段分段线性函数）？

任何帮助是极大的赞赏

最佳答案

您可以考虑为此使用lmfit（https://lmfit.github.io/lmfit-py）。
Lmfit为曲线拟合提供了更高级别的接口，并使拟合参数成为一流的python对象。除其他事项外，这很容易修复某些参数，并以比scipy.optimize.curve_fit所使用的Python语言更Python化的方式设置参数的界限。特别是对于您的问题，lmfit参数还支持使用数学表达式作为所有参数的约束表达式。

要将模型函数piecewise_linear()转换为模型以使用lmfit进行曲线拟合，您需要执行以下操作

from lmfit import Model

# make a model
mymodel = Model(piecewise_linear)

# create parameters and set initial values
# note that parameters are *named* from the
# names of arguments of your model function
params = mymodel.make_params(t0=0, t1=1, b=3, m1=2, m2=2)

# now, you can place bounds on parameters, maybe like
params['b'].min = 0
params['m1'].min = 0

# but what you want is an inequality constraint, so
#   1. add a new parameter 'tdiff'
#   2. constrain t1 = t0 + tdiff
#   3. set a minimum value of 0 for tdiff
params.add('tdiff', value=1, min=0)
params['t1'].expr = 't0 + tdiff'

# now perform the fit
result = mymodel.fit(yy, params, t=t)

# print out results
print(result.fit_report())

您可以在lmfit文档或其他SO问题中阅读如何从拟合结果中提取其他信息。