Gaussian Model 3D with lmfit
Introduction
The objective of this notebook is to show how to use the Gaussian Model 3D model to perform some fits using lmfit.
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.
[1]:
# Units of parameters for selected QENS model and experimental data
dict_physical_units = {'omega': "meV",
'q': "1/Angstrom",
'D': "meV.Angstrom^2",
'variance_ux': "Angstrom",
'scale': "unit_of_signal.meV",
'center': "meV"}
Import libraries
[2]:
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets
import lmfit
import QENSmodels
Plot fitting model
The widget below shows the peak shape function imported from QENSmodels where the function’s parameters can be varied.
[3]:
# Dictionary of initial values
ini_parameters = {'q': 1., 'scale': 5., 'center': 0., 'D': 1., 'variance_ux': 1.}
def interactive_fct(q, scale, center, D, variance_ux):
xs = np.linspace(-10, 10, 100)
fig1, ax1 = plt.subplots()
ax1.plot(xs, QENSmodels.sqwGaussianModel3D(xs, q, scale, center, D, variance_ux))
ax1.set_xlabel('x')
ax1.grid()
# Define sliders for modifiable parameters and their range of variations
q_slider = ipywidgets.FloatSlider(value=ini_parameters['q'],
min=0.1, max=10., step=0.1,
description='q',
continuous_update=False)
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)
D_slider = ipywidgets.FloatSlider(value=ini_parameters['D'],
min=0.1, max=10, step=0.1,
description='D',
continuous_update=False)
variance_ux_slider = ipywidgets.FloatSlider(value=ini_parameters['variance_ux'],
min=0.1, max=10, step=0.1,
description='variance_ux',
continuous_update=False)
grid_sliders = ipywidgets.HBox([ipywidgets.VBox([q_slider, scale_slider, center_slider]),
ipywidgets.VBox([D_slider, variance_ux_slider])])
# Define function to reset all parameters' values to the initial ones
def reset_values(b):
"""
Reset the interactive plots to inital values
"""
q_slider.value = ini_parameters['q']
scale_slider.value = ini_parameters['scale']
center_slider.value = ini_parameters['center']
D_slider.value = ini_parameters['D']
variance_ux_slider.value = ini_parameters['variance_ux']
# 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,
{'q': q_slider,
'scale': scale_slider,
'center': center_slider,
'D': D_slider,
'variance_ux': variance_ux_slider})
display(grid_sliders, interactive_plot, reset_button)
Creating the reference data
[4]:
nb_points = 200
xx = np.linspace(-2, 2, nb_points)
added_noise = np.random.normal(0, 1, nb_points)
gaussian_model_3d_noisy = QENSmodels.sqwGaussianModel3D(
xx,
q=2.4,
scale=1,
center=0,
D=0.1,
variance_ux=5.4
) * (1 + 0.1 * added_noise) + 0.01 * added_noise
# plot initial mode
fig0 = plt.figure()
plt.plot(xx, gaussian_model_3d_noisy)
plt.grid();
Setting and fitting
[5]:
model = lmfit.Model(QENSmodels.sqwGaussianModel3D)
print(f'Names of parameters: {model.param_names}\nIndependent variable(s): {model.independent_vars}')
initial_parameters_values = {'scale': 1.22, 'center': 0.01, 'D': 0.2, 'variance_ux': 4}
# Define boundaries for parameters to be refined
model.set_param_hint('scale', min=0)
model.set_param_hint('center', min=-5, max=5)
model.set_param_hint('D', min=0)
model.set_param_hint('variance_ux', min=0)
# Fix some of the parameters
model.set_param_hint('q', vary=False)
# Fit
result = model.fit(
gaussian_model_3d_noisy,
w=xx,
q=1.,
scale=initial_parameters_values['scale'],
center=initial_parameters_values['center'],
D=initial_parameters_values['D'],
variance_ux=initial_parameters_values['variance_ux']
)
Names of parameters: ['q', 'scale', 'center', 'D', 'variance_ux']
Independent variable(s): ['w']
[6]:
# Plot Initial model and reference data
fig1 = plt.figure()
plt.plot(xx, gaussian_model_3d_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();
Plotting results
using methods implemented in lmfit
[7]:
# display result
print('Result of fit:\n',result.fit_report())
# plot fitting results using lmfit functionality
result.plot();
Result of fit:
[[Model]]
Model(sqwGaussianModel3D)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 56
# data points = 200
# variables = 4
chi-square = 0.25275757
reduced chi-square = 0.00128958
Akaike info crit = -1326.72836
Bayesian info crit = -1313.53510
R-squared = 0.95803193
[[Variables]]
q: 1 (fixed)
scale: 0.99310908 +/- 0.01786048 (1.80%) (init = 1.22)
center: -0.00226924 +/- 0.00720230 (317.39%) (init = 0.01)
D: 0.55409556 +/- 0.02542888 (4.59%) (init = 0.2)
variance_ux: 37.9474121 +/- 53.5252613 (141.05%) (init = 4)
[[Correlations]] (unreported correlations are < 0.100)
C(D, variance_ux) = -0.901
C(scale, D) = 0.814
C(scale, variance_ux) = -0.630
Other option: plot fitting results and reference data using matplotlib.pyplot
[8]:
fig2 = plt.figure()
plt.plot(xx, gaussian_model_3d_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();
[ ]: