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()
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()
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')
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()