.. include:: references.txt Saba: Sherpa-Astropy Bridge =========================== The Saba package provides a bridge between the convenient model definition language provided in the `astropy.modeling` package and the powerful fitting capabilities of the Sherpa_ modeling and fitting package. In particular, Sherpa has a selection of robust optimization algorithms coupled with configurable fit statistics. Once the model fit is complete Sherpa has three different ways to estimate parameter confidence intervals, including methods that allow for coupled non-gaussian errors. Finally, Sherpa has an MCMC sampler that can be used to generate draws from the probability distribution assuming given priors. Once Saba and Sherpa are installed, the Saba package exposes the above Sherpa functionality within the `astropy.modeling.fitting` package via a single |SherpaFitter| class which acts as a fitting backend within astropy. If using the latest version of astropy (development or >= 1.3), a plugin registry system automatically makes the |SherpaFitter| class available within the `astropy.modeling.fitting` module without requiring an explicit import. ``Saba`` is the Sherpa people's word for "bridge". Installation ------------ The following installation notes apply to the development version of Saba and assume use of the conda + Anaconda package system. Prerequisites ^^^^^^^^^^^^^ * Python 2.7 (support for Python 3.5+ is in work) * numpy * astropy * sherpa .. doctest-skip:: conda install numpy To make use of the entry points plugin registry which automatically makes the |SherpaFitter| class available within `astropy.modeling.fitting` install ``astropy`` version >= 1.3. Otherwise one can just use the latest stable ``astropy`` via:: conda install astropy Next install Sherpa_ using the conda ``sherpa`` channel. Note that Sherpa currently needs to be installed after astropy on Mac OSX. .. doctest-skip:: conda install -c sherpa sherpa Finally install ``saba`` using pip:: pip install saba Getting started --------------- If you are not already familiar with `astropy.modeling`, now is a good time to review the introductory documentation there along with the `astropy.modeling.fitting` module details. To start with Saba let's import the |SherpaFitter| class which is the interface with Sherpa's fitting routines. |SherpaFitter| is available in one of two ways, either directly from `saba` or through `astropy.modeling.fitting` through the plugin registry system. The latter method is preferred but requires `astropy` version >= 1.3 or the latest development (master) version. Use:: from saba import SherpaFitter or .. doctest-skip:: from astropy.modeling.fitting import SherpaFitter Initialization ^^^^^^^^^^^^^^ To initialize a fitter we provide string values to define the ``statistic``, ``optimizer`` and ``estmethod`` (error estimation method). The available values for those can be found in the docstring of |SherpaFitter| and relate to objects within `sherpa.stats`, `sherpa.optmethods` and `sherpa.estmethods`. .. doctest-skip:: sfit = SherpaFitter(statistic='chi2', optimizer='levmar', estmethod='confidence') Now that we have a fitter instance we need something to fit. So let's import an astropy model, specifically `~astropy.modeling.functional_models.Gaussian1D`. A full description astropy's model and capabilities can be found in the ``astropy`` `Instantiating and Evaluating Models `_ section. .. doctest-skip:: from astropy.modeling.models import Gaussian1D We also need some data so let's make some data with some added noise. .. plot:: :include-source: import numpy as np from astropy.modeling.models import Gaussian1D, Gaussian2D import matplotlib.pyplot as plt np.random.seed(0x1337) true = Gaussian1D(amplitude=3, mean=0.9, stddev=0.5) err = 0.05 step = 0.1 x = np.arange(-3, 3, step) y = true(x) + err * np.random.uniform(-1, 1, size=len(x)) yerrs = err * np.random.uniform(0.2, 1, size=len(x)) binwidth = step * np.ones(x.shape) fit_model = true.copy() # ofset fit model from true fit_model.amplitude = 2 fit_model.mean = 0 fit_model.stddev = 0.2 plt.plot(x, true(x), label="True") plt.errorbar(x, y, xerr=binwidth, yerr=yerrs, ls="", label="Data") plt.plot(x, fit_model(x), label="Starting fit model") plt.legend(loc=2, frameon=False) plt.xlim((-3, 3)) plt.show() Fitting ^^^^^^^ Now we have some data let's fit it and hopefully we get something similar to "True" back. The ``sfit`` fitter object has already been initialized (as would be done for other `astropy.modeling.fitting` fitters) so we just call it with some data and an astropy model and we get the fitted model returned. .. doctest-skip:: fitted_model = sfit(fit_model, x, y, xbinsize=binwidth, err=yerrs) plt.plot(x, true(x), label="True") plt.errorbar(x, y, xerr=binwidth, yerr=yerrs, ls="", label="Data") plt.plot(x, fit_model(x), label="Starting fit model") plt.plot(x, fitted_model(x), label="Fitted model") plt.legend(loc=2, frameon=False) plt.xlim((-3, 3)) Once again plotting the data. .. plot:: from saba import SherpaFitter from astropy.modeling.models import Gaussian1D, Gaussian2D import numpy as np import matplotlib.pyplot as plt sfitter = SherpaFitter(statistic='chi2', optimizer='levmar', estmethod='confidence') sfitter.est_config['max_rstat'] = 4 np.random.seed(0x1337) true = Gaussian1D(amplitude=3, mean=0.9, stddev=0.5) err = 0.05 step = 0.1 x = np.arange(-3, 3, step) y = true(x) + err * np.random.uniform(-1, 1, size=len(x)) yerrs = err * np.random.uniform(0.2, 1, size=len(x)) binwidth = step * np.ones(x.shape) fit_model = true.copy() # ofset fit model from true fit_model.amplitude = 2 fit_model.mean = 0 fit_model.stddev = 0.2 fitted_model = sfitter(fit_model, x, y, xbinsize=binwidth, err=yerrs) plt.plot(x, true(x), label="True") plt.errorbar(x, y, xerr=binwidth, yerr=yerrs, ls="", label="Data") plt.plot(x, fit_model(x), label="Starting fit model") plt.plot(x, fitted_model(x), label="Fitted model") plt.legend(loc=2, frameon=False) plt.xlim((-3, 3)) plt.show() Now we have a fit we can look at the outputs by doing: .. doctest-skip:: print(sfit.fit_info) .. doctest-skip:: datasets = None itermethodname = none methodname = levmar statname = chi2 succeeded = True parnames = ('wrap_.amplitude', 'wrap_.mean', 'wrap_.stddev') parvals = (3.0646789274093185, 0.77853851419777986, 0.50721937454701504) statval = 82.7366242121 istatval = 553.030876852 dstatval = 470.29425264 numpoints = 30 dof = 27 qval = 1.44381192266e-07 rstat = 3.06431941526 message = successful termination nfev = 84 Note that the ``fit_info`` attribute is custom to the |SherpaFitter| class and provides a direct link to the internal fitting results from the Sherpa fit process. Uncertainty estimation ^^^^^^^^^^^^^^^^^^^^^^ One of the main drivers for Saba is to get access the uncertainty estimation methods provided by Sherpa_. This is done through the `~saba.SherpaFitter.est_errors` method which uses the Sherpa's `~sherpa.fit.Fit.est_errors` method. To get the errors make a call such as: .. doctest-skip:: param_errors = sfit.est_errors(sigma=3) # Note that sigma can be an input In return we get a tuple of (``parameter_name``, ``best_fit_value``, ``lower_value`` , ``upper_value``). For the sake of plotting them we make models for the upper and lower values, and then output the values while we're at it. .. doctest-skip:: min_model = fitted_model.copy() max_model = fitted_model.copy() for pname, pval, pmin, pmax in list(zip(*param_errors)): print(pname, pval, pmin, pmax) getattr(min_model, pname).value = pval + pmin getattr(max_model, pname).value = pval + pmax .. doctest-skip:: ('amplitude', 3.0646789274093185, -0.50152026852144349, 0.56964617033348119) ('mean', 0.77853851419777986, -0.096264447380365548, 0.10293940565584792) ('stddev', 0.50721937454701504, -0.098092469817728456, 0.11585973498734969) .. plot:: from saba import SherpaFitter from astropy.modeling.models import Gaussian1D, Gaussian2D import numpy as np import matplotlib.pyplot as plt sfitter = SherpaFitter(statistic='chi2', optimizer='levmar', estmethod='confidence') sfitter.est_config['max_rstat'] = 4 np.random.seed(0x1337) true = Gaussian1D(amplitude=3, mean=0.9, stddev=0.5) err = 0.05 step = 0.1 x = np.arange(-3, 3, step) y = true(x) + err * np.random.uniform(-1, 1, size=len(x)) yerrs = err * np.random.uniform(0.2, 1, size=len(x)) binwidth = step * np.ones(x.shape) fit_model = true.copy() # ofset fit model from true fit_model.amplitude = 2 fit_model.mean = 0 fit_model.stddev = 0.2 fitted_model = sfitter(fit_model, x, y, xbinsize=binwidth, err=yerrs) param_errors = sfitter.est_errors(sigma=3) min_model = fitted_model.copy() max_model = fitted_model.copy() for pname, pval, pmin, pmax in list(zip(*param_errors)): getattr(min_model, pname).value = pval + pmin getattr(max_model, pname).value = pval + pmax plt.plot(x, true(x), label="True") plt.errorbar(x, y, xerr=binwidth, yerr=yerrs, ls="") plt.plot(x, fitted_model(x), label="Fitted model") plt.plot(x, min_model(x), label="min model", ls="--") plt.plot(x, max_model(x), label="max model", ls="--") plt.legend(loc=2, frameon=False) plt.xlim((-3, 3)) plt.show() Using Saba ---------- .. toctree:: :maxdepth: 2 examples_complex examples_2d examples_mcmc API/Reference ------------- .. toctree:: :maxdepth: 2 api.rst Credit ------ The development of this package was made possible by the generous support of the `Google Summer of Code `_ program in 2016 under the `OpenAstronomy `_ by `Michele Costa `_ with the support and advice of mentors `Tom Aldcroft `_, `Omar Laurino `_, `Moritz Guenther `_, and `Doug Burke `_.