分段函数lmfit

Question

我试图定义一个分段函数通过在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函数也试过内联。

我得到一个错误“数组的一个以上的要素是模糊的真值。”当我得到的错误，这是有意义为什么会发生，因为（参数2 <X）将产生一个逻辑阵列。不过，我似乎无法找到一种方法，以分段方式与给定的情况下，把它装上lmfit.minimize（）函数来定义函数。我已经看到了在Matlab做了回答，其中的nlinfit功能似乎来评估数据元素与明智的没有问题（我试图寻找Python是否具有等效的操作定义元素方面的计算，如*或+，但似乎并不为明确存在）。

lmfit也似乎操作有些不同相比nlinfit，因为我们总是有我们的残差回报（模型 - Y）的一次函数给出，这我不知道可能是另外一个问题，而nlinfit输出结果。

因此，要重申，我的主要问题是，如果有定义分段函数，使得它可以比较参数数据集的方法。

任何帮助或解释，将不胜感激，谢谢！

Answer 1

在地方(param2 < x)（其中param2是float和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

我也应该提醒你注意这种方法具有变量是一个分段函数的边界的潜在问题。

在非线性拟合，变量总是浮点（连续的，非离散）的值。作为配合进行，它将使中的值小的调整，并看到小的变化如何改变的结果。在你的方法，参数“两节”作为两片之间的过渡为线偏移 - 这是好的。

如果参数只能作为过渡，它可能无法正常工作。考虑，比如说，x=np.array([0, 1., 2., 3., 4., ..., 20.0])。然后有two = 10.5和two=10.4会给出相同的结果。在这种情况下，适合将无法改变two的价值：它会尝试一个非常小的变化，看到在结果没有变化，放弃。

所以，要么确保two也在自己的真实模型中使用的其他地方（假设你的真实模型比给出的例子更复杂），或者考虑使用更温和的过渡，而不是件很难改变。我觉得宽〜x点之间的间距往往工作的误差函数。根据您的问题的性质，你可以尝试这样的事：

 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

希望帮助。