Sherpa: Modeling and Fitting in Python¶
Sherpa is a modeling and fitting application for Python. It contains a powerful language for combining simple models into complex expressions that can be fit to the data using a variety of statistics and optimization methods. It is easily extensible to include user models, statistics and optimization methods.
What can you do with Sherpa?¶
- fit 1D (multiple) data including: spectra, surface brightness profiles, light curves, general ASCII arrays
- fit 2D images/surfaces in Poisson/Gaussian regime
- build complex model expressions
- import and use your own models
- use appropriate statistics for modeling Poisson or Gaussian data
- import the new statistics, with priors if required by analysis
- visualize a parameter space with simulations or using 1D/2D cuts of the parameter space
- calculate confidence levels on the best fit model parameters
- choose a robust optimization method for the fit: Levenberg-Marquardt, Nelder-Mead Simplex or Monte Carlo/Differential Evolution.
To install Sherpa see Install Section. For an example of a Sherpa session check run-sherpa section. More examples are given in a few Sherpa IPython Notebooks:
For detailed documentation see: http://cxc.harvard.edu/sherpa
Sherpa 4.8.2 was released on September 22, 2016. This is the first release for both Python 2.7 and Python 3.5 (Beta release). The 4.8.2 release is available for Python 2.7 and for the first time (as a beta release with some caveats) for Python 3.5 The Python 3.5 release is considered a Beta designed to maintain backwards compatibility with Python 2.7.
The complete 4.8.2 Release Notes can be found here.
The code has been tested on Linux and Mac OSX (versions > 10.8)
Sherpa can be installed from a binary distribution or built from sources. The binary distribution is suited for people wanting to have Sherpa up and running as soon as possible in its standard form. The binaries are built and tested on Linux and Mac OSX (>=10.8)
Source installation is available for platforms incompatible with the binary builds, or for users wanting to customize the way Sherpa is built and installed.
The code and installation instructions are available on GitHub.
Sherpa expects a user configuration file sherpa-standalone.rc to be in the $HOME directory. If this file is not present Sherpa will use the default internal configuration file with the IO and Plotting back-ends set to pyfits and pylab.
matplotlib comes with a configuration file matplotlibrc. For smooth behavior with Sherpa, be sure to indicate interactive=True in ~/.matplotlib/matplotlibrc.
You can import Sherpa into your ipython session:
(conda)$ ipython --pylab Python 2.7.12 |Continuum Analytics, Inc.| (default, Jul 2 2016, 17:43:17) Type "copyright", "credits" or "license" for more information. IPython 5.1.0 -- An enhanced Interactive Python. ? -> Introduction and overview of IPython's features. %quickref -> Quick reference. help -> Python's own help system. object? -> Details about 'object', use 'object??' for extra details. Using matplotlib backend: MacOSX In : from sherpa.astro.ui import * WARNING: imaging routines will not be available, failed to import sherpa.image.ds9_backend due to 'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH'
The standard warnings are issued if you do not have ds9 models in your path. The image with ds9 will not be available. See the Dependencies section below.
Now to simulate a simple shape (a parabola with errors):
In : x = np.arange(-5, 5.1) In : y = x*x + 23.2 + np.random.normal(size=x.size) In : e = np.ones(x.size)
The data can now be loaded into Sherpa:
In : load_arrays(1, x, y, e) In : plot_data()
For this example we know what model to use, so pick a polynomial and free-up some of the parameters:
In : set_source(polynom1d.poly) In : print(poly) polynom1d.poly Param Type Value Min Max Units ----- ---- ----- --- --- ----- poly.c0 thawed 1 -3.40282e+38 3.40282e+38 poly.c1 frozen 0 -3.40282e+38 3.40282e+38 poly.c2 frozen 0 -3.40282e+38 3.40282e+38 poly.c3 frozen 0 -3.40282e+38 3.40282e+38 poly.c4 frozen 0 -3.40282e+38 3.40282e+38 poly.c5 frozen 0 -3.40282e+38 3.40282e+38 poly.c6 frozen 0 -3.40282e+38 3.40282e+38 poly.c7 frozen 0 -3.40282e+38 3.40282e+38 poly.c8 frozen 0 -3.40282e+38 3.40282e+38 poly.offset frozen 0 -3.40282e+38 3.40282e+38 In : thaw(poly.c1, poly.c2)
With everything set up, the data can be fit using the standard optimization method levmar and chi2 statistics:
In : fit() Dataset = 1 Method = levmar Statistic = chi2 Initial fit statistic = 12190 Final fit statistic = 5.40663 at function evaluation 8 Data points = 11 Degrees of freedom = 8 Probability [Q-value] = 0.713361 Reduced statistic = 0.675829 Change in statistic = 12184.6 poly.c0 22.2341 poly.c1 0.109262 poly.c2 1.06812 In : plot_fit_resid()
and an estimate of 1 sigma parameters uncertainties is:
In : conf() poly.c0 lower bound: -0.455477 poly.c1 lower bound: -0.0953463 poly.c0 upper bound: 0.455477 poly.c2 lower bound: -0.0341394 poly.c1 upper bound: 0.0953463 poly.c2 upper bound: 0.0341394 Dataset = 1 Confidence Method = confidence Iterative Fit Method = None Fitting Method = levmar Statistic = chi2gehrels confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- poly.c0 22.2341 -0.455477 0.455477 poly.c1 0.109262 -0.0953463 0.0953463 poly.c2 1.06812 -0.0341394 0.0341394
Data I/O support and plotting need astropy (pyFITS is deprecated) and matplotlib. Imaging requires both ds9 and XPA.
|[mpl]||Hunter, JD (2007). Matplotlib: A 2D graphics environment. Computing in Science and Engineering. 9: 90-95. http://matplotlib.sourceforge.net.|
Caveats for Python 3.5¶
show_all does not work on Python3 with PHA. E.g., the following code throws an exception:from sherpa.astro.ui.utils import Session from sherpa.astro.data import DataPHA session = Session() session.load_arrays(1, [1, 2, 3], [1, 2, 3], DataPHA) session.show_all()
the parallel_map function throws an exception on Python 3 when only one core is available on the system running Sherpa. This impacts a number of functions taking advantage of tasks parallelization on multi-core systems. This functions include calc_flux, sample_flux, a number of functions for plotting parameter projections, and the grid_search optimization algorithm.
string representation of some (but not all) dataset classes is broken on Python 3. For instance, the following code will result in an exception:from sherpa.astro.ui import * load_pha('3c273.pi') print(get_data())
sample_flux does not work on Python 3 when called with the scales argument, e.g.from sherpa.astro import ui load_pha('sherpa-test-data/sherpatest/3c273.pi') set_model('polynom1d.p1') fit() covar() scal = get_covar_results().parmaxes sample_flux(p1, 0.5, 1, num=5, correlated=False, scales=scal)