Usage details¶
Now that you have the basics let’s move on to some more complex usage of the fitter interface.
First a quick preamble to do some imports and create our SherpaFitter
object.
from saba import SherpaFitter
sfit = SherpaFitter(statistic='chi2', optimizer='levmar', estmethod='confidence')
from astropy.modeling.models import Gaussian1D
import numpy as np
np.random.seed(0x1337)
Parameter constraints¶
Parameter constraints can be set in astropy models, and those constraints are
taken into account during the fitting procedure. Let’s take a quick look at
that. Firstly, let’s make a compound model by adding two Gaussian1D
instances.
After that, let’s set the following constraint on the parameter amplitude_1
(the amplitude of the right hand side Gaussian1D
):
1.2*amplitude_0
.
In addition, let’s create some synthetic data:
from astropy.modeling.models import Gaussian1D
import numpy as np
np.random.seed(0x1337)
double_gaussian = Gaussian1D(
amplitude=10, mean=-1.5, stddev=0.5) + Gaussian1D(amplitude=1, mean=0.9,
stddev=0.5)
def tiedfunc(obj): # a function used for tying amplitude_1
return 1.2 * obj.amplitude_0
double_gaussian.amplitude_1.tied = tiedfunc
double_gaussian.amplitude_1.value = double_gaussian.amplitude_1.tied(
double_gaussian)
err = 0.8
step = 0.2
x = np.arange(-3, 3, step)
y = double_gaussian(x) + err * np.random.uniform(-1, 1, size=len(x))
yerrs = err * np.random.uniform(0.2, 1, size=len(x))
# please note these are binsize/2 not true errors!
binsize = (step / 2) * np.ones(x.shape)
plt.errorbar(x, y, xerr=binsize, yerr=yerrs, ls="", label="data")
# once again xerrs are binsize/2 not true errors!
plt.plot(x, double_gaussian(x), label="True")
plt.legend(loc=(0.78, 0.8), frameon=False)
plt.xlim((-3, 3))
(Source code, png, hires.png, pdf)
Note
without astropy PR #5129 we need to do this
double_gaussian.amplitude_1.value = double_gaussian.amplitude_1.tied(double_gaussian)
Let’s set some more parameter constraints to the model and fit the data. We can print the sherpa models to check things are doing what they should.
fit_gg = double_gaussian.copy()
fit_gg.mean_0.value = -0.5
# sets the lower bound so we can force the parameter against it
fit_gg.mean_0.min = -1.25
fit_gg.mean_1.value = 0.8
fit_gg.stddev_0.value = 0.9
fit_gg.stddev_0.fixed = True
Fitting Config¶
A SherpaFitter
object has the opt_config
property which holds the
configuration details for the optimization routine. Its docstring contains
information about the the properties of the optimizer.
We can see those configuration details by using print(sfit.opt_config)
which outputs:
{'epsfcn': 1.1920928955078125e-07,
'factor': 100.0,
'ftol': 1.1920928955078125e-07,
'gtol': 1.1920928955078125e-07,
'maxfev': None,
'verbose': 0,
'xtol': 1.1920928955078125e-07}
Similarly for print(sfit.opt_config.__doc__)
:
Levenberg-Marquardt optimization method.
The Levenberg-Marquardt method is an interface to the MINPACK
subroutine lmdif to find the local minimum of nonlinear least
squares functions of several variables by a modification of the
Levenberg-Marquardt algorithm [1]_.
Attributes
----------
ftol : number
The function tolerance to terminate the search for the minimum;
the default is sqrt(DBL_EPSILON) ~ 1.19209289551e-07, where
DBL_EPSILON is the smallest number x such that `1.0 != 1.0 +
x`. The conditions are satisfied when both the actual and
predicted relative reductions in the sum of squares are, at
most, ftol.
xtol : number
The relative error desired in the approximate solution; default
is sqrt( DBL_EPSILON ) ~ 1.19209289551e-07, where DBL_EPSILON
is the smallest number x such that `1.0 != 1.0 + x`. The
conditions are satisfied when the relative error between two
consecutive iterates is, at most, `xtol`.
...
The parameters can be changed as
sfit.opt_config['ftol'] = 1e-5
print(sfit.opt_config)
{'epsfcn': 1.1920928955078125e-07,
'factor': 100.0,
'ftol': 1e-05,
'gtol': 1.1920928955078125e-07,
'maxfev': None,
'verbose': 0,
'xtol': 1.1920928955078125e-07}
Fitting this model is similar as showing previously. For the sake of comparison let’s also fit and unconstrained model:
fitted_gg = sfit(fit_gg, x, y, xbinsize=binsize, err=yerrs)
sfit_unconstrained = SherpaFitter(statistic='chi2', optimizer='levmar',
estmethod='covariance')
free_gg = sfit_unconstrained(double_gaussian.copy(), x, y,
xbinsize=binsize, err=yerrs)
plt.figure(figsize=(10, 5))
plt.plot(x, double_gaussian(x), label="True")
plt.errorbar(x, y, xerr=binsize, yerr=yerrs, ls="", label="data")
plt.plot(x, fit_gg(x), label="Pre fit")
plt.plot(x, fitted_gg(x), label="Fitted")
plt.plot(x, free_gg(x), label="Free")
plt.subplots_adjust(right=0.8)
plt.legend(loc=(1.01, 0.55), frameon=False)
plt.xlim((-3, 3))
(Source code, png, hires.png, pdf)
The fitter keeps a copy of the converted model so we can use it to compare the constrained and unconstrained model setups:
Note
wrap\_.amplitude_1
should be linked
, sherpa notation of astropy’s tied
wrap\_.stddev_0
should be frozen
, sherpa notation for fixed
and finally wrap\_.mean_0
’s value should have moved to its minimum while fitting
“wrap_” is just perpended to the model name (we didn’t set one so it’s blank) on conversion to the sherpa Model
.
print("##Fit with constraints")
print(sfit._fitmodel.sherpa_model)
print("##Fit without constraints")
print(sfit_unconstrained._fitmodel.sherpa_model)
##Fit with constraints
Param Type Value Min Max Units
----- ---- ----- --- --- -----
wrap_.amplitude_0 thawed 5.58947 -3.40282e+38 3.40282e+38
wrap_.mean_0 thawed -1.25 -1.25 3.40282e+38
wrap_.stddev_0 frozen 0.9 -3.40282e+38 3.40282e+38
wrap_.amplitude_1 linked 6.70736 expr: (1.2 * wrap_.amplitude_0)
wrap_.mean_1 thawed 0.869273 -3.40282e+38 3.40282e+38
wrap_.stddev_1 thawed 0.447021 -3.40282e+38 3.40282e+38
##Fit without constraints
Param Type Value Min Max Units
----- ---- ----- --- --- -----
wrap_.amplitude_0 thawed 6.95483 -3.40282e+38 3.40282e+38
wrap_.mean_0 thawed -1.59091 -3.40282e+38 3.40282e+38
wrap_.stddev_0 thawed 0.545582 -3.40282e+38 3.40282e+38
wrap_.amplitude_1 linked 8.34579 expr: (1.2 * wrap_.amplitude_0)
wrap_.mean_1 thawed 0.785016 -3.40282e+38 3.40282e+38
wrap_.stddev_1 thawed 0.46393 -3.40282e+38 3.40282e+38
Error Estimation Configuration¶
As with the optmethods
before we are able to adjust the configuration of the estmethods
. Some of the properties can be passed through est_errors
as keyword arguments such as the sigma
however for access to all options we have the est_config
property.
print(sfit.est_config)
sfit.est_config['numcores'] = 5
sfit.est_config['max_rstat'] = 4
print(sfit.est_config)
{'eps': 0.01,
'fast': False,
'max_rstat': 3,
'maxfits': 5,
'maxiters': 200,
'numcores': 8,
'openinterval': False,
'parallel': True,
'remin': 0.01,
'sigma': 1,
'soft_limits': False,
'tol': 0.2,
'verbose': False}
{'eps': 0.01,
'fast': False,
'max_rstat': 3,
'maxfits': 5,
'maxiters': 200,
'numcores': 5,
'openinterval': False,
'parallel': True,
'remin': 0.01,
'sigma': 1,
'soft_limits': False,
'tol': 0.2,
'verbose': False}
Multiple models or multiple datasets¶
We have three scenarios we can handle:
Fitting
N
datasets withN
modelsFitting a single dataset with
N
modelsFitting
N
datasets with a single model
If N > 1
for any of the scenarios then calling the fitter will return a list of models. Firstly we look at a single dataset with the two models as above.
We quickly copy the two models above and supply them to the fitter as a list - hopefully we get the same result.
fit_gg = double_gaussian.copy()
fit_gg.mean_0.value = -0.5
fit_gg.mean_0.min = -1.25
fit_gg.mean_1.value = 0.8
fit_gg.stddev_0.value = 0.9
fit_gg.stddev_0.fixed = True
fm1, fm2 = sfit([fit_gg, double_gaussian.copy()], x, y, xbinsize=binsize, err=yerrs)
(Source code, png, hires.png, pdf)
We also can fit multiple datasets with a single model so let’s make a second dataset:
second_gg = double_gaussian.copy()
second_gg.mean_0 = -2
second_gg.mean_1 = 0.5
second_gg.amplitude_0 = 8
second_gg.amplitude_1 = 5
second_gg.stddev_0 = 0.4
second_gg.stddev_1 = 0.8
y2 = second_gg(x) + err * np.random.uniform(-1, 1, size=len(x))
y2errs = err * np.random.uniform(0.2, 1, size=len(x))
Here we supply lists for each of the data parameters. You can also use None
for when you don’t have something like a missing binsizes - a lack of binsizes is a contrived example but a lack of y
errors is not suitable for a chi:sup:2 fit and you don’t want to make a new fitter.
fit_gg = double_gaussian.copy()
fit_gg.mean_0 = -2.3
fit_gg.mean_1 = 0.7
fit_gg.amplitude_0 = 2
fit_gg.amplitude_1 = 3
fit_gg.stddev_0 = 0.3
fit_gg.stddev_1 = 0.5
fm1, fm2 = sfit(fit_gg, x=[x, x], y=[y, y2], xbinsize=[binsize, None], err=[yerrs, y2errs])
(Source code, png, hires.png, pdf)
Background Data¶
It is also possible specify background data which is required for several of the fit statistics.
This is done by supplying a background array using the bkg
keyword. If there is a scaling of the background relative to the source data then you can use the bkg_scale
keyword
y[y<0]=0
cfit = SherpaFitter(statistic='cstat', optimizer='levmar', estmethod='covariance')
cfit(fit_gg, x=x, y=y, xbinsize=binsize, err=yerrs, bkg=y, bkg_scale=0.3)
(Source code, png, hires.png, pdf)