Piecewise Function lmfit

Question

I am trying to define a piecewise function to be fitted by lmfit library in Python. The issue I am having is a parameter I have defined for the function will not evaluate alongside the data I am submitting.

I have one example of a case somewhat similar to mine here. However, the vectorize function the answer describes wasn't producing values I wanted, and when reading the documentation, it didn't seem to be the answer to my solution. I also used scipy.optimize.leastsq, but I got the same issue with lmfit described below.

I have a my residual function defined such as

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))

I also tried defining the function such that

  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

I have also tried inline using a lambda function.

I am getting an error 'The truth value of an array with more than one element is ambiguous.' When I got the error, it made sense why it happened, because (param2 < x) would produce a logical array. However, I can't seem to find a way to define the function in a piecewise fashion with the given case to get it fitted with the lmfit.minimize() function. I have seen the answer done in Matlab, in which it's nlinfit function seems to evaluate the data element-wise without issue (I tried searching if Python has an equivalent operation to define element-wise computation such as .* or .+, but that doesn't seem to exist as explicitly).

lmfit also seems to operate a bit differently compared to nlinfit, because we have to always have our residuals return (model - y) while nlinfit outputs the result once the function is given, which I am not sure could be another issue.

So to reiterate, my main question is if there is a method of defining the piecewise function such that it can compare the parameter to the data set.

Any help or explanation would be appreciated, thank you!

Answer 1

In place of (param2 < x) (where param2 is a float and x is an numpy array), you want to use numpy.where. You might try:

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

I should also warn you about a potential problem with this approach to having a variable be a boundary for a piecewise function.

In non-linear fits, variables are always floating point (continuous, non-discrete) values. As the fit proceeds, it will make small adjustments in the values and see how that small change alters the result. In your approach, the parameter 'two' is used as both the transition between pieces and the offset for the line -- that is good.

If a parameter is used only as the transition, it may not work. Consider, say, x=np.array([0, 1., 2., 3., 4., ..., 20.0]). Having two = 10.5 and two=10.4 would then give the same result. In that case, the fit would not be able to alter the value of two: it would try a very small change, see no change in the result and give up.

So, either make sure that two is also used elsewhere in your real model (assuming your real model is more complicated than the example given), or consider using a more gentle transition rather than a hard change in pieces. I find an error-function of width ~spacing between x points often works. Depending on the nature of your problem, you might try something like this:

 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

Hope that helps.