分段函数 lmfit

Piecewise Function lmfit

我正在尝试定义一个分段函数以供 Python 中的 lmfit 库拟合。我遇到的问题是我为函数定义的参数不会与我提交的数据一起评估。

我有一个例子和我的有点相似here。但是,答案描述的矢量化函数并没有产生我想要的值,并且在阅读文档时,它似乎不是我解决方案的答案。我也使用了 scipy.optimize.leastsq,但我在下面描述的 lmfit 中遇到了同样的问题。

我定义了一个残差函数,例如

from lmfit import minimize, Parameters, Model

def residual(params, y, x):
    param1 = params['one']
    param2 = params['two']
    if(param2 < x):
        p = 1
    else:
        p = param1*x + param2
    return p - y 

params = Parameters()
params.add('one', value=1)
params.add('two', value=2)
out = minimize(residual, params,args=(y,x))

我也试过这样定义函数

  def f(param1,param2,x):
    if(param2 < x):
        p = 1
    else:
        p = param1*x + param2
    return p

  def residual(params, y, x):
    param1 = params['one']
    param2 = params['two']
    return f(param1,param2,x) - y

我也尝试过使用 lambda 函数进行内联。

我收到一个错误 'The truth value of an array with more than one element is ambiguous.' 当我收到错误时,它发生的原因是有道理的,因为 (param2 < x) 会生成一个逻辑数组。但是,我似乎无法找到一种方法来使用给定的案例以分段方式定义函数以使其适合 lmfit.minimize() 函数。我已经看到在 Matlab 中完成的答案,其中它的 nlinfit 函数似乎可以毫无问题地按元素评估数据(我尝试搜索 Python 是否具有等效的操作来定义按元素计算,例如 .* 或 . +,但这似乎并不明确存在)。

lmfit 与 nlinfit 相比,操作似乎也有点不同,因为我们必须始终有残差 return(模型 - y),而 nlinfit 在给出函数后输出结果,而我不是当然可能是另一个问题。

重申一下,我的主要问题是是否有一种方法可以定义分段函数,以便将参数与数据集进行比较。

如有任何帮助或解释,我们将不胜感激,谢谢!

代替 (param2 < x)(其中 param2 是一个浮点数,x 是一个 numpy 数组),您想使用 numpy.where。你可以试试:

def residual(params, y, x):
    param1 = params['one']
    param2 = params['two']
    p = param1 * x + param2
    p[np.where(param2 < x)] = 1.0 
    return p - y

我还应该警告您这种将变量作为分段函数边界的方法可能存在的问题。

在非线性拟合中,变量始终是浮点(连续、非离散)值。随着拟合的进行,它将对值进行小的调整,并查看该小的变化如何改变结果。在您的方法中,参数 'two' 用作片段之间的过渡和线的偏移量——这很好。

如果一个参数只用作转换,它可能不起作用。考虑一下,比如说,x=np.array([0, 1., 2., 3., 4., ..., 20.0])two = 10.5two=10.4 会给出相同的结果。在那种情况下,拟合将无法更改 two 的值:它会尝试一个非常小的更改,看到结果没有变化并放弃。

因此,要么确保 two 也在您的真实模型的其他地方使用(假设您的真实模型比给定的示例更复杂),要么考虑使用更温和的过渡而不是硬更改成碎片。我发现 x 点之间宽度 ~spacing 的误差函数通常有效。根据问题的性质,您可以尝试这样的操作:

 from scipy.special import erf, erfc
 def residual(params, y, x):
    param1 = params['one']
    param2 = params['two']
    dx = (max(x) - min(x))/(len(x)-1)
    xhi = (erf((x-param2)/dx) + 1)/2.0
    xlo = (erfc((x-param2)/dx) + 1)/2.0
    p = xlo*1.0 + xhi*(param1*x + param2)
    # note: did you really want?
    # p = xlo*param + xhi*(param1*x + param2)
    # p = param2 + xhi*param1*x
    return p - y

希望对您有所帮助。