Curve fitting in ERLabPy largely relies on lmfit.Along with some convenient models for common fitting tasks, ERLabPy provides a powerfulaccessor that streamlines curve fitting on multidimensional xarray objects.
ERLabPy also provides optional integration of lmfit models with iminuit, which is a Python interface to the MinuitC++ library developed at CERN.
In this tutorial, we will start with the basics of curve fitting using lmfit, introducesome models that are available in ERLabPy, and demonstrate curve fitting with themodelfit accessor tofit multidimensional xarray objects. Finally, we will show how to use iminuit with lmfit models.
If you are familiar with lmfit, you can skip the first section.
[1]:
import erlab.plotting.erplot as epltimport matplotlib.pyplot as pltimport numpy as npimport xarray as xrfrom erlab.io.exampledata import generate_data
Let’s start by defining a model function and the data to fit.
[3]:
def poly1(x, a, b): return a * x + b# Generate some toy datax = np.linspace(0, 10, 20)y = poly1(x, 1, 2)# Add some noise with fixed seed for reproducibilityrng = np.random.default_rng(1)yerr = np.full_like(x, 0.5)y = rng.normal(y, yerr)
Fitting a simple function¶
A lmfit model can be created by calling lmfit.Model class with the modelfunction and the independent variable(s) as arguments.
[4]:
import lmfitmodel = lmfit.Model(poly1)params = model.make_params(a=1.0, b=2.0)result = model.fit(y, x=x, params=params, weights=1 / yerr)result.plot()result
[4]:
Fit Result
Model: Model(poly1)
| fitting method | leastsq |
| # function evals | 7 |
| # data points | 20 |
| # variables | 2 |
| chi-square | 5.86642273 |
| reduced chi-square | 0.32591237 |
| Akaike info crit. | -20.5297448 |
| Bayesian info crit. | -18.5382803 |
| R-squared | 0.99155992 |
| name | value | standard error | relative error | initial value | min | max | vary |
|---|---|---|---|---|---|---|---|
| a | 0.96713173 | 0.02103116 | (2.17%) | 1.0 | -inf | inf | True |
| b | 2.18308498 | 0.12301076 | (5.63%) | 2.0 | -inf | inf | True |
| Parameter1 | Parameter 2 | Correlation |
|---|---|---|
| a | b | -0.8549 |
By passing dictionaries to make_params, we can set the initial values of the parameters and also set the bounds for the parameters.
[5]:
model = lmfit.Model(poly1)params = model.make_params( a={"value": 1.0, "min": 0.0}, b={"value": 2.0, "vary": False},)result = model.fit(y, x=x, params=params, weights=1 / yerr)_ = result.plot()
result is a lmfit.model.ModelResult object that contains the results of thefit. The best-fit parameters can be accessed through the result.params attribute.
Note
Since all weights are the same in this case, it has little effect on the fitresults. However, if we are confident that we have a good estimate of yerr, wecan pass scale_covar=True to the fit method to obtain accurateuncertainties.
[6]:
result.params
[6]:
| name | value | standard error | relative error | initial value | min | max | vary |
|---|---|---|---|---|---|---|---|
| a | 0.99389030 | 0.01125608 | (1.13%) | 1.0 | 0.00000000 | inf | True |
| b | 2.00000000 | 0.00000000 | (0.00%) | 2.0 | -inf | inf | False |
[7]:
result.params["a"].value, result.params["a"].stderr
[7]:
(0.9938903037183116, 0.011256084507121714)
The parameters can also be retrieved in a form that allows easy error propagation calculation, enabled by the uncertainties package.
[8]:
a_uvar = result.uvars["a"]print(a_uvar)print(a_uvar**2)
0.994+/-0.0110.988+/-0.022
Fitting with composite models¶
Before fitting, let us generate a Gaussian peak on a linear background.
[9]:
# Generate toy datax = np.linspace(0, 10, 50)y = -0.1 * x + 2 + 3 * np.exp(-((x - 5) ** 2) / (2 * 1**2))# Add some noise with fixed seed for reproducibilityrng = np.random.default_rng(5)yerr = np.full_like(x, 0.3)y = rng.normal(y, yerr)# Plot the dataplt.errorbar(x, y, yerr, fmt="o")
[9]:
<ErrorbarContainer object of 3 artists>
A composite model can be created by adding multiple models together.
[10]:
from lmfit.models import GaussianModel, LinearModelmodel = GaussianModel() + LinearModel()params = model.make_params(slope=-0.1, center=5.0, sigma={"value": 0.1, "min": 0})params
[10]:
| name | value | initial value | min | max | vary | expression |
|---|---|---|---|---|---|---|
| amplitude | 1.00000000 | None | -inf | inf | True | |
| center | 5.00000000 | 5.0 | -inf | inf | True | |
| sigma | 0.10000000 | 0.1 | 0.00000000 | inf | True | |
| slope | -0.10000000 | -0.1 | -inf | inf | True | |
| intercept | 0.00000000 | None | -inf | inf | True | |
| fwhm | 0.23548200 | None | -inf | inf | False | 2.3548200*sigma |
| height | 3.98942300 | None | -inf | inf | False | 0.3989423*amplitude/max(1e-15, sigma) |
[11]:
result = model.fit(y, x=x, params=params, weights=1 / yerr)result.plot()result
[11]:
Fit Result
Model: (Model(gaussian) + Model(linear))
| fitting method | leastsq |
| # function evals | 43 |
| # data points | 50 |
| # variables | 5 |
| chi-square | 35.1640791 |
| reduced chi-square | 0.78142398 |
| Akaike info crit. | -7.59989616 |
| Bayesian info crit. | 1.96021887 |
| R-squared | 0.94677517 |
| name | value | standard error | relative error | initial value | min | max | vary | expression |
|---|---|---|---|---|---|---|---|---|
| amplitude | 7.16590435 | 0.37616176 | (5.25%) | 1.0 | -inf | inf | True | |
| center | 4.99152727 | 0.04220625 | (0.85%) | 5.0 | -inf | inf | True | |
| sigma | 0.94300480 | 0.04687009 | (4.97%) | 0.1 | 0.00000000 | inf | True | |
| slope | -0.11380431 | 0.01318512 | (11.59%) | -0.1 | -inf | inf | True | |
| intercept | 2.00633810 | 0.08441486 | (4.21%) | 0.0 | -inf | inf | True | |
| fwhm | 2.22060656 | 0.11037063 | (4.97%) | 0.23548200000000002 | -inf | inf | False | 2.3548200*sigma |
| height | 3.03156714 | 0.11942779 | (3.94%) | 3.989423 | -inf | inf | False | 0.3989423*amplitude/max(1e-15, sigma) |
| Parameter1 | Parameter 2 | Correlation |
|---|---|---|
| slope | intercept | -0.7822 |
| amplitude | sigma | +0.7041 |
| amplitude | intercept | -0.4390 |
| sigma | intercept | -0.3091 |
| center | slope | -0.2592 |
| center | intercept | +0.2027 |
How about multiple gaussian peaks? Since the parameter names overlap between the models, we must use the prefix argument to distinguish between them.
[12]:
model = GaussianModel(prefix="p0_") + GaussianModel(prefix="p1_") + LinearModel()model.make_params()
[12]:
| name | value | initial value | min | max | vary | expression |
|---|---|---|---|---|---|---|
| p0_amplitude | 1.00000000 | None | -inf | inf | True | |
| p0_center | 0.00000000 | None | -inf | inf | True | |
| p0_sigma | 1.00000000 | None | 0.00000000 | inf | True | |
| p1_amplitude | 1.00000000 | None | -inf | inf | True | |
| p1_center | 0.00000000 | None | -inf | inf | True | |
| p1_sigma | 1.00000000 | None | 0.00000000 | inf | True | |
| slope | 1.00000000 | None | -inf | inf | True | |
| intercept | 0.00000000 | None | -inf | inf | True | |
| p0_fwhm | 2.35482000 | None | -inf | inf | False | 2.3548200*p0_sigma |
| p0_height | 0.39894230 | None | -inf | inf | False | 0.3989423*p0_amplitude/max(1e-15, p0_sigma) |
| p1_fwhm | 2.35482000 | None | -inf | inf | False | 2.3548200*p1_sigma |
| p1_height | 0.39894230 | None | -inf | inf | False | 0.3989423*p1_amplitude/max(1e-15, p1_sigma) |
For more information, see the lmfit documentation.
Fitting with pre-defined models¶
Creating composite models with different prefixes every time can be cumbersome, soERLabPy provides some pre-defined models in erlab.analysis.fit.models.
Fitting multiple peaks¶
One example is MultiPeakModel, whichworks like a composite model of multiple Gaussian or Lorentzian peaks.
By supplying keyword arguments, you can specify the number of peaks, their shapes,whether to multiply with a Fermi-Dirac distribution, the shape of the background, andwhether to convolve the result with experimental resolution.
For a detailed explanation of all the arguments, see its documentation<erlab.analysis.fit.models.MultiPeakModel>.
The model can be constructed as follows:
[13]:
from erlab.analysis.fit.models import MultiPeakModelmodel = MultiPeakModel(npeaks=1, peak_shapes=["gaussian"], fd=False, convolve=False)params = model.make_params(p0_center=5.0, p0_width=0.2, p0_height=3.0)params
[13]:
| name | value | initial value | min | max | vary | expression |
|---|---|---|---|---|---|---|
| p0_center | 5.00000000 | 5.0 | -inf | inf | True | |
| p0_width | 0.20000000 | 0.2 | 0.00000000 | inf | True | |
| p0_height | 3.00000000 | 3.0 | 0.00000000 | inf | True | |
| const_bkg | 0.00000000 | None | -inf | inf | True | |
| lin_bkg | 0.00000000 | None | -inf | inf | True | |
| p0_sigma | 0.08493218 | None | -inf | inf | False | p0_width / (2 * sqrt(2 * log(2))) |
| p0_amplitude | 0.10164911 | None | -inf | inf | False | p0_height * p0_sigma / sqrt(2 * pi) |
[14]:
result = model.fit(y, x=x, params=params, weights=1 / yerr)_ = result.plot()
We can also plot components.
[15]:
comps = result.eval_components()plt.errorbar(x, y, yerr, fmt="o", zorder=-1, alpha=0.3)plt.plot(x, result.eval(), label="Best fit")plt.plot(x, comps["1Peak_p0"], "--", label="Peak")plt.plot(x, comps["1Peak_bkg"], "--", label="Background")plt.legend()
[15]:
<matplotlib.legend.Legend at 0x7f584c481110>
Now, let us try fitting MDCs cut from simulated data with multiple Lorentzian peaks, convolved with a common instrumental resolution.
[16]:
data = generate_data(bandshift=-0.2, count=5e8, seed=1).Tcut = data.qsel(ky=0.3)cut.qplot(colorbar=True)
[16]:
<matplotlib.image.AxesImage at 0x7f584c480610>
[17]:
mdc = cut.qsel(eV=0.0)mdc.qplot()
[17]:
[<matplotlib.lines.Line2D at 0x7f584c1ca590>]
First, we define the model and set the initial parameters.
[18]:
from erlab.analysis.fit.models import MultiPeakModelmodel = MultiPeakModel(npeaks=4, peak_shapes=["lorentzian"], fd=False, convolve=True)params = model.make_params( p0_center=-0.6, p1_center=-0.45, p2_center=0.45, p3_center=0.6, p0_width=0.02, p1_width=0.02, p2_width=0.02, p3_width=0.02, lin_bkg={"value": 0.0, "vary": False}, const_bkg=0.0, resolution=0.03,)params
[18]:
| name | value | initial value | min | max | vary | expression |
|---|---|---|---|---|---|---|
| p0_center | -0.60000000 | -0.6 | -inf | inf | True | |
| p0_width | 0.02000000 | 0.02 | 0.00000000 | inf | True | |
| p0_height | -inf | None | 0.00000000 | inf | True | |
| p1_center | -0.45000000 | -0.45 | -inf | inf | True | |
| p1_width | 0.02000000 | 0.02 | 0.00000000 | inf | True | |
| p1_height | -inf | None | 0.00000000 | inf | True | |
| p2_center | 0.45000000 | 0.45 | -inf | inf | True | |
| p2_width | 0.02000000 | 0.02 | 0.00000000 | inf | True | |
| p2_height | -inf | None | 0.00000000 | inf | True | |
| p3_center | 0.60000000 | 0.6 | -inf | inf | True | |
| p3_width | 0.02000000 | 0.02 | 0.00000000 | inf | True | |
| p3_height | -inf | None | 0.00000000 | inf | True | |
| const_bkg | 0.00000000 | 0.0 | -inf | inf | True | |
| lin_bkg | 0.00000000 | 0.0 | -inf | inf | False | |
| resolution | 0.03000000 | 0.03 | 0.00000000 | inf | True | |
| p0_sigma | 0.01000000 | None | -inf | inf | False | p0_width / 2 |
| p0_amplitude | -inf | None | -inf | inf | False | p0_height * p0_sigma * pi |
| p1_sigma | 0.01000000 | None | -inf | inf | False | p1_width / 2 |
| p1_amplitude | -inf | None | -inf | inf | False | p1_height * p1_sigma * pi |
| p2_sigma | 0.01000000 | None | -inf | inf | False | p2_width / 2 |
| p2_amplitude | -inf | None | -inf | inf | False | p2_height * p2_sigma * pi |
| p3_sigma | 0.01000000 | None | -inf | inf | False | p3_width / 2 |
| p3_amplitude | -inf | None | -inf | inf | False | p3_height * p3_sigma * pi |
Then, we can fit the model to the data:
[19]:
result = model.fit(mdc, x=mdc.kx, params=params)result.plot()result
[19]:
Fit Result
Model: Model(4Peak)
| fitting method | leastsq |
| # function evals | 305 |
| # data points | 250 |
| # variables | 14 |
| chi-square | 566.514610 |
| reduced chi-square | 2.40048564 |
| Akaike info crit. | 232.510488 |
| Bayesian info crit. | 281.810941 |
| R-squared | 0.99988540 |
| name | value | standard error | relative error | initial value | min | max | vary | expression |
|---|---|---|---|---|---|---|---|---|
| p0_center | -0.59650469 | 2.5399e-05 | (0.00%) | -0.6 | -inf | inf | True | |
| p0_width | 0.02403002 | 1.6855e-04 | (0.70%) | 0.02 | 0.00000000 | inf | True | |
| p0_height | 1161.03463 | 6.49245237 | (0.56%) | 0.0 | 0.00000000 | inf | True | |
| p1_center | -0.44489439 | 4.6412e-04 | (0.10%) | -0.45 | -inf | inf | True | |
| p1_width | 0.01804181 | 0.00148733 | (8.24%) | 0.02 | 0.00000000 | inf | True | |
| p1_height | 68.4523785 | 4.18812978 | (6.12%) | 0.0 | 0.00000000 | inf | True | |
| p2_center | 0.44353079 | 5.1714e-04 | (0.12%) | 0.45 | -inf | inf | True | |
| p2_width | 0.02372127 | 0.00164247 | (6.92%) | 0.02 | 0.00000000 | inf | True | |
| p2_height | 57.0667112 | 2.80977111 | (4.92%) | 0.0 | 0.00000000 | inf | True | |
| p3_center | 0.59606494 | 2.5489e-05 | (0.00%) | 0.6 | -inf | inf | True | |
| p3_width | 0.02443832 | 1.6821e-04 | (0.69%) | 0.02 | 0.00000000 | inf | True | |
| p3_height | 1152.96218 | 6.30752762 | (0.55%) | 0.0 | 0.00000000 | inf | True | |
| const_bkg | 0.27454009 | 0.13298720 | (48.44%) | 0.0 | -inf | inf | True | |
| lin_bkg | 0.00000000 | 0.00000000 | 0.0 | -inf | inf | False | ||
| resolution | 0.02634117 | 1.6289e-04 | (0.62%) | 0.03 | 0.00000000 | inf | True | |
| p0_sigma | 0.01201501 | 8.4274e-05 | (0.70%) | 0.01 | -inf | inf | False | p0_width / 2 |
| p0_amplitude | 43.8247214 | 0.08237775 | (0.19%) | 0.0 | -inf | inf | False | p0_height * p0_sigma * pi |
| p1_sigma | 0.00902091 | 7.4366e-04 | (8.24%) | 0.01 | -inf | inf | False | p1_width / 2 |
| p1_amplitude | 1.93994144 | 0.06140643 | (3.17%) | 0.0 | -inf | inf | False | p1_height * p1_sigma * pi |
| p2_sigma | 0.01186064 | 8.2124e-04 | (6.92%) | 0.01 | -inf | inf | False | p2_width / 2 |
| p2_amplitude | 2.12637915 | 0.06652226 | (3.13%) | 0.0 | -inf | inf | False | p2_height * p2_sigma * pi |
| p3_sigma | 0.01221916 | 8.4107e-05 | (0.69%) | 0.01 | -inf | inf | False | p3_width / 2 |
| p3_amplitude | 44.2594703 | 0.08273375 | (0.19%) | 0.0 | -inf | inf | False | p3_height * p3_sigma * pi |
| Parameter1 | Parameter 2 | Correlation |
|---|---|---|
| p0_width | p0_height | -0.9807 |
| p3_width | p3_height | -0.9801 |
| p1_width | p1_height | -0.9455 |
| p2_width | p2_height | -0.9151 |
| p0_height | resolution | +0.9068 |
| p3_height | resolution | +0.9054 |
| p0_width | resolution | -0.8802 |
| p3_width | resolution | -0.8785 |
| p0_height | p3_height | +0.8238 |
| p0_width | p3_height | -0.8029 |
| p0_height | p3_width | -0.8028 |
| p0_width | p3_width | +0.7866 |
| p0_width | const_bkg | -0.4211 |
| p3_width | const_bkg | -0.4178 |
| p0_height | const_bkg | +0.3711 |
| p3_height | const_bkg | +0.3683 |
| const_bkg | resolution | +0.3524 |
| p2_width | const_bkg | -0.2401 |
| p1_width | const_bkg | -0.2398 |
| p1_height | const_bkg | +0.1708 |
| p2_height | const_bkg | +0.1437 |
Fitting xarray objects¶
ERLabPy provides accessors for xarray objects that allows you to fit data with lmfitmodels: xarray.DataArray.modelfit orxarray.Dataset.modelfit,depending on whether you want to fit a single DataArray or multiple DataArrays in aDataset. The syntax is similar to xarray.DataArray.curvefit() andxarray.Dataset.curvefit(). The accessor returns a xarray.Dataset with thebest-fit parameters and the fit statistics. The example below illustrates how to use theaccessor to conduct the same fit as in the previous example.
[20]:
result_ds = mdc.modelfit("kx", model, params=params)result_ds
[20]:
<xarray.Dataset> Size: 14kBDimensions: (param: 23, cov_i: 23, cov_j: 23, fit_stat: 8, kx: 250)Coordinates: ky float64 8B 0.2967 eV float64 8B -0.000301 * kx (kx) float64 2kB -0.89 -0.8829 ... 0.8829 0.89 * param (param) <U12 1kB 'p0_center' ... 'p3_amplitude' * fit_stat (fit_stat) <U6 192B 'nfev' 'nvarys' ... 'aic' 'bic' * cov_i (cov_i) <U12 1kB 'p0_center' ... 'p3_amplitude' * cov_j (cov_j) <U12 1kB 'p0_center' ... 'p3_amplitude'Data variables: modelfit_results object 8B <lmfit.model.ModelResult object at 0x7f5... modelfit_coefficients (param) float64 184B -0.5965 0.02403 ... 44.26 modelfit_stderr (param) float64 184B 2.54e-05 0.0001685 ... 0.08273 modelfit_covariance (cov_i, cov_j) float64 4kB 6.451e-10 ... nan modelfit_stats (fit_stat) float64 64B 305.0 14.0 ... 232.5 281.8 modelfit_data (kx) float64 2kB 2.337 2.061 2.908 ... 3.202 4.133 modelfit_best_fit (kx) float64 2kB 2.336 2.437 2.545 ... 2.495 2.392
modelfit_results contains the underlying lmfit.model.ModelResult objects.The coefficients and their errors are stored in modelfit_coefficients andmodelfit_stderr, respectively. The goodness-of-fit statistics are storedin modelfit_stats.
Fitting across multiple dimensions¶
Note
There is a dedicated module for Fermi edge fitting and correction, describedhere. The following example is for illustrative purposes.
The accessor comes in handy when you have to fit a single model to multiple data pointsacross arbitrary dimensions. Let’s demonstrate this with a simulated cut that representsa curved Fermi edge at 100 K, with an energy broadening of 20 meV.
[21]:
from erlab.io.exampledata import generate_gold_edgefrom erlab.analysis.fit.models import FermiEdgeModelgold = generate_gold_edge(temp=100, Eres=0.02, count=5e5, seed=1)gold.qplot(cmap="Greys")
[21]:
<matplotlib.image.AxesImage at 0x7f5847397e50>
We first select ± 0.2 eV around the Fermi level and fit the model across the energy axis for every EDC.
[22]:
gold_selected = gold.sel(eV=slice(-0.2, 0.2))result_ds = gold_selected.modelfit( coords="eV", model=FermiEdgeModel(), params={"temp": {"value": 100.0, "vary": False}}, guess=True,)result_ds
[22]:
<xarray.Dataset> Size: 358kBDimensions: (alpha: 200, param: 7, cov_i: 7, cov_j: 7, fit_stat: 8, eV: 75)Coordinates: * alpha (alpha) float64 2kB -15.0 -14.85 -14.7 ... 14.85 15.0 * eV (eV) float64 600B -0.1977 -0.1923 ... 0.193 0.1983 * param (param) <U10 280B 'center' 'temp' ... 'dos0' 'dos1' * fit_stat (fit_stat) <U6 192B 'nfev' 'nvarys' ... 'aic' 'bic' * cov_i (cov_i) <U10 280B 'center' 'temp' ... 'dos0' 'dos1' * cov_j (cov_j) <U10 280B 'center' 'temp' ... 'dos0' 'dos1'Data variables: modelfit_results (alpha) object 2kB <lmfit.model.ModelResult object... modelfit_coefficients (alpha, param) float64 11kB -0.01774 100.0 ... -3.53 modelfit_stderr (alpha, param) float64 11kB 0.005598 0.0 ... 4.111 modelfit_covariance (alpha, cov_i, cov_j) float64 78kB 3.134e-05 ... 16.9 modelfit_stats (alpha, fit_stat) float64 13kB 106.0 6.0 ... -26.34 modelfit_data (alpha, eV) float64 120kB 5.1 7.021 7.453 ... 0.0 0.0 modelfit_best_fit (alpha, eV) float64 120kB 6.642 6.565 ... -0.0711
Let’s plot the fitted parameters as a function of angle!
[23]:
gold.qplot(cmap="Greys")plt.errorbar( gold_selected.alpha, result_ds.modelfit_coefficients.sel(param="center"), result_ds.modelfit_stderr.sel(param="center"), fmt=".",)
[23]:
<ErrorbarContainer object of 3 artists>
We can also conduct a global fit to all the EDCs with a multidimensional model. First,we normalize the data with the averaged intensity of each EDC and then fit the data toFermiEdge2dModel.
[24]:
from erlab.analysis.fit.models import FermiEdge2dModelgold_norm = gold_selected / gold_selected.mean("eV")result_ds = gold_norm.T.modelfit( coords=["eV", "alpha"], model=FermiEdge2dModel(), params={"temp": {"value": 100.0, "vary": False}}, guess=True,)result_ds
[24]:
<xarray.Dataset> Size: 244kBDimensions: (param: 8, cov_i: 8, cov_j: 8, fit_stat: 8, eV: 75, alpha: 200)Coordinates: * eV (eV) float64 600B -0.1977 -0.1923 ... 0.193 0.1983 * alpha (alpha) float64 2kB -15.0 -14.85 -14.7 ... 14.85 15.0 * param (param) <U10 320B 'c0' 'c1' ... 'offset' 'resolution' * fit_stat (fit_stat) <U6 192B 'nfev' 'nvarys' ... 'aic' 'bic' * cov_i (cov_i) <U10 320B 'c0' 'c1' ... 'offset' 'resolution' * cov_j (cov_j) <U10 320B 'c0' 'c1' ... 'offset' 'resolution'Data variables: modelfit_results object 8B <lmfit.model.ModelResult object at 0x7f5... modelfit_coefficients (param) float64 64B 0.03604 8.09e-06 ... 0.04345 modelfit_stderr (param) float64 64B 0.0004463 2.749e-05 ... 0.001243 modelfit_covariance (cov_i, cov_j) float64 512B 1.992e-07 ... 1.545e-06 modelfit_stats (fit_stat) float64 64B 82.0 7.0 ... -3.98e+04 modelfit_data (eV, alpha) float64 120kB 2.027 1.905 ... 0.0 0.0 modelfit_best_fit (eV, alpha) float64 120kB 1.965 1.965 ... 0.02488
Let’s plot the fit results and the residuals.
[25]:
best_fit = result_ds.modelfit_best_fit.transpose(*gold_norm.dims)fig, axs = eplt.plot_slices( [gold_norm, best_fit, best_fit - gold_norm], figsize=(4, 5), cmap=["viridis", "viridis", "bwr"], norm=[plt.Normalize(), plt.Normalize(), eplt.CenteredPowerNorm(1.0, vcenter=0)], colorbar="all", hide_colorbar_ticks=False, colorbar_kw={"width": 7},)eplt.set_titles(axs, ["Data", "FermiEdge2dModel", "Residuals"])
Providing multiple initial guesses¶
You can also provide different initial guesses and bounds for different coordinates. Todemonstrate, let’s create some data containing multiple MDCs, each with a different peakposition.
[26]:
# Define angle coordinates for 2D dataalpha = np.linspace(-5.0, 5.0, 100)beta = np.linspace(-1.0, 1.0, 3)# Center of the peaks along betacenter = np.array([-2.0, 0.0, 2.0])[:, np.newaxis]# Gaussian peak on a linear backgroundy = -0.1 * alpha + 2 + 3 * np.exp(-((alpha - center) ** 2) / (2 * 1**2))# Add some noise with fixed seed for reproducibilityrng = np.random.default_rng(5)yerr = np.full_like(y, 0.1)y = rng.normal(y, yerr)# Transform to DataArraydarr = xr.DataArray(y, dims=["beta", "alpha"], coords={"beta": beta, "alpha": alpha})darr.qplot()
[26]:
<matplotlib.image.AxesImage at 0x7f58472a80d0>
We can provide different initial guesses for the peak positions along the betadimension by passing a dictionary of DataArrays to the params argument.
[27]:
result_ds = darr.modelfit( coords="alpha", model=GaussianModel() + LinearModel(), params={ "center": xr.DataArray([-2, 0, 2], coords=[darr.beta]), "slope": -0.1, },)result_ds
[27]:
<xarray.Dataset> Size: 7kBDimensions: (beta: 3, param: 5, cov_i: 5, cov_j: 5, fit_stat: 8, alpha: 100)Coordinates: * beta (beta) float64 24B -1.0 0.0 1.0 * alpha (alpha) float64 800B -5.0 -4.899 -4.798 ... 4.899 5.0 * param (param) <U9 180B 'amplitude' 'center' ... 'intercept' * fit_stat (fit_stat) <U6 192B 'nfev' 'nvarys' ... 'aic' 'bic' * cov_i (cov_i) <U9 180B 'amplitude' 'center' ... 'intercept' * cov_j (cov_j) <U9 180B 'amplitude' 'center' ... 'intercept'Data variables: modelfit_results (beta) object 24B <lmfit.model.ModelResult object ... modelfit_coefficients (beta, param) float64 120B 7.324 -2.0 ... 1.993 modelfit_stderr (beta, param) float64 120B 0.1349 0.011 ... 0.01735 modelfit_covariance (beta, cov_i, cov_j) float64 600B 0.01819 ... 0.00... modelfit_stats (beta, fit_stat) float64 192B 31.0 5.0 ... -450.2 modelfit_data (beta, alpha) float64 2kB 2.453 2.402 ... 1.404 1.64 modelfit_best_fit (beta, alpha) float64 2kB 2.553 2.553 ... 1.526 1.504
Let’s overlay the fitted peak positions on the data.
[28]:
result_ds.modelfit_data.qplot()result_center = result_ds.sel(param="center")plt.plot(result_center.modelfit_coefficients, result_center.beta, '.-')
[28]:
[<matplotlib.lines.Line2D at 0x7f584c64cf50>]
The same can be done with all parameter attributes that can be passed tolmfit.create_params() (e.g., vary, min, max, etc.). For example:
[29]:
result_ds = darr.modelfit( coords="alpha", model=GaussianModel() + LinearModel(), params={ "center": { "value": xr.DataArray([-2, 0, 2], coords=[darr.beta]), "min": -5.0, "max": xr.DataArray([0, 2, 5], coords=[darr.beta]), }, "slope": -0.1, },)result_ds
[29]:
<xarray.Dataset> Size: 7kBDimensions: (beta: 3, param: 5, cov_i: 5, cov_j: 5, fit_stat: 8, alpha: 100)Coordinates: * beta (beta) float64 24B -1.0 0.0 1.0 * alpha (alpha) float64 800B -5.0 -4.899 -4.798 ... 4.899 5.0 * param (param) <U9 180B 'amplitude' 'center' ... 'intercept' * fit_stat (fit_stat) <U6 192B 'nfev' 'nvarys' ... 'aic' 'bic' * cov_i (cov_i) <U9 180B 'amplitude' 'center' ... 'intercept' * cov_j (cov_j) <U9 180B 'amplitude' 'center' ... 'intercept'Data variables: modelfit_results (beta) object 24B <lmfit.model.ModelResult object ... modelfit_coefficients (beta, param) float64 120B 7.324 -2.0 ... 1.993 modelfit_stderr (beta, param) float64 120B 0.1349 0.011 ... 0.01735 modelfit_covariance (beta, cov_i, cov_j) float64 600B 0.01819 ... 0.00... modelfit_stats (beta, fit_stat) float64 192B 31.0 5.0 ... -450.2 modelfit_data (beta, alpha) float64 2kB 2.453 2.402 ... 1.404 1.64 modelfit_best_fit (beta, alpha) float64 2kB 2.553 2.553 ... 1.526 1.504
Visualizing fits¶
If hvplot is installed, we can visualize the fitresults interactively with the qshow accessor.
Note
If you are viewing this documentation online, the plot will not be interactive. Run the code locally to try it out.
[30]:
result_ds.qshow(plot_components=True)
[30]:
Parallelization¶
The accessors are tightly integrated with xarray, so passing a dask array willparallelize the fitting process. See Parallel Computing with Dask for more information.
For non-dask objects, you can achieve joblib-based parallelization:
For non-dask Datasets, basic parallelization across multiple data variables can beachieved with the
parallelargument to xarray.Dataset.modelfit.For parallelizing fits on a single DataArray along a dimension with many points, thexarray.DataArray.parallel_fit accessor can be used. This methodis similar to xarray.DataArray.modelfit, but requires the name ofthe dimension to parallelize over instead of the coordinates to fit along. Forexample, to parallelize the fit in the previous example, you can use the followingcode:
gold_selected.parallel_fit( dim="alpha", model=FermiEdgeModel(), params={"temp": {"value": 100.0, "vary": False}}, guess=True,)
Note
Note that the initial run will take a long time due to the overhead of creatingparallel workers. Subsequent calls will run faster, since joblib’s default backendwill try to reuse the workers.
The accessor has some intrinsic overhead due to post-processing. If you need thebest performance, handle the parallelization yourself with joblib and
lmfit.Model.fit.
Saving and loading fits¶
Since the resulting dataset contains no Python objects, it can be saved and loadedreliably as a netCDF file using xarray.Dataset.to_netcdf() andxarray.open_dataset(), respectively. If you wish to obtain thelmfit.model.ModelResult objects from the fit, you can use the output_resultargument.
Warning
Saving full model results that includes the model functions can be difficult.Instead of saving the fit results, it is recommended to save the code that canreproduce the fit. See the relevant lmfit documentation formore information.
Fermi edge fitting¶
Functions related to the Fermi edge are available in erlab.analysis.gold. To fita polynomial to a Fermi edge, you can use erlab.analysis.gold.poly().
[32]:
import erlab.analysis as eraimport erlab.plotting.erplot as epltgold = generate_gold_edge(temp=100, seed=1)result = era.gold.poly( gold, angle_range=(-15, 15), eV_range=(-0.2, 0.2), temp=100.0, vary_temp=False, degree=2, plot=True,)
The resulting polynomial can be used to correct the Fermi edge witherlab.analysis.gold.correct_with_edge().
[33]:
era.correct_with_edge(gold, result).qplot(cmap="Greys")eplt.fermiline()
[33]:
<matplotlib.lines.Line2D at 0x7f583f08a2d0>
Also check out the interactive Fermi edge fitting tool,erlab.interactive.goldtool().
Using iminuit¶
Note
This part requires iminuit.
iminuit is a powerful Python interface to theMinuit C++ library developed atCERN. To learn more, see the iminuit documentation.
ERLabPy provides a thin wrapper around iminuit.Minuit that allows you to uselmfit models with iminuit. The example below conducts the same fit as the previous one,but using iminuit.
[34]:
from erlab.analysis.fit.minuit import Minuitfrom erlab.analysis.fit.models import MultiPeakModelmodel = MultiPeakModel(npeaks=4, peak_shapes=["lorentzian"], fd=False, convolve=True)m = Minuit.from_lmfit( model, mdc, mdc.kx, p0_center=-0.6, p1_center=-0.45, p2_center=0.45, p3_center=0.6, p0_width=0.02, p1_width=0.02, p2_width=0.02, p3_width=0.02, p0_height=1500, p1_height=50, p2_height=50, p3_height=1500, lin_bkg={"value": 0.0, "vary": False}, const_bkg=0.0, resolution=0.03,)m.migrad()m.minos()m.hesse()
[34]:
| Migrad | |
|---|---|
| FCN = 566.5 (χ²/ndof = 2.4) | Nfcn = 5120 |
| EDM = 4.13e-06 (Goal: 0.0002) | time = 0.9 sec |
| Valid Minimum | Below EDM threshold (goal x 10) |
| No parameters at limit | Below call limit |
| Hesse ok | Covariance accurate |
| Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
|---|---|---|---|---|---|---|---|---|
| 0 | p0_center | -596.505e-3 | 0.016e-3 | -0.016e-3 | 0.016e-3 | |||
| 1 | p0_width | 24.03e-3 | 0.11e-3 | -0.11e-3 | 0.11e-3 | 0 | ||
| 2 | p0_height | 1.161e3 | 0.004e3 | -0.004e3 | 0.004e3 | 0 | ||
| 3 | p1_center | -444.89e-3 | 0.30e-3 | -0.30e-3 | 0.30e-3 | |||
| 4 | p1_width | 18e-3 | 1e-3 | -1e-3 | 1e-3 | 0 | ||
| 5 | p1_height | 68.5 | 2.8 | -2.7 | 2.9 | 0 | ||
| 6 | p2_center | 443.53e-3 | 0.34e-3 | -0.34e-3 | 0.34e-3 | |||
| 7 | p2_width | 0.0237 | 0.0010 | -0.0010 | 0.0010 | 0 | ||
| 8 | p2_height | 57.1 | 1.8 | -1.7 | 1.8 | 0 | ||
| 9 | p3_center | 596.065e-3 | 0.016e-3 | -0.016e-3 | 0.017e-3 | |||
| 10 | p3_width | 24.44e-3 | 0.11e-3 | -0.11e-3 | 0.11e-3 | 0 | ||
| 11 | p3_height | 1.153e3 | 0.004e3 | -0.004e3 | 0.004e3 | 0 | ||
| 12 | const_bkg | 0.27 | 0.09 | -0.09 | 0.09 | |||
| 13 | lin_bkg | 0.0 | 0.1 | yes | ||||
| 14 | resolution | 26.34e-3 | 0.10e-3 | -0.10e-3 | 0.10e-3 | 0 |
| p0_center | p0_width | p0_height | p1_center | p1_width | p1_height | p2_center | p2_width | p2_height | p3_center | p3_width | p3_height | const_bkg | resolution | |||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Error | -0.016e-3 | 0.016e-3 | -0.11e-3 | 0.11e-3 | -4 | 4 | -0.3e-3 | 0.3e-3 | -1e-3 | 1e-3 | -2.7 | 2.9 | -0.34e-3 | 0.34e-3 | -0.001 | 0.001 | -1.7 | 1.8 | -0.016e-3 | 0.017e-3 | -0.11e-3 | 0.11e-3 | -4 | 4 | -0.09 | 0.09 | -0.1e-3 | 0.1e-3 |
| Valid | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True | True |
| At Limit | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False |
| Max FCN | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False |
| New Min | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False | False |
| p0_center | p0_width | p0_height | p1_center | p1_width | p1_height | p2_center | p2_width | p2_height | p3_center | p3_width | p3_height | const_bkg | lin_bkg | resolution | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| p0_center | 2.68e-10 | 0.01e-9 (0.004) | -309.64e-9 (-0.005) | 0.02e-9 (0.004) | -0.47e-9 (-0.029) | 1.09497e-6 (0.024) | -0 | -0.04e-9 (-0.002) | 24.04e-9 | -0 | 0.01e-9 (0.004) | -283.08e-9 (-0.004) | 9.90e-9 (0.007) | 0 | -0.01e-9 (-0.005) |
| p0_width | 0.01e-9 (0.004) | 1.17e-08 | -442.668e-6 (-0.981) | 0.001e-6 (0.036) | 0.001e-6 (0.007) | -10.032e-6 (-0.033) | -0.002e-6 (-0.045) | 0.006e-6 (0.054) | -13.035e-6 (-0.068) | 0.01e-9 (0.008) | 0.009e-6 (0.785) | -352.815e-6 (-0.801) | -3.894e-6 (-0.419) | 0 | -0.010e-6 (-0.879) |
| p0_height | -309.64e-9 (-0.005) | -442.668e-6 (-0.981) | 17.4 | -41.28e-6 (-0.033) | -63.7e-6 (-0.015) | 1 (0.043) | 66.39e-6 (0.047) | -173.7e-6 (-0.040) | 0.5 (0.062) | -523.22e-9 (-0.008) | -362.570e-6 (-0.802) | 14 (0.823) | 0.132 (0.369) | 0 | 394.794e-6 (0.906) |
| p1_center | 0.02e-9 (0.004) | 0.001e-6 (0.036) | -41.28e-6 (-0.033) | 8.75e-08 | 0 (0.016) | -7.69e-6 (-0.009) | -0 (-0.001) | 0 (0.004) | -1.73e-6 (-0.003) | 0 | 0.001e-6 (0.026) | -31.38e-6 (-0.026) | -0.59e-6 (-0.023) | 0 | -0.001e-6 (-0.028) |
| p1_width | -0.47e-9 (-0.029) | 0.001e-6 (0.007) | -63.7e-6 (-0.015) | 0 (0.016) | 9.86e-07 | -2.6249e-3 (-0.949) | -0 | 0.1e-6 (0.057) | -53.5e-6 (-0.030) | 0.04e-9 (0.003) | 0.007e-6 (0.067) | -212.3e-6 (-0.053) | -20.8e-6 (-0.244) | 0e-6 | -0.005e-6 (-0.045) |
| p1_height | 1.09497e-6 (0.024) | -10.032e-6 (-0.033) | 1 (0.043) | -7.69e-6 (-0.009) | -2.6249e-3 (-0.949) | 7.76 | 2.16e-6 (0.002) | -111.5e-6 (-0.039) | 0.1 (0.023) | -98.11e-9 (-0.002) | -24.143e-6 (-0.080) | 1 (0.072) | 0.042 (0.177) | 0 | 20.462e-6 (0.070) |
| p2_center | -0 | -0.002e-6 (-0.045) | 66.39e-6 (0.047) | -0 (-0.001) | -0 | 2.16e-6 (0.002) | 1.15e-07 | 0.03e-6 (0.084) | -44.92e-6 (-0.075) | 0.05e-9 (0.009) | -0.002e-6 (-0.065) | 83.33e-6 (0.060) | 0.30e-6 (0.010) | 0 | 0.002e-6 (0.052) |
| p2_width | -0.04e-9 (-0.002) | 0.006e-6 (0.054) | -173.7e-6 (-0.040) | 0 (0.004) | 0.1e-6 (0.057) | -111.5e-6 (-0.039) | 0.03e-6 (0.084) | 1.07e-06 | -1.6728e-3 (-0.911) | 0.57e-9 (0.033) | -0.001e-6 (-0.011) | -0.7e-6 | -20.8e-6 (-0.234) | 0 | -0.003e-6 (-0.031) |
| p2_height | 24.04e-9 | -13.035e-6 (-0.068) | 0.5 (0.062) | -1.73e-6 (-0.003) | -53.5e-6 (-0.030) | 0.1 (0.023) | -44.92e-6 (-0.075) | -1.6728e-3 (-0.911) | 3.15 | -752.74e-9 (-0.026) | -4.459e-6 (-0.023) | 0.3 (0.036) | 0.021 (0.136) | 0.0 | 11.425e-6 (0.062) |
| p3_center | -0 | 0.01e-9 (0.008) | -523.22e-9 (-0.008) | 0 | 0.04e-9 (0.003) | -98.11e-9 (-0.002) | 0.05e-9 (0.009) | 0.57e-9 (0.033) | -752.74e-9 (-0.026) | 2.71e-10 | 0.01e-9 (0.005) | -373.21e-9 (-0.006) | -16.96e-9 (-0.012) | 0 | -0.01e-9 (-0.008) |
| p3_width | 0.01e-9 (0.004) | 0.009e-6 (0.785) | -362.570e-6 (-0.802) | 0.001e-6 (0.026) | 0.007e-6 (0.067) | -24.143e-6 (-0.080) | -0.002e-6 (-0.065) | -0.001e-6 (-0.011) | -4.459e-6 (-0.023) | 0.01e-9 (0.005) | 1.18e-08 | -432.171e-6 (-0.980) | -3.878e-6 (-0.417) | 0 | -0.010e-6 (-0.879) |
| p3_height | -283.08e-9 (-0.004) | -352.815e-6 (-0.801) | 14 (0.823) | -31.38e-6 (-0.026) | -212.3e-6 (-0.053) | 1 (0.072) | 83.33e-6 (0.060) | -0.7e-6 | 0.3 (0.036) | -373.21e-9 (-0.006) | -432.171e-6 (-0.980) | 16.5 | 0.128 (0.367) | 0 | 384.873e-6 (0.905) |
| const_bkg | 9.90e-9 (0.007) | -3.894e-6 (-0.419) | 0.132 (0.369) | -0.59e-6 (-0.023) | -20.8e-6 (-0.244) | 0.042 (0.177) | 0.30e-6 (0.010) | -20.8e-6 (-0.234) | 0.021 (0.136) | -16.96e-9 (-0.012) | -3.878e-6 (-0.417) | 0.128 (0.367) | 0.00736 | 0.000 | 3.146e-6 (0.351) |
| lin_bkg | 0 | 0 | 0 | 0 | 0e-6 | 0 | 0 | 0 | 0.0 | 0 | 0 | 0 | 0.000 | 0 | 0 |
| resolution | -0.01e-9 (-0.005) | -0.010e-6 (-0.879) | 394.794e-6 (0.906) | -0.001e-6 (-0.028) | -0.005e-6 (-0.045) | 20.462e-6 (0.070) | 0.002e-6 (0.052) | -0.003e-6 (-0.031) | 11.425e-6 (0.062) | -0.01e-9 (-0.008) | -0.010e-6 (-0.879) | 384.873e-6 (0.905) | 3.146e-6 (0.351) | 0 | 1.09e-08 |
You can also use the interactive fitting interface provided by iminuit.
Note
Requires ipywidgets to beinstalled.
If you are viewing this documentation online, changing the sliders won’t changethe plot. run the code locally to try it out.
[35]:
m.interactive()
[35]:
