Search code examples
pythonnumpyscipy-optimizelmfit

lmfit does not fit the only parameter on a simple example that returns a single scalar value

A residue function (res) does the sum of y values for y>thr (threshold). It returns the residue between this sum and the target.

In this example, a target is calculated for y>70 and I would like to find, after minimization, y=70 starting from y=60.

import numpy as np
import lmfit

x=np.linspace(0,10,51)
y=(x+4)*(x-14)*-1

def res(params,y,target):
    parvals = params.valuesdict()
    sum=y[y>parvals['thr']].sum()
    return [(target-sum)**2]

target=y[y>70].sum()

pars=lmfit.Parameters()
pars.add('thr', value=60, min=y.min(), max=y.max())

minner = lmfit.Minimizer(res, pars, fcn_args=(y, target))
result = minner.minimize()

Why does not the fit work : 70 is not returned but 60 (the initial value) is returned ?

Thanks for answer.

Solution

  • The fit does not work because you are using a continuous variable (pars['thr']) as a discrete value [y>parvals['thr']]. When the fit is running, it will try to calculate the change in the result by making small changes in the variable values. If you add a print() to your function:

    def res(params,y,target):
    parvals = params.valuesdict()
    print(parvals['thr'])
    sum=y[y>parvals['thr']].sum()
    return [(target-sum)**2]

    you will get

    60.0
60.0
60.0
60.00000089406967
60.0

    As the fit tries to find which way and by how much thr should change. And it will see that changing thr has no impact on the result.

    Basically, you need to change your use of thr from being a discrete variable into being a continuous variable. One easy way to do this is with the erf function (or another sigmoidal function), scaling to go from ~0 to ~1 over some small but not infinitesimal interval, for example

    from scipy.special import erf

def res(params,y,target):
    parvals = params.valuesdict()   
    step = (1 + erf(y-parvals['thr']))/2.0
    sum = (y*step).sum()
    return target-sum   # was [(target-sum)**2]

    That step array will have non-zero/non-one values near the threshold giving just enough signal for the fit to decide where to move.

    Also, note that returning [(target-sum)**2] will work, but that just returning the residual target-sum will allow the fit to see not only the magnitude but also the sign of the misfit, and allow the fit to converge faster.

    With those changes, you should get the correct value for thr of 70 with about 11 function evaluations.