Construction and fitting of modular lineshape models

In this tutorial, a basic example is shown how qspec.models can be used to construct lineshape models for data fitting.

First, we import all the required modules and generate our lineshape model. We can print all parameters of the model and set them to new values. Additionally, we can fix certain parameters for fitting by changing the list model.fixes. Here, we define sigma as a Prior, which is basically a known parameter with uncertainty.

import numpy as np
import qspec.models as mod
import matplotlib.pyplot as plt

# Generate the model used in this tutorial
model = mod.Offset(  # Add a y-axis shift (off0e0)

    mod.NPeak(  # Add a single (n_peaks=1)
        # x-axis shift (x0) and an intensity (p0)

        mod.Voigt(), n_peaks=1)) # Add a Voigt lineshape
        # with Lorentzian (Gamma) and Gaussian (sigma) widths

print(f'\nParameter names: {model.names}\n')
# >>> Parameter names: ['Gamma', 'sigma', 'x0', 'p0', 'off0e0']

# Change all parameter values and if they stay fixed during fitting at once
model.set_vals([20., 6., 0., 10., 3.])
model.set_fixes([False, '5(0.7)', False, False, False])

Now, we generate some data for the fit. The data is normally distributed around the model with its current parameters.

# Generate some random data for this tutorial
x = np.linspace(-80., 80., 81)  # x-values
sigma_y = np.full_like(x, 0.5)  # y-uncertainties
y = np.random.normal(model(x, *model.vals), sigma_y)
# Random y-values around the model

We change the initial parameters of the peak position and intensity to later compare the model before and after fitting. During the fit, model.vals is set to the optimized parameters. Therefore, we save our initial values to p_init. We set report=True to get a summary of the fit results.

# Change the initial value of the peak position and intensity
model.set_val(2, 15.)  # Parameter 2 (x0), the peak position
model.set_val(3, 5.)  # Parameter 3 (p0), the peak intensity
p_init = model.vals.copy()  # Copy the initial values for the plot

# Fit the constructed model to the data
popt, pcov, info = mod.fit(
    model, x, y, sigma_y=sigma_y, report=True)

Optimized parameters:
0:   Gamma = 23.07256228925477 +/- 1.6786599666472717
1:   sigma = 4.697853592625128 +/- 0.7099638812366352
2:   x0 = 0.3717374185209716 +/- 0.29514512340842586
3:   p0 = 10.584453120309163 +/- 0.2375261208105466
4:   off0e0 = 2.799096283098853 +/- 0.11049924526028845

Cov. Matrix:
0:    1.00   -0.72   -0.00    0.03   -0.68
1:   -0.72    1.00    0.00   -0.34    0.32
2:   -0.00    0.00    1.00   -0.00    0.00
3:    0.03   -0.34   -0.00    1.00   -0.31
4:   -0.68    0.32    0.00   -0.31    1.00

Red. chi2 = 1.17
Fit successful.

Finally, we plot everything.

# Plot the data, the initial model and the fitted model
plt.errorbar(x, y, yerr=sigma_y, fmt='.k', label='Data')
plt.plot(x, model(x, *p_init), '-C0', label='Initial model')
plt.plot(x, model(x, *popt), '-C1', label='Fitted model')

# Improve the plot
plt.legend()
plt.xlabel('x'), plt.ylabel('y')
plt.subplots_adjust(left=0.08, bottom=0.09, top=0.99, right=0.99)
plt.show()

←  Previous

Next  →