Lorentzian + background with lmfit

Introduction

The objective of this notebook is to show how to combine the models of the QENSlibrary. Here, we use the Lorentzian profile and a flat background, created from background_polynomials, to perform some fits.

lmfit is used for fitting.

Physical units

For information about unit conversion, please refer to the jupyter notebook called Convert_units.ipynb in the tools folder.

The dictionary of units defined in the cell below specify the units of the refined parameters adapted to the convention used in the experimental datafile.

# Units of parameters for selected QENS model and experimental data
dict_physical_units = {'omega': "1/ps",
                       'scale': "unit_of_signal/ps",
                       'center': "1/ps",
                       'hwhm': "1/ps",}

Import libraries

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets
import lmfit
import QENSmodels

Plot fitting model

The widget below shows the lorentzian peak shape function with a constant background imported from QENSmodels where the functions’ parameters Scale, Center, FWHM and background can be varied.

# Dictionary of initial values
ini_parameters = {'scale': 5, 'center': 0, 'hwhm': 3, 'background': 0.}

def interactive_fct(scale, center, hwhm, background):
    """
    Plot to be updated when ipywidgets sliders are modified
    """
    xs = np.linspace(-10, 10, 100)

    fig1, ax1 = plt.subplots()
    ax1.plot(xs,
             QENSmodels.lorentzian(xs, scale, center, hwhm) +\
             QENSmodels.background_polynomials(xs, background))
    ax1.set_xlabel('x')
    ax1.grid()

# Define sliders for modifiable parameters and their range of variations

scale_slider = ipywidgets.FloatSlider(value=ini_parameters['scale'],
                                      min=0.1, max=10, step=0.1,
                                      description='scale',
                                      continuous_update=False)

center_slider = ipywidgets.IntSlider(value=ini_parameters['center'],
                                     min=-10, max=10, step=1,
                                     description='center',
                                     continuous_update=False)

hwhm_slider = ipywidgets.FloatSlider(value=ini_parameters['hwhm'],
                                     min=0.1, max=10, step=0.1,
                                     description='hwhm',
                                     continuous_update=False)

background_slider = ipywidgets.FloatSlider(value=ini_parameters['background'],
                                       min=0.1, max=10, step=0.1,
                                       description='background',
                                       continuous_update=False)

grid_sliders = ipywidgets.HBox([ipywidgets.VBox([scale_slider, center_slider]),
                                ipywidgets.VBox([hwhm_slider, background_slider])])

# Define function to reset all parameters' values to the initial ones
def reset_values(b):
    """
    Reset the interactive plots to inital values
    """
    scale_slider.value = ini_parameters['scale']
    center_slider.value = ini_parameters['center']
    hwhm_slider.value = ini_parameters['hwhm']
    background_slider.value = ini_parameters['background']


# Define reset button and occurring action when clicking on it
reset_button = ipywidgets.Button(description = "Reset")
reset_button.on_click(reset_values)

# Display the interactive plot
interactive_plot = ipywidgets.interactive_output(interactive_fct,
                                         {'scale': scale_slider,
                                          'center': center_slider,
                                          'hwhm': hwhm_slider,
                                          'background': background_slider})

display(grid_sliders, interactive_plot, reset_button)

Create the reference data

# Create array of reference data: noisy lorentzian with background
nb_points = 100
xx = np.linspace(-5, 5, nb_points)
added_noise = np.random.normal(0, 1, nb_points)
lorentzian_noisy = QENSmodels.lorentzian(
    xx,
    scale=0.89,
    center=-0.025,
    hwhm=0.45
) * (1 + 0.1 * added_noise) + 0.5 * (1 + 0.02 * added_noise)

Setting and fitting

def flat_background(x, A0):
    """
    Define flat background to be added to fitting model
    """
    return QENSmodels.background_polynomials(x, A0)

gmodel = lmfit.Model(QENSmodels.lorentzian) + lmfit.Model(flat_background)
print(f'Names of parameters: {gmodel.param_names}\nIndependent variable(s): {gmodel.independent_vars}')

initial_parameters_values = {'scale': 1, 'center':0.2, 'hwhm': 0.5, 'A0': 0.33}

# Fit
result = gmodel.fit(
    lorentzian_noisy,
    x=xx,
    scale=initial_parameters_values['scale'],
    center=initial_parameters_values['center'],
    hwhm=initial_parameters_values['hwhm'],
    A0=initial_parameters_values['A0']
)

Names of parameters: ['scale', 'center', 'hwhm', 'A0']
Independent variable(s): ['x']

# Plot initial model and reference data
fig0 = plt.figure()
plt.plot(xx, lorentzian_noisy, 'b-', label='reference data')
plt.plot(xx, result.init_fit, 'k--', label='model with initial guesses')
plt.xlabel('x')
plt.title('Initial model and reference data')
plt.grid()
plt.legend();

Plot results

using methods implemented in lmfit

# display result
print('Result of fit:\n', result.fit_report())

# plot fitting result using lmfit functionality
result.plot()

Result of fit:
 [[Model]]
    (Model(lorentzian) + Model(flat_background))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 31
    # data points      = 100
    # variables        = 4
    chi-square         = 0.08314541
    reduced chi-square = 8.6610e-04
    Akaike info crit   = -701.233444
    Bayesian info crit = -690.812764
    R-squared          = 0.96564274
[[Variables]]
    scale:   0.91785782 +/- 0.02883869 (3.14%) (init = 1)
    center: -0.02611501 +/- 0.01023971 (39.21%) (init = 0.2)
    hwhm:    0.43150919 +/- 0.01672896 (3.88%) (init = 0.5)
    A0:      0.50059193 +/- 0.00389563 (0.78%) (init = 0.33)
[[Correlations]] (unreported correlations are < 0.100)
    C(scale, hwhm) = 0.790
    C(scale, A0)   = -0.655
    C(hwhm, A0)    = -0.499

Other option: plot fitting result and reference data using matplotlib.pyplot

fig1 = plt.figure()
plt.plot(xx, lorentzian_noisy, 'b-', label='reference data')
plt.plot(xx, result.best_fit, 'r.', label='fitting result')
plt.legend()
plt.xlabel('x')
plt.title('Fit result and reference data')
plt.grid();

Print values and errors of refined parameters:

for item in ['hwhm', 'center', 'scale']:
    print(f"{item}: {result.params[item].value} +/- {result.params[item].stderr} {dict_physical_units[item]}")

hwhm: 0.4315091873010395 +/- 0.016728957918933012 1/ps
center: -0.02611500946307296 +/- 0.010239707149920764 1/ps
scale: 0.9178578224996412 +/- 0.028838688012300575 unit_of_signal/ps

---

Generated by nbsphinx from a Jupyter notebook.