Plotting and modeling data with lmfit - Fit doesn't match data. What am I doing wrong?

Question

I have some data I'm trying to model with lmfit's Model.

Specifically, I'm measuring superconducting resistors. I'm trying fit the experimental data (resistance vs. temperature) to a model which incorporates the critical temperature Tc (material dependent), the resistance below Tc (nominally 0), and the resistance above Tc (structure dependent).

Here's a simplified version (with simulated data) of the code I'm using to plot my data, along with the output plot.

I'm not getting any errors but, as you can see, I'm also not getting a fit that matches my data.

What am I doing wrong? This is my first time using lmfit and Model, so I may be making a newbie mistake. I thought I was following the lmfit example but, as I said, I'm obviously doing something wrong.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from lmfit import Model


def main():
    x = np.linspace(0, 12, 50)
    x_ser = pd.Series(x)  # Simulated temperature data

    y1 = [0] * 20
    y2 = [10] * 30
    y1_ser = pd.Series(y1)  # Simulated resistance data below Tc
    y2_ser = pd.Series(y2)  # Simulated resistance data above Tc (
    y_ser = y1_ser.append(y2_ser, ignore_index=True)

    xcrit_model = Model(data_equation)
    params = xcrit_model.make_params(y1_guess=0, y2_guess=12, xcrit_guess=9)
    print('params: {}'.format(params))
    result = xcrit_model.fit(y_ser, params, x=x_ser)
    print(result.fit_report())

    plt.plot(x_ser, y_ser, 'bo', label='simulated data')
    plt.plot(x_ser, result.init_fit, 'k.', label='initial fit')
    plt.plot(x_ser, result.best_fit, 'r:', label='best fit')
    plt.legend()
    plt.show()


def data_equation(x, y1_guess, y2_guess, xcrit_guess):
    x_lt_xcrit = x[x < xcrit_guess]
    x_ge_xcrit = x[x >= xcrit_guess]
    y1 = [y1_guess] * x_lt_xcrit.size
    y1_ser = pd.Series(data=y1)
    y2 = [y2_guess] * x_ge_xcrit.size
    y2_ser = pd.Series(data=y2)

    y = y1_ser.append(y2_ser, ignore_index=True)

    return y

if __name__ == '__main__':
    main()

Answer 1

lmfit (and basically all similar solvers) work with continuous variables and investigate how they alter the result by making tiny changes in the parameter values and seeing how that effects this fit.

But your xcrit_guess parameter is used only as a discrete variable. If its value changes from 9.0000 to 9.00001, the fit will not change at all.

So, basically, don't do:

x_lt_xcrit = x[x < xcrit_guess]
x_ge_xcrit = x[x >= xcrit_guess]

Instead, you should use a smoother sigmoidal step function. In fact, lmfit has one of these built-in. So you might try something like this (note, there is no point in converting numpy.arrays to pandas.Series - the code will just turn these back to numpy arrays anyway):

import numpy as np
from lmfit.models import StepModel
import matplotlib.pyplot as plt

x = np.linspace(0, 12, 50)
y = 9.5*np.ones(len(x))
y[:26] = 0.0
y = y + np.random.normal(size=len(y), scale=0.0002)

xcrit_model = StepModel(form='erf')
params = xcrit_model.make_params(amplitude=4, center=5, sigma=1)

result = xcrit_model.fit(y, params, x=x)
print(result.fit_report())

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