lmfit not exploring parameter space

Question

I'm trying to use lmfit to find the best fit parameters of a function for some random data using the Model and Parameters classes. However, it doesn't seem to be exploring the parameter space very much. It does ~10 function evaluations and then returns a terrible fit.

Here is the code:

import numpy as np
from lmfit.model import Model
from lmfit.parameter import Parameters
import matplotlib.pyplot as plt

def dip(x, loc, wid, dep):
    """Make a line with a dip in it"""
    # Array of ones
    y = np.ones_like(x)

    # Define start and end points of dip
    start = np.abs(x - (loc - (wid/2.))).argmin()
    end = np.abs(x - (loc + (wid/2.))).argmin()

    # Set depth of the dip
    y[start:end] *= dep

    return y

def fitter(x, loc, wid, dep, scatter=0.001, sigma=3):
    """Find the parameters of the dip function in random data"""
    # Make the lmfit model
    model = Model(dip)

    # Make random data and print input values
    rand_loc = abs(np.random.normal(loc, scale=0.02))
    rand_wid = abs(np.random.normal(wid, scale=0.03))
    rand_dep = abs(np.random.normal(dep, scale=0.005))
    print('rand_loc: {}\nrand_wid: {}\nrand_dep: {}\n'.format(rand_loc, rand_wid, rand_dep))
    data = dip(x, rand_loc, rand_wid, rand_dep) + np.random.normal(0, scatter, x.size)

    # Make parameter ranges
    params = Parameters()
    params.add('loc', value=loc, min=x.min(), max=x.max())
    params.add('wid', value=wid, min=0, max=x.max()-x.min())
    params.add('dep', value=dep, min=scatter*10, max=0.8)

    # Fit the data
    result = model.fit(data, x=x, params)
    print(result.fit_report())

    # Plot it
    plt.plot(x, data, 'bo')
    plt.plot(x, result.init_fit, 'k--', label='initial fit')
    plt.plot(x, result.best_fit, 'r-', label='best fit')
    plt.legend(loc='best')
    plt.show()

And then I run:

fitter(np.linspace(55707.97, 55708.1, 100), loc=55708.02, wid=0.04, dep=0.98)

Which returns (for example, since it's randomized data):

rand_loc: 55707.99659784677
rand_wid: 0.02015076619874132
rand_dep: 0.9849809461153651

[[Model]]
    Model(dip)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 9
    # data points      = 100
    # variables        = 3
    chi-square         = 0.00336780
    reduced chi-square = 3.4720e-05
    Akaike info crit   = -1023.86668
    Bayesian info crit = -1016.05117
##  Warning: uncertainties could not be estimated:
    loc:  at initial value
    wid:  at initial value
[[Variables]]
    loc:  55708.0200 (init = 55708.02)
    wid:  0.04000000 (init = 0.04)
    dep:  0.99754082 (init = 0.98)

Any idea why it executes so few function evaluations returning a bad fit? Any assistance with this would be greatly appreciated!

Answer 1

This is a similar question to fitting step function with variation in the step location with scipy optimize curve_fit. See https://stackoverflow.com/a/59504874/5179748.

Basically, the solvers in scipy.optimize/lmfit assume that parameters are continuous -- not discrete -- variables. They make small changes to the parameters to see what change that makes in the result. A small change in your loc and wid parameters will have no effect on the result, as argmin() will always return an integer value.

You might find that using a Rectangle Model with a finite width (see https://lmfit.github.io/lmfit-py/builtin_models.html#rectanglemodel) will be helpful. I changed your example a bit, but it should be enough to get you started:

import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import RectangleModel, ConstantModel

def dip(x, loc, wid, dep):
    """Make a line with a dip in it"""
    # Array of ones
    y = np.ones_like(x)

    # Define start and end points of dip
    start = np.abs(x - (loc - (wid/2.))).argmin()
    end = np.abs(x - (loc + (wid/2.))).argmin()

    # Set depth of the dip
    y[start:end] *= dep
    return y

x = np.linspace(0, 1, 201)
data = dip(x, 0.3, 0.09, 0.98) + np.random.normal(0, 0.001, x.size)

model = RectangleModel() + ConstantModel()
params = model.make_params(c=1.0, amplitude=-0.01, center1=.100, center2=0.7, sigma1=0.15)

params['sigma2'].expr = 'sigma1' # force left and right widths to be the same size
params['c'].vary = False         # force offset = 1.0 : value away from "dip"


result = model.fit(data, params, x=x)
print(result.fit_report())

plt.plot(x, data, 'bo')
plt.plot(x, result.init_fit, 'k--', label='initial fit')
plt.plot(x, result.best_fit, 'r-', label='best fit')
plt.legend(loc='best')
plt.show()