Comparing Gaussian, Lorentzian, and Voigt 1D models

Overview

In this example we will try to fit a peaked profile with a range of 1D models. If you have read the lmfit documentation then this example should seem familiar!

Setting up

The following sections will load in classes from Sherpa as needed, but it is assumed that the following module has been loaded:

>>> import matplotlib.pyplot as plt

Loading the data

The data can be retrieved from the lmfit GitHub page. It is a two-column ASCII file that can be read in with NumPy or the sherpa.io.read_data() routine:

>>> from sherpa.io import read_data
>>> d = read_data('test_peak.dat')
>>> print(d)
name      = test_peak.dat
x         = Float64[401]
y         = Float64[401]
staterror = None
syserror  = None
None

This could be displayed with matplotlib directly, but let’s create a DataPlot to display it:

>>> from sherpa.plot import DataPlot
>>> dplot = DataPlot()
>>> dplot.prepare(d)
>>> dplot.plot()
../_images/data.png

Setting up the fits

We can then fit this with a variety of models, using the LeastSq statistic and NelderMead optimiser. For each model - NormGauss1D, Lorentz1D, PseudoVoigt1D, and Voigt1D - we create a Fit instance:

>>> from sherpa.models.basic import NormGauss1D
>>> from sherpa.astro.models import Lorentz1D, PseudoVoigt1D, Voigt1D
>>> from sherpa.stats import LeastSq
>>> from sherpa.optmethods import NelderMead
>>> from sherpa.fit import Fit
>>> models = [NormGauss1D, Lorentz1D, PseudoVoigt1D, Voigt1D]
>>> stat = LeastSq()
>>> method = NelderMead()
>>> fits = {}
>>> for model in models:
...     mdl = model()
...     fit = Fit(d, mdl, stat, method)
...     fits[mdl.name] = fit

The fits dictionary is indexed by the model name which, in this case, uses the default model names (as none were given when initializing the model), which is the lower-case value of the model. So the keys are 'normgauss1d', 'lorentz1d', 'pseudovoigt1d', and 'voigt1d'.

Fitting the NormGauss1D and Lorentz1D models

The fits could have been done in this loop, but let’s run them separately. First the NormGauss1D and Lorentz1D models, which bracket the data, using the format() method to display the fit results:

>>> result = fits['normgauss1d'].fit()
>>> print(result.format())
Method                = neldermead
Statistic             = leastsq
Initial fit statistic = 4222.73
Final fit statistic   = 29.9943 at function evaluation 439
Data points           = 401
Degrees of freedom    = 398
Change in statistic   = 4192.74
   normgauss1d.fwhm   2.90157
   normgauss1d.pos   9.24277
   normgauss1d.ampl   30.3135
>>> result = fits['lorentz1d'].fit()
>>> print(result.format())
Method                = neldermead
Statistic             = leastsq
Initial fit statistic = 4213.59
Final fit statistic   = 53.7535 at function evaluation 401
Data points           = 401
Degrees of freedom    = 398
Change in statistic   = 4159.84
   lorentz1d.fwhm   2.30967
   lorentz1d.pos   9.24439
   lorentz1d.ampl   38.9727

For reference, these results are essentially the same as obtained by lmfit.

We can create ModelPlot objects for each model (extracting the model from the fit structure as I forgot to save it earlier):

>>> from sherpa.plot import ModelPlot
>>> mplot1 = ModelPlot()
>>> mplot2 = ModelPlot()
>>> mplot1.prepare(d, fits['normgauss1d'].model)
>>> mplot2.prepare(d, fits['lorentz1d'].model)

These can be used to compare the two model plots (I added in some matplotlib commands to set the plot legend):

>>> dplot.plot(alpha=0.5)
>>> ax = plt.gca()
>>> ax.lines[-1].set_label('data')
>>> mplot1.overplot()
>>> ax.lines[-1].set_label('NormGauss1D')
>>> mplot2.overplot()
>>> ax.lines[-1].set_label('Lorentz1D')
>>> plt.legend()
../_images/models-normgauss1d-lorentz1d.png

Fitting the PseudoVoigt1D and Voigt1D models

So neither model describe the data at the peak or the wings. Let’s see if the extra freedom given by the PseudoVoigt1D and Voigt1D models can help:

>>> result = fits['pseudovoigt1d'].fit()
>>> print(result.format())
Method                = neldermead
Statistic             = leastsq
Initial fit statistic = 4220.01
Final fit statistic   = 10.8366 at function evaluation 851
Data points           = 401
Degrees of freedom    = 397
Change in statistic   = 4209.18
   pseudovoigt1d.frac   0.533597
   pseudovoigt1d.fwhm   2.71357
   pseudovoigt1d.pos   9.2436
   pseudovoigt1d.ampl   34.4984

The Voigt1D fit here should be compared to the “gamma unconstrained” version from the lmfit example.

>>> result = fits['voigt1d'].fit()
>>> print(result.format())
Method                = neldermead
Statistic             = leastsq
Initial fit statistic = 4211.75
Final fit statistic   = 10.9302 at function evaluation 738
Data points           = 401
Degrees of freedom    = 397
Change in statistic   = 4200.82
   voigt1d.fwhm_g   2.10801
   voigt1d.fwhm_l   1.0508
   voigt1d.pos    9.24375
   voigt1d.ampl   34.1915

The first thing to note is that the final statistic value, \(\sim 11\), is much lower than either of the previous models, so they should be a better fit, and that the PseudoVoigt1D model is “better” than the Voigt1D model, but not by much. Since the two models appear similar let’s create two plots:

>>> from sherpa.plot import SplitPlot
>>> mplot3 = ModelPlot()
>>> mplot4 = ModelPlot()
>>> mplot3.prepare(d, fits['voigt1d'].model)
>>> mplot4.prepare(d, fits['pseudovoigt1d'].model)
>>> splot = SplitPlot()
>>> splot.addplot(dplot, alpha=0.5)
>>> splot.overlayplot(mplot3)
>>> ax = plt.gca()
>>> ax.set_title('Voigt1D')
>>> ax.set_xlabel('')
>>> splot.addplot(dplot, alpha=0.5)
>>> splot.overlayplot(mplot4)
>>> ax = plt.gca()
>>> ax.set_title('PseudoVoigt1D')
../_images/models-voigts.png

Comparing the models

We can also overlay the models (here only showing one of the Voigt profiles since on this plot they appear identical):

>>> mplot1.plot()
>>> ax = plt.gca()
>>> ax.lines[-1].set_label('NormGauss1D')
>>> mplot2.overplot()
>>> ax.lines[-1].set_label('Lorentz1D')
>>> mplot3.overplot()
>>> ax.lines[-1].set_label('Voigt1D')
>>> plt.legend()
../_images/models.png