Welcome to Sherpa’s documentation¶
Welcome to the Sherpa documentation. Sherpa is a Python package for modeling and fitting data. It was originally developed by the Chandra X-ray Center for use in analysing X-ray data (both spectral and imaging) from the Chandra X-ray telescope, but it is designed to be a general-purpose package, which can be enhanced with domain-specific tasks (such as X-ray Astronomy). Sherpa contains an expressive and powerful modeling language, coupled with a range of statistics and robust optimisers.
Sherpa is released under the GNU General Public License v3.0, and is compatible with Python versions 2.7, 3.5, 3.6, and 3.7. Information on recent releases and citation information for Sherpa is available using the Digital Object Identifier (DOI) 10.5281/zenodo.593753.
Installation¶
Quick overview¶
For those users who have already read this page, and need a quick refresher (or prefer to act first, and read documentation later), the following commands can be used to install Sherpa, depending on your environment and set up.
conda install -c sherpa sherpa
pip install sherpa
python setup.py install
Requirements¶
Sherpa has the following requirements:
- Python 2.7, 3.5, 3.6, or 3.7
- NumPy (the exact lower limit has not been determined, but it is likely to be 1.7.0 or later)
- Linux or OS-X (patches to add Windows support are welcome)
Sherpa can take advantage of the following Python packages if installed:
- astropy: for reading and writing files in FITS format. The minimum required version of astropy is version 1.3, although only versions 2 and higher are used in testing.
- matplotlib: for visualisation of one-dimensional data or models, one- or two- dimensional error analysis, and the results of Monte-Carlo Markov Chain runs. There are no known incompatabilities with matplotlib, but there has only been limited testing. Please report any problems you find.
The Sherpa build can be configured to create the
sherpa.astro.xspec module, which provides the models and utility
functions from the XSPEC.
The supported versions of XSPEC are 12.10.1 (patch level a or later),
12.10.0, 12.9.1, and 12.9.0.
Interactive display and manipulation of two-dimensional images is available if the DS9 image viewer and the XPA commands are installed. It is expected that any recent version of DS9 can be used.
Releases and version numbers¶
The Sherpa release policy has a major release at the start of the year, corresponding to the code that is released in the previous December as part of the CIAO release, followed by several smaller releases throughout the year.
Information on the Sherpa releases is available from the Zenodo page for Sherpa, using the Digital Object Identifier (DOI) 10.5281/zenodo.593753.
What version of Sherpa is installed?¶
The version number and git commit id of Sherpa can be retrieved from
the sherpa._version module using the following command:
% python -c 'import sherpa._version; print(sherpa._version.get_versions())'
{'version': '4.10.0', 'full': 'c7732043124b08d5e949b9a95c2eb6833e009421'}
Citing Sherpa¶
Information on citing Sherpa can be found from the CITATION document in the Sherpa repository, or from the Sherpa Zenodo page.
Installing a pre-compiled version of Sherpa¶
Additional useful Python packages include astropy, matplotlib,
and ipython-notebook.
Using the Anaconda python distribution¶
The Chandra X-ray Center provides releases of Sherpa that can be
installed using
Anaconda
from the sherpa channel. First check
to see what the latest available version is by using:
conda install -c sherpa sherpa --dry-run
and then, if there is a version available and there are no significant upgrades to the dependencies, Sherpa can be installed using:
conda install -c sherpa sherpa
It is strongly suggested that Sherpa is installed into a named conda environment (i.e. not the default environment).
Using pip¶
Sherpa is also available from PyPI at https://pypi.python.org/pypi/sherpa and can be installed with the command:
pip install sherpa
The NumPy package must already have been installed for this to work.
Building from source¶
Prerequisites¶
The prerequisites for building from source are:
- Python versions: 2.7, 3.5, 3.6, 3.7
- Python packages:
setuptools,numpy - System:
gcc,g++,make,flex,bison(the aim is to support recent versions of these tools; please report problems to the Sherpa issue tracker).
It is highly recommended that matplotlib and astropy be installed before building Sherpa, to avoid skipping a number of tests in the test suite.
The full Sherpa test suite requires the mock package (Python 2.7 only), pytest, and pytest-xvfb. These packages should be installed automatically for you by the test suite if they do not already exist.
Note
As of the Sherpa 4.10.1 release, a Fortran compiler is no-longer required to build Sherpa.
Obtaining the source package¶
The source code can be obtained as a release package from
Zenodo - e.g.
the Sherpa 4.10.0 release -
or from
the Sherpa repository on GitHub,
either a release version,
such as the
4.10.0 tag,
or the master branch (which is not guaranteed to be stable).
For example:
git clone git://github.com/sherpa/sherpa.git
cd sherpa
git checkout 4.10.0
will use the 4.10.0 tag.
Configuring the build¶
The Sherpa build is controlled by the setup.cfg file in the
root of the Sherpa source tree. These configuration options
include:
FFTW¶
Sherpa ships with the fftw library source
code and builds it by default. To use a different version, change
the fftw options in the sherpa_config section of the
setup.cfg file. The options to change are:
fftw=local
fftw-include_dirs=/usr/local/include
fftw-lib-dirs=/use/local/lib
fftw-libraries=fftw3
The fftw option must be set to local and then the remaining
options changed to match the location of the local installation.
XSPEC¶
Note
The version number of XSPEC must be specified using the
xspec_version configuration option, as described below. This is
a change from previous releases of Sherpa, but is required in order
to support changes made in XSPEC 12.10.0.
Sherpa can be built to use the Astronomy models provided by
XSPEC versions 12.10.1 (patch level a or later), 12.10.0,
12.9.1, and 12.9.0. To enable XSPEC support, several changes must be
made to the xspec_config section of the setup.cfg file. The
available options (with default values) are:
with-xspec = False
xspec_version = 12.9.0
xspec_lib_dirs = None
xspec_include_dirs = None
xspec_libraries = XSFunctions XSModel XSUtil XS
cfitsio_lib_dirs = None
cfitsio_libraries = cfitsio
ccfits_lib_dirs = None
ccfits_libraries = CCfits
wcslib_lib_dirs = None
wcslib_libraries = wcs
gfortran_lib_dirs = None
gfortran_libraries = gfortran
To build the sherpa.astro.xspec module, the
with-xspec option must be set to True and the
xspec_version option set to the correct version string (the XSPEC
patch level must not be included), and then the
remaining options depend on the version of XSPEC and whether
the XSPEC model library or the full XSPEC system has been installed.
In the examples below, the $HEADAS value must be replaced
by the actual path to the HEADAS installation, and the versions of
the libraries - such as CCfits_2.5 - may need to be changed to
match the contents of the XSPEC installation.
If the full XSPEC 12.10.1 system has been built then use:
with-xspec = True xspec_version = 12.10.1 xspec_lib_dirs = $HEADAS/lib xspec_include_dirs = $HEADAS/include xspec_libraries = XSFunctions XSUtil XS hdsp_6.25 ccfits_libraries = CCfits_2.5 wcslib_libraries = wcs-5.19.1
where the version numbers were taken from version 6.25 of HEASOFT and may need updating with a newer release.
If the full XSPEC 12.10.0 system has been built then use:
with-xspec = True xspec_version = 12.10.0 xspec_lib_dirs = $HEADAS/lib xspec_include_dirs = $HEADAS/include xspec_libraries = XSFunctions XSModel XSUtil XS hdsp_3.0 ccfits_libraries = CCfits_2.5 wcslib_libraries = wcs-5.16
If the full XSPEC 12.9.x system has been built then use:
with-xspec = True xspec_version = 12.9.1 xspec_lib_dirs = $HEADAS/lib xspec_include_dirs = $HEADAS/include xspec_libraries = XSFunctions XSModel XSUtil XS ccfits_libraries = CCfits_2.5 wcslib_libraries = wcs-5.16
changing
12.9.1to12.9.0as appropriate.If the model-only build of XSPEC has been installed, then the configuration is similar, but the library names may not need version numbers and locations, depending on how the
cfitsio,CCfits, andwcslibraries were installed.Note that XSPEC 12.10.0 introduces a new
--enable-xs-models-onlyflag when building HEASOFT which simplifies the installation of these extra libraries, but can cause problems for the Sherpa build.
A common problem is to set one or both of the xspec_lib_dirs
and xspec_lib_include options to the value of $HEADAS instead of
$HEADAS/lib and $HEADAS/include (after expanding out the
environment variable). Doing so will cause the build to fail with
errors about being unable to find various XSPEC libraries such as
XSFunctions and XSModel.
The gfortran options should be adjusted if there are problems
using the XSPEC module.
In order for the XSPEC module to be used from Python, the
HEADAS environment variable must be set before the
sherpa.astro.xspec module is imported.
The Sherpa test suite includes an extensive set of tests of this module, but a quick check of an installed version can be made with the following command:
% python -c 'from sherpa.astro import xspec; print(xspec.get_xsversion())'
12.10.1b
Warning
The --enable-xs-models-only flag with XSPEC 12.10.0 is known
to cause problems for Sherpa. It is strongly recommended that
either that the full XSPEC distribution is built, or that the
XSPEC installation from CIAO 4.11 is used.
Other options¶
The remaining options in the setup.cfg file allow Sherpa to be
built in specific environments, such as when it is built as part
of the CIAO analysis system. Please
see the comments in the setup.cfg file for more information on
these options.
Building and Installing¶
Note
It is highly recommended that some form of virtual environment, such as a conda environment or that provided by Virtualenv, be used when building and installing Sherpa.
A standard installation¶
From the root of the Sherpa source tree, Sherpa can be built by saying:
python setup.py build
and installed with one of:
python setup.py install
python setup.py install --user
A development build¶
The develop option should be used when developing Sherpa (such as
adding new functionality or fixing a bug):
python setup.py develop
Tests can then be run with the test option:
python setup.py test
The test command is a wrapper that calls pytest under the hood,
and includes the develop command.
You can pass additional arguments to pytest with the -a or
--pytest-args arguments. As an example, a single test can be run
using the syntax:
python setup.py test -a sherpa/astro/datastack/tests/test_datastack.py::test_load::test_case3
Note
If you run both install and develop or test in the same
Python environment you end up with two competing installations of
Sherpa which result in unexpected behavior. If this happens, simply
run pip uninstall sherpa as many times as necessary, until you
get an error message that no more Sherpa installations are
available. At this point you can re-install Sherpa.
The same issue may occur if you install a Sherpa binary release and then try to build Sherpa from source in the same environment.
The
Sherpa test data suite
can be installed to reduce the number of tests
that are skipped with the following (this is only for those builds
which used git to access the source code):
git submodule init
git submodule update
When both the DS9 image viewer and XPA toolset are installed, the test suite will include tests that check that DS9 can be used from Sherpa. This causes several copies of the DS9 viewer to be created, which can be distracting, as it can cause loss of mouse focus (depending on how X-windows is set up). This can be avoided by installing the X virtual-frame buffer (Xvfb).
Note
Although the standard Python setuptools approach is used to build
Sherpa, there may be issues when using some of the other build
targets, such as build_ext. Please report these to the
Sherpa issues page.
Building the documentation¶
Building the documentation requires the Sherpa source code and several additional packages:
- Python 3.5 or greater
- Sphinx, version 1.3 or later
- The
sphinx_rtd_theme - NumPy, six, and sphinx_astropy
- Graphviz (for the inheritance diagrams)
With these installed, the documentation can be built with the
build_sphinx target:
python setup.py build_sphinx
This can be done without building Sherpa (either an installation or development version), since Mock objects are used to represent compiled and optional components.
The documentation should be placed in build/sphinx/html/index.html,
although this may depend on what version of Sphinx is used.
It is also possible to build the documentation from within the docs/
directory:
cd docs
make html
This places the documentation in _build/html/index.html.
Testing the Sherpa installation¶
A very-brief “smoke” test can be run from the command-line with
the sherpa_smoke executable:
sherpa_smoke
WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available
----------------------------------------------------------------------
Ran 7 tests in 0.456s
OK (skipped=5)
or from the Python prompt:
>>> import sherpa
>>> sherpa.smoke()
WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available
----------------------------------------------------------------------
Ran 7 tests in 0.447s
OK (skipped=5)
This provides basic validation that Sherpa has been installed
correctly, but does not run many functional tests. The screen output
will include additional warning messages if the astropy or
matplotlib packages are not installed, or Sherpa was built
without support for the XSPEC model library.
The Sherpa installation also includes the sherpa_test command-line
tool which will run through the Sherpa test suite (the number of
tests depends on what optional packages are available and how
Sherpa was configured when built):
sherpa_test
The sherpa Anaconda channel contains the sherpatest package, which
provides a number of data files in ASCII and FITS formats. This is
only useful when developing Sherpa, since the package is large. It
will automatically be picked up by the sherpa_test script
once it is installed.
A quick guide to modeling and fitting in Sherpa¶
Here are some examples of using Sherpa to model and fit data. It is based on some of the examples used in the astropy.modelling documentation.
Getting started¶
The following modules are assumed to have been imported:
>>> import numpy as np
>>> import matplotlib.pyplot as plt
The basic process, which will be followed below, is:
- create a data object
- define the model
- select the statistic
- select the optimisation routine
- fit the data
- extract the parameter values
- Calculating error values
Although presented as a list, it is not necessarily a linear process, in that the order can be different to that above, and various steps can be repeated. The above list also does not include any visualization steps needed to inform and validate any choices.
Fitting a one-dimensional data set¶
The following data - where x is the independent axis and
y the dependent one - is used in this example:
>>> np.random.seed(0)
>>> x = np.linspace(-5., 5., 200)
>>> ampl_true = 3
>>> pos_true = 1.3
>>> sigma_true = 0.8
>>> err_true = 0.2
>>> y = ampl_true * np.exp(-0.5 * (x - pos_true)**2 / sigma_true**2)
>>> y += np.random.normal(0., err_true, x.shape)
>>> plt.plot(x, y, 'ko');
The aim is to fit a one-dimensional gaussian to this data and to recover
estimates of the true parameters of the model, namely the position
(pos_true), amplitude (ampl_true), and width (sigma_true).
The err_true term adds in random noise (using a
Normal distribution)
to ensure the data is not perfectly-described by the model.
Creating a data object¶
Rather than pass around the arrays to be fit, Sherpa has the
concept of a “data object”, which stores the independent and dependent
axes, as well as any related metadata. For this example, the
class to use is Data1D, which requires
a string label (used to identify the data), the independent
axis, and then dependent axis:
>>> from sherpa.data import Data1D
>>> d = Data1D('example', x, y)
>>> print(d)
name = example
x = Float64[200]
y = Float64[200]
staterror = None
syserror = None
At this point no errors are being used in the fit, so the staterror
and syserror fields are empty. They can be set either when the
object is created or at a later time.
Plotting the data¶
The sherpa.plot module provides a number of classes that
create pre-canned plots. For example, the
sherpa.plot.DataPlot class can be used to display the data.
The steps taken are normally:
- create the object;
- call the
prepare()method with the appropriate arguments, in this case the data object; - and then call the
plot()method.
Sherpa has two plotting backends: matplotlib, which is used by default for the standalone version, and ChIPS, which is used by CIAO. Limited support for customizing these plots - such as always drawing the Y axis with a logarithmic scale - is provided, but extensive changes will require calling the plotting back-end directly.
As an example of the DataPlot output:
>>> from sherpa.plot import DataPlot
>>> dplot = DataPlot()
>>> dplot.prepare(d)
>>> dplot.plot()
It is not required to use these classes and in the following, plots will be created either via these classes or directly via matplotlib.
Define the model¶
In this example a single model is used - a one-dimensional
gaussian provided by the Gauss1D
class - but more complex examples are possible: these
include multiple components,
sharing models between data sets, and
adding user-defined models.
A full description of the model language and capabilities is provided in
Creating model instances:
>>> from sherpa.models.basic import Gauss1D
>>> g = Gauss1D()
>>> print(g)
gauss1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
gauss1d.fwhm thawed 10 1.17549e-38 3.40282e+38
gauss1d.pos thawed 0 -3.40282e+38 3.40282e+38
gauss1d.ampl thawed 1 -3.40282e+38 3.40282e+38
It is also possible to restrict the range of a parameter, toggle parameters so that they are fixed or fitted, and link parameters together.
The sherpa.plot.ModelPlot class can be used to visualize
the model. The prepare() method
takes both a data object and the model to plot:
>>> from sherpa.plot import ModelPlot
>>> mplot = ModelPlot()
>>> mplot.prepare(d, g)
>>> mplot.plot()
There is also a sherpa.plot.FitPlot class which will
combine the two plot results,
but it is often just-as-easy to combine them directly:
>>> dplot.plot()
>>> mplot.overplot()
The model parameters can be changed - either manually or automatically - to try and start the fit off closer to the best-fit location, but for this example we shall leave the initial parameters as they are.
Select the statistics¶
In order to optimise a model - that is, to change the model parameters until the best-fit location is found - a statistic is needed. The statistic calculates a numerical value for a given set of model parameters; this is a measure of how well the model matches the data, and can include knowledge of errors on the dependent axis values. The optimiser (chosen below) attempts to find the set of parameters which minimises this statistic value.
For this example, since the dependent axis (y)
has no error estimate, we shall pick the least-square statistic
(LeastSq), which calculates the
numerical difference of the model to the data for each point:
>>> from sherpa.stats import LeastSq
>>> stat = LeastSq()
Select the optimisation routine¶
The optimiser is the part that determines how to minimise the statistic
value (i.e. how to vary the parameter values of the model to find
a local minimum). The main optimisers provided by Sherpa are
NelderMead
(also known as Simplex) and
LevMar
(Levenberg-Marquardt). The latter is often quicker, but less robust,
so we start with it (the optimiser can be changed and the data re-fit):
>>> from sherpa.optmethods import LevMar
>>> opt = LevMar()
>>> print(opt)
name = levmar
ftol = 1.19209289551e-07
xtol = 1.19209289551e-07
gtol = 1.19209289551e-07
maxfev = None
epsfcn = 1.19209289551e-07
factor = 100.0
verbose = 0
Fit the data¶
The Fit class is used to bundle up the
data, model, statistic, and optimiser choices. The
fit() method runs the optimiser, and
returns a
FitResults instance, which
contains information on how the fit performed. This
infomation includes the
succeeded
attribute, to determine whether the fit converged, as well
as information on the fit (such as the start and end
statistic values) and best-fit parameter values. Note that
the model expression can also be queried for the new
parameter values.
>>> from sherpa.fit import Fit
>>> gfit = Fit(d, g, stat=stat, method=opt)
>>> print(gfit)
data = example
model = gauss1d
stat = LeastSq
method = LevMar
estmethod = Covariance
To actually fit the data, use the
fit() method, which - depending
on the data, model, or statistic being used - can take some
time:
>>> gres = gfit.fit()
>>> print(gres.succeeded)
True
One useful method for interactive analysis is
format(), which returns
a string representation of the fit results, as shown below:
>>> print(gres.format())
Method = levmar
Statistic = leastsq
Initial fit statistic = 180.71
Final fit statistic = 8.06975 at function evaluation 30
Data points = 200
Degrees of freedom = 197
Change in statistic = 172.641
gauss1d.fwhm 1.91572 +/- 0.165982
gauss1d.pos 1.2743 +/- 0.0704859
gauss1d.ampl 3.04706 +/- 0.228618
Note
The LevMar optimiser calculates the
covariance matrix at the best-fit location, and the errors from
this are reported in the output from the call to the
fit() method. In this particular case -
which uses the LeastSq statistic -
the error estimates do not have much meaning. As discussed
below, Sherpa can make use of error estimates
on the data
to calculate meaningful parameter errors.
The sherpa.plot.FitPlot class will display the data
and model. The prepare() method
requires data and model plot objects; in this case the previous
versions can be re-used, although the model plot needs to be
updated to reflect the changes to the model parameters:
>>> from sherpa.plot import FitPlot
>>> fplot = FitPlot()
>>> mplot.prepare(d, g)
>>> fplot.prepare(dplot, mplot)
>>> fplot.plot()
As the model can be evaluated directly, this plot can also be created manually:
>>> plt.plot(d.x, d.y, 'ko', label='Data')
>>> plt.plot(d.x, g(d.x), linewidth=2, label='Gaussian')
>>> plt.legend(loc=2);
Extract the parameter values¶
The fit results include a large number of attributes, many of which
are not relevant here (as the fit was done with no error values).
The following relation is used to convert from the full-width
half-maximum value, used by the Gauss1D
model, to the Gaussian sigma value used to create the data:
\(\rm{FWHM} = 2 \sqrt{2ln(2)} \sigma\):
>>> print(gres)
datasets = None
itermethodname = none
methodname = levmar
statname = leastsq
succeeded = True
parnames = ('gauss1d.fwhm', 'gauss1d.pos', 'gauss1d.ampl')
parvals = (1.915724111406394, 1.2743015983545247, 3.0470560360944017)
statval = 8.069746329529591
istatval = 180.71034547759984
dstatval = 172.64059914807027
numpoints = 200
dof = 197
qval = None
rstat = None
message = successful termination
nfev = 30
>>> conv = 2 * np.sqrt(2 * np.log(2))
>>> ans = dict(zip(gres.parnames, gres.parvals))
>>> print("Position = {:.2f} truth= {:.2f}".format(ans['gauss1d.pos'], pos_true))
Position = 1.27 truth= 1.30
>>> print("Amplitude= {:.2f} truth= {:.2f}".format(ans['gauss1d.ampl'], ampl_true))
Amplitude= 3.05 truth= 3.00
>>> print("Sigma = {:.2f} truth= {:.2f}".format(ans['gauss1d.fwhm']/conv, sigma_true))
Sigma = 0.81 truth= 0.80
The model, and its parameter values, can also be queried directly, as they have been changed by the fit:
>>> print(g)
gauss1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
gauss1d.fwhm thawed 1.91572 1.17549e-38 3.40282e+38
gauss1d.pos thawed 1.2743 -3.40282e+38 3.40282e+38
gauss1d.ampl thawed 3.04706 -3.40282e+38 3.40282e+38
>>> print(g.pos)
val = 1.2743015983545247
min = -3.4028234663852886e+38
max = 3.4028234663852886e+38
units =
frozen = False
link = None
default_val = 0.0
default_min = -3.4028234663852886e+38
default_max = 3.4028234663852886e+38
Including errors¶
For this example, the error on each bin is assumed to be the same, and equal to the true error:
>>> dy = np.ones(x.size) * err_true
>>> de = Data1D('with-errors', x, y, staterror=dy)
>>> print(de)
name = with-errors
x = Float64[200]
y = Float64[200]
staterror = Float64[200]
syserror = None
The statistic is changed from least squares to
chi-square (Chi2), to take advantage
of this extra knowledge (i.e. the Chi-square statistic includes
the error value per bin when calculating the statistic value):
>>> from sherpa.stats import Chi2
>>> ustat = Chi2()
>>> ge = Gauss1D('gerr')
>>> gefit = Fit(de, ge, stat=ustat, method=opt)
>>> geres = gefit.fit()
>>> print(geres.format())
Method = levmar
Statistic = chi2
Initial fit statistic = 4517.76
Final fit statistic = 201.744 at function evaluation 30
Data points = 200
Degrees of freedom = 197
Probability [Q-value] = 0.393342
Reduced statistic = 1.02408
Change in statistic = 4316.01
gerr.fwhm 1.91572 +/- 0.0331963
gerr.pos 1.2743 +/- 0.0140972
gerr.ampl 3.04706 +/- 0.0457235
>>> if not geres.succeeded: print(geres.message)
Since the error value is independent of bin, then the fit results
should be the same here (that is, the parameters in g are the
same as ge):
>>> print(g)
gauss1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
gauss1d.fwhm thawed 1.91572 1.17549e-38 3.40282e+38
gauss1d.pos thawed 1.2743 -3.40282e+38 3.40282e+38
gauss1d.ampl thawed 3.04706 -3.40282e+38 3.40282e+38
>>> print(ge)
gerr
Param Type Value Min Max Units
----- ---- ----- --- --- -----
gerr.fwhm thawed 1.91572 1.17549e-38 3.40282e+38
gerr.pos thawed 1.2743 -3.40282e+38 3.40282e+38
gerr.ampl thawed 3.04706 -3.40282e+38 3.40282e+38
The difference is that more of the fields
in the result structure are populated: in particular the
rstat and
qval fields, which give the
reduced statistic and the probability of obtaining this statistic value
respectively.:
>>> print(geres)
datasets = None
itermethodname = none
methodname = levmar
statname = chi2
succeeded = True
parnames = ('gerr.fwhm', 'gerr.pos', 'gerr.ampl')
parvals = (1.9157241114064163, 1.2743015983545292, 3.047056036094392)
statval = 201.74365823823976
istatval = 4517.758636940002
dstatval = 4316.014978701763
numpoints = 200
dof = 197
qval = 0.3933424667915623
rstat = 1.0240794834428415
message = successful termination
nfev = 30
Error analysis¶
The default error estimation routine is
Covariance, which will be replaced by
Confidence for this example:
>>> from sherpa.estmethods import Confidence
>>> gefit.estmethod = Confidence()
>>> print(gefit.estmethod)
name = confidence
sigma = 1
eps = 0.01
maxiters = 200
soft_limits = False
remin = 0.01
fast = False
parallel = True
numcores = 4
maxfits = 5
max_rstat = 3
tol = 0.2
verbose = False
openinterval = False
Running the error analysis can take time, for particularly complex
models. The default behavior is to use all the available CPU cores
on the machine, but this can be changed with the numcores
attribute. Note that a message is displayed to the screen when each
bound is calculated, to indicate progress:
>>> errors = gefit.est_errors()
gerr.fwhm lower bound: -0.0326327
gerr.fwhm upper bound: 0.0332578
gerr.pos lower bound: -0.0140981
gerr.pos upper bound: 0.0140981
gerr.ampl lower bound: -0.0456119
gerr.ampl upper bound: 0.0456119
The results can be displayed:
>>> print(errors.format())
Confidence Method = confidence
Iterative Fit Method = None
Fitting Method = levmar
Statistic = chi2
confidence 1-sigma (68.2689%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
gerr.fwhm 1.91572 -0.0326327 0.0332578
gerr.pos 1.2743 -0.0140981 0.0140981
gerr.ampl 3.04706 -0.0456119 0.0456119
The ErrorEstResults instance returned by
est_errors() contains the parameter
values and limits:
>>> print(errors)
datasets = None
methodname = confidence
iterfitname = none
fitname = levmar
statname = chi2
sigma = 1
percent = 68.26894921370858
parnames = ('gerr.fwhm', 'gerr.pos', 'gerr.ampl')
parvals = (1.9157241114064163, 1.2743015983545292, 3.047056036094392)
parmins = (-0.0326327431233302, -0.014098074065578947, -0.045611913713536456)
parmaxes = (0.033257800216357714, 0.014098074065578947, 0.045611913713536456)
nfits = 29
The data can be accessed, e.g. to create a dictionary where the keys are the parameter names and the values represent the parameter ranges:
>>> dvals = zip(errors.parnames, errors.parvals, errors.parmins,
... errors.parmaxes)
>>> pvals = {d[0]: {'val': d[1], 'min': d[2], 'max': d[3]}
for d in dvals}
>>> pvals['gerr.pos']
{'min': -0.014098074065578947, 'max': 0.014098074065578947, 'val': 1.2743015983545292}
Screen output¶
The default behavior - when not using the default
Covariance method - is for
est_errors() to print out the parameter
bounds as it finds them, which can be useful in an interactive session
since the error analysis can be slow. This can be controlled using
the Sherpa logging interface.
A single parameter¶
It is possible to investigate the error surface of a single
parameter using the
IntervalProjection class. The following shows
how the error surface changes with the position of the gaussian. The
prepare() method are given
the range over which to vary the parameter (the range is chosen to
be close to the three-sigma limit from the confidence analysis above,
ahd the dotted line is added to indicate the three-sigma
limit above the best-fit for a single parameter):
>>> from sherpa.plot import IntervalProjection
>>> iproj = IntervalProjection()
>>> iproj.prepare(min=1.23, max=1.32, nloop=41)
>>> iproj.calc(gefit, ge.pos)
This can take some time, depending on the complexity of the model and number of steps requested. The resulting data looks like:
>>> iproj.plot()
>>> plt.axhline(geres.statval + 9, linestyle='dotted');
The curve is stored in the
IntervalProjection object (in fact, these
values are created by the call to
calc() and so can be accesed without
needing to create the plot):
>>> print(iproj)
x = [ 1.23 , 1.2323, 1.2345, 1.2368, 1.239 , 1.2412, 1.2435, 1.2457, 1.248 ,
1.2503, 1.2525, 1.2548, 1.257 , 1.2592, 1.2615, 1.2637, 1.266 , 1.2683,
1.2705, 1.2728, 1.275 , 1.2772, 1.2795, 1.2817, 1.284 , 1.2863, 1.2885,
1.2908, 1.293 , 1.2953, 1.2975, 1.2997, 1.302 , 1.3043, 1.3065, 1.3088,
1.311 , 1.3133, 1.3155, 1.3177, 1.32 ]
y = [ 211.597 , 210.6231, 209.6997, 208.8267, 208.0044, 207.2325, 206.5113,
205.8408, 205.2209, 204.6518, 204.1334, 203.6658, 203.249 , 202.883 ,
202.5679, 202.3037, 202.0903, 201.9279, 201.8164, 201.7558, 201.7461,
201.7874, 201.8796, 202.0228, 202.2169, 202.462 , 202.758 , 203.105 ,
203.5028, 203.9516, 204.4513, 205.0018, 205.6032, 206.2555, 206.9585,
207.7124, 208.5169, 209.3723, 210.2783, 211.235 , 212.2423]
min = 1.23
max = 1.32
nloop = 41
delv = None
fac = 1
log = False
A contour plot of two parameters¶
The RegionProjection class supports
the comparison of two parameters. The contours indicate the one,
two, and three sigma contours.
>>> from sherpa.plot import RegionProjection
>>> rproj = RegionProjection()
>>> rproj.prepare(min=[2.8, 1.75], max=[3.3, 2.1], nloop=[21, 21])
>>> rproj.calc(gefit, ge.ampl, ge.fwhm)
As with the interval projection, this step can take time.
>>> rproj.contour()
As with the single-parameter case, the statistic values for the grid are
stored in the RegionProjection object by the
calc() call,
and so can be accesed without needing to create the contour plot. Useful
fields include x0 and x1 (the two parameter values),
y (the statistic value), and levels (the values used for the
contours):
>>> lvls = rproj.levels
>>> print(lvls)
[ 204.03940717 207.92373254 213.57281632]
>>> nx, ny = rproj.nloop
>>> x0, x1, y = rproj.x0, rproj.x1, rproj.y
>>> x0.resize(ny, nx)
>>> x1.resize(ny, nx)
>>> y.resize(ny, nx)
>>> plt.imshow(y, origin='lower', cmap='viridis_r', aspect='auto',
... extent=(x0.min(), x0.max(), x1.min(), x1.max()))
>>> plt.colorbar()
>>> plt.xlabel(rproj.xlabel)
>>> plt.ylabel(rproj.ylabel)
>>> cs = plt.contour(x0, x1, y, levels=lvls)
>>> lbls = [(v, r"${}\sigma$".format(i+1)) for i, v in enumerate(lvls)]
>>> plt.clabel(cs, lvls, fmt=dict(lbls));
Fitting two-dimensional data¶
Sherpa has support for two-dimensional data - that is data defined
on the independent axes x0 and x1. In the example below a
contiguous grid is used, that is the pixel size is constant, but
there is no requirement that this is the case.
>>> np.random.seed(0)
>>> x1, x0 = np.mgrid[:128, :128]
>>> y = 2 * x0**2 - 0.5 * x1**2 + 1.5 * x0 * x1 - 1
>>> y += np.random.normal(0, 0.1, y.shape) * 50000
Creating a data object¶
To support irregularly-gridded data, the multi-dimensional
data classes require
that the coordinate arrays and data values are one-dimensional.
For example, the following code creates a
Data2D object:
>>> from sherpa.data import Data2D
>>> x0axis = x0.ravel()
>>> x1axis = x1.ravel()
>>> yaxis = y.ravel()
>>> d2 = Data2D('img', x0axis, x1axis, yaxis, shape=(128, 128))
>>> print(d2)
name = img
x0 = Int64[16384]
x1 = Int64[16384]
y = Float64[16384]
shape = (128, 128)
staterror = None
syserror = None
Define the model¶
Creating the model is the same as the one-dimensional case; in this
case the Polynom2D class is used
to create a low-order polynomial:
>>> from sherpa.models.basic import Polynom2D
>>> p2 = Polynom2D('p2')
>>> print(p2)
p2
Param Type Value Min Max Units
----- ---- ----- --- --- -----
p2.c thawed 1 -3.40282e+38 3.40282e+38
p2.cy1 thawed 0 -3.40282e+38 3.40282e+38
p2.cy2 thawed 0 -3.40282e+38 3.40282e+38
p2.cx1 thawed 0 -3.40282e+38 3.40282e+38
p2.cx1y1 thawed 0 -3.40282e+38 3.40282e+38
p2.cx1y2 thawed 0 -3.40282e+38 3.40282e+38
p2.cx2 thawed 0 -3.40282e+38 3.40282e+38
p2.cx2y1 thawed 0 -3.40282e+38 3.40282e+38
p2.cx2y2 thawed 0 -3.40282e+38 3.40282e+38
Control the parameters being fit¶
To reduce the number of parameters being fit, the frozen attribute
can be set:
>>> for n in ['cx1', 'cy1', 'cx2y1', 'cx1y2', 'cx2y2']:
...: getattr(p2, n).frozen = True
...:
>>> print(p2)
p2
Param Type Value Min Max Units
----- ---- ----- --- --- -----
p2.c thawed 1 -3.40282e+38 3.40282e+38
p2.cy1 frozen 0 -3.40282e+38 3.40282e+38
p2.cy2 thawed 0 -3.40282e+38 3.40282e+38
p2.cx1 frozen 0 -3.40282e+38 3.40282e+38
p2.cx1y1 thawed 0 -3.40282e+38 3.40282e+38
p2.cx1y2 frozen 0 -3.40282e+38 3.40282e+38
p2.cx2 thawed 0 -3.40282e+38 3.40282e+38
p2.cx2y1 frozen 0 -3.40282e+38 3.40282e+38
p2.cx2y2 frozen 0 -3.40282e+38 3.40282e+38
Fit the data¶
Fitting is no different (the same statistic and optimisation objects used earlier could have been re-used here):
>>> f2 = Fit(d2, p2, stat=LeastSq(), method=LevMar())
>>> res2 = f2.fit()
>>> if not res2.succeeded: print(res2.message)
>>> print(res2)
datasets = None
itermethodname = none
methodname = levmar
statname = leastsq
succeeded = True
parnames = ('p2.c', 'p2.cy2', 'p2.cx1y1', 'p2.cx2')
parvals = (-80.28947555488139, -0.48174521913599017, 1.5022711710872119, 1.9894112623568638)
statval = 400658883390.66907
istatval = 6571471882611.967
dstatval = 6170812999221.298
numpoints = 16384
dof = 16380
qval = None
rstat = None
message = successful termination
nfev = 45
>>> print(p2)
p2
Param Type Value Min Max Units
----- ---- ----- --- --- -----
p2.c thawed -80.2895 -3.40282e+38 3.40282e+38
p2.cy1 frozen 0 -3.40282e+38 3.40282e+38
p2.cy2 thawed -0.481745 -3.40282e+38 3.40282e+38
p2.cx1 frozen 0 -3.40282e+38 3.40282e+38
p2.cx1y1 thawed 1.50227 -3.40282e+38 3.40282e+38
p2.cx1y2 frozen 0 -3.40282e+38 3.40282e+38
p2.cx2 thawed 1.98941 -3.40282e+38 3.40282e+38
p2.cx2y1 frozen 0 -3.40282e+38 3.40282e+38
p2.cx2y2 frozen 0 -3.40282e+38 3.40282e+38
Display the model¶
The model can be visualized by evaluating it over a grid of points and then displaying it:
>>> m2 = p2(x0axis, x1axis).reshape(128, 128)
>>> def pimg(d, title):
... plt.imshow(d, origin='lower', interpolation='nearest',
... vmin=-1e4, vmax=5e4, cmap='viridis')
... plt.axis('off')
... plt.colorbar(orientation='horizontal',
... ticks=[0, 20000, 40000])
... plt.title(title)
...
>>> plt.figure(figsize=(8, 3))
>>> plt.subplot(1, 3, 1);
>>> pimg(y, "Data")
>>> plt.subplot(1, 3, 2)
>>> pimg(m2, "Model")
>>> plt.subplot(1, 3, 3)
>>> pimg(y - m2, "Residual")
Note
The sherpa.image model provides support for interactive
image visualization, but this only works if the
DS9 image viewer is installed.
For the examples in this document, matplotlib plots will be
created to view the data directly.
Simultaneous fits¶
Sherpa allows multiple data sets to be fit at the same time, although there is only really a benefit if there is some model component or value that is shared between the data sets). In this example we have a dataset containing a lorentzian signal with a background component, and another with just the background component. Fitting both together can improve the constraints on the parameter values.
First we start by simulating the data, where the
Polynom1D
class is used to model the background as a straight line, and
Lorentz1D
for the signal:
>>> from sherpa.models import Polynom1D
>>> from sherpa.astro.models import Lorentz1D
>>> tpoly = Polynom1D()
>>> tlor = Lorentz1D()
>>> tpoly.c0 = 50
>>> tpoly.c1 = 1e-2
>>> tlor.pos = 4400
>>> tlor.fwhm = 200
>>> tlor.ampl = 1e4
>>> x1 = np.linspace(4200, 4600, 21)
>>> y1 = tlor(x1) + tpoly(x1) + np.random.normal(scale=5, size=x1.size)
>>> x2 = np.linspace(4100, 4900, 11)
>>> y2 = tpoly(x2) + np.random.normal(scale=5, size=x2.size)
>>> print("x1 size {} x2 size {}".format(x1.size, x2.size))
x1 size 21 x2 size 11
There is no requirement that the data sets have a common grid, as can be seen in a raw view of the data:
>>> plt.plot(x1, y1)
>>> plt.plot(x2, y2)
The fits are set up as before; a data object is needed for each data set, and model instances are created:
>>> d1 = Data1D('a', x1, y1)
>>> d2 = Data1D('b', x2, y2)
>>> fpoly, flor = Polynom1D(), Lorentz1D()
>>> fpoly.c1.thaw()
>>> flor.pos = 4500
To help the fit, we use a simple algorithm to estimate the starting point for the source amplitude, by evaluating the model on the data grid and calculating the change in the amplitude needed to make it match the data:
>>> flor.ampl = y1.sum() / flor(x1).sum()
For simultaneous fits the same optimisation and statistic needs to be used for each fit (this is an area we are looking to improve):
>>> from sherpa.optmethods import NelderMead
>>> stat, opt = LeastSq(), NelderMead()
Set up the fits to the individual data sets:
>>> f1 = Fit(d1, fpoly + flor, stat, opt)
>>> f2 = Fit(d2, fpoly, stat, opt)
and a simultaneous (i.e. to both data sets) fit:
>>> from sherpa.data import DataSimulFit
>>> from sherpa.models import SimulFitModel
>>> sdata = DataSimulFit('all', (d1, d2))
>>> smodel = SimulFitModel('all', (fpoly + flor, fpoly))
>>> sfit = Fit(sdata, smodel, stat, opt)
Note that there is a simulfit() method that
can be used to fit using multiple sherpa.fit.Fit objects,
which wraps the above (using individual fit objects allows some
of the data to be fit first, which may help reduce the parameter
space needed to be searched):
>>> res = sfit.fit()
>>> print(res)
datasets = None
itermethodname = none
methodname = neldermead
statname = leastsq
succeeded = True
parnames = ('polynom1d.c0', 'polynom1d.c1', 'lorentz1d.fwhm', 'lorentz1d.pos', 'lorentz1d.ampl')
parvals = (36.829217311393585, 0.012540257025027028, 249.55651534213359, 4402.7031194359088, 12793.559398547319)
statval = 329.6525419378109
istatval = 3813284.1822045334
dstatval = 3812954.52966
numpoints = 32
dof = 27
qval = None
rstat = None
message = Optimization terminated successfully
nfev = 1152
The values of the numpoints and dof fields show that both datasets
have been used in the fit.
The data can then be viewed (in this case a separate grid is used, but the data objects could be used to define the grid):
>>> plt.plot(x1, y1, label='Data 1')
>>> plt.plot(x2, y2, label='Data 2')
>>> x = np.arange(4000, 5000, 10)
>>> plt.plot(x, (fpoly + flor)(x), linestyle='dotted', label='Fit 1')
>>> plt.plot(x, fpoly(x), linestyle='dotted', label='Fit 2')
>>> plt.legend();
How do you do error analysis? Well, can call sfit.est_errors(), but
that will fail with the current statistic (LeastSq), so need to
change it. The error is 5, per bin, which has to be set up:
>>> print(sfit.calc_stat_info())
name =
ids = None
bkg_ids = None
statname = leastsq
statval = 329.6525419378109
numpoints = 32
dof = 27
qval = None
rstat = None
>>> d1.staterror = np.ones(x1.size) * 5
>>> d2.staterror = np.ones(x2.size) * 5
>>> sfit.stat = Chi2()
>>> check = sfit.fit()
How much did the fit change?:
>>> check.dstatval
0.0
Note that since the error on each bin is the same value, the best-fit
value is not going to be different to the LeastSq result (so dstatval
should be 0):
>>> print(sfit.calc_stat_info())
name =
ids = None
bkg_ids = None
statname = chi2
statval = 13.186101677512438
numpoints = 32
dof = 27
qval = 0.988009259609
rstat = 0.48837413620416437
>>> sres = sfit.est_errors()
>>> print(sres)
datasets = None
methodname = covariance
iterfitname = none
fitname = neldermead
statname = chi2
sigma = 1
percent = 68.2689492137
parnames = ('polynom1d.c0', 'polynom1d.c1', 'lorentz1d.fwhm', 'lorentz1d.pos', 'lorentz1d.ampl')
parvals = (36.829217311393585, 0.012540257025027028, 249.55651534213359, 4402.7031194359088, 12793.559398547319)
parmins = (-4.9568824809960628, -0.0011007470586726147, -6.6079122387075824, -2.0094070026087474, -337.50275154547768)
parmaxes = (4.9568824809960628, 0.0011007470586726147, 6.6079122387075824, 2.0094070026087474, 337.50275154547768)
nfits = 0
Error estimates on a single parameter are as above:
>>> iproj = IntervalProjection()
>>> iproj.prepare(min=6000, max=18000, nloop=101)
>>> iproj.calc(sfit, flor.ampl)
>>> iproj.plot()
Sherpa and CIAO¶
The Sherpa package was developed by the Chandra X-ray Center (CXC) as a general purpose fitting and modeling tool, with specializations for handling X-ray Astronomy data. It is provided as part of the CIAO analysis package, where the code is the same as that available from the Sherpa GitHub page, with the following modifications:
- the plotting and I/O backends use the CIAO versions, namely ChIPS and Crates, rather than matplotlib and astropy;
- a set of customized IPython routines are provided as part of CIAO that automatically load Sherpa and related packages, as well as tweak the IPython look and feel;
- and the CIAO version of Sherpa includes the optional XSPEC model
library (
sherpa.astro.xspec).
The online documentation provided for Sherpa as part of CIAO,
namely http://cxc.harvard.edu/sherpa/, can be used with the
standalone version of Sherpa, but note that the focus of this
documentation is the
session-based API
provided by the
sherpa.astro.ui and sherpa.ui modules.
These are wrappers around the Object-Oriented
interface described in this document, and data management
and utility routines.
What data is to be fit?¶
The Sherpa Data class is used to
carry around the data to be fit: this includes the independent
axis (or axes), the dependent axis (the data), and any
necessary metadata. Although the design of Sherpa supports
multiple-dimensional data sets, the current classes only
support one- and two-dimensional data sets.
Overview¶
The following modules are assumed to have been imported:
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from sherpa.stats import LeastSq
>>> from sherpa.optmethods import LevMar
>>> from sherpa import data
Names¶
The first argument to any of the Data classes is the name
of the data set. This is used for display purposes only,
and can be useful to identify which data set is in use.
It is stored in the name attribute of the object, and
can be changed at any time.
The independent axis¶
The independent axis - or axes - of a data set define the
grid over which the model is to be evaluated. It is referred
to as x, x0, x1, … depending on the dimensionality
of the data (for
binned datasets there are lo
and hi variants).
Although dense multi-dimensional data sets can be stored as arrays with dimensionality greater than one, the internal representation used by Sherpa is often a flattened - i.e. one-dimensional - version.
The dependent axis¶
This refers to the data being fit, and is referred to as y.
Unbinned data¶
Unbinned data sets - defined by classes which do not end in
the name Int - represent point values; that is, the the data
value is the value at the coordinates given by the independent
axis.
Examples of unbinned data classes are
Data1D and Data2D.
>>> np.random.seed(0)
>>> x = np.arange(20, 40, 0.5)
>>> y = x**2 + np.random.normal(0, 10, size=x.size)
>>> d1 = data.Data1D('test', x, y)
>>> print(d1)
name = test
x = Float64[40]
y = Float64[40]
staterror = None
syserror = None
Binned data¶
Binned data sets represent values that are defined over a range,
such as a histogram.
The integrated model classes end in Int: examples are
Data1DInt
and Data2DInt.
It can be a useful optimisation to treat a binned data set as an unbinned one, since it avoids having to estimate the integral of the model over each bin. It depends in part on how the bin size compares to the scale over which the model changes.
>>> z = np.random.gamma(20, scale=0.5, size=1000)
>>> (y, edges) = np.histogram(z)
>>> d2 = data.Data1DInt('gamma', edges[:-1], edges[1:], y)
>>> print(d2)
name = gamma
xlo = Float64[10]
xhi = Float64[10]
y = Int64[10]
staterror = None
syserror = None
>>> plt.bar(d2.xlo, d2.y, d2.xhi - d2.xlo);
Errors¶
There is support for both statistical and systematic
errors by either using the staterror and syserror
parameters when creating the data object, or by changing the
staterror and
syserror attributes of the object.
Filtering data¶
Sherpa supports filtering data sets; that is, temporarily removing parts of the data (perhaps because there are problems, or to help restrict parameter values).
The mask attribute indicates
whether a filter has been applied: if it returns True then
no filter is set, otherwise it is a bool array
where False values indicate those elements that are to be
ignored. The ignore() and
notice() methods are used to
define the ranges to exclude or include. For example, the following
hides those values where the independent axis values are between
21.2 and 22.8:
>>> d1.ignore(21.2, 22.8)
>>> d1.x[np.invert(d1.mask)]
array([ 21.5, 22. , 22.5])
After this, a fit to the data will ignore these values, as shown
below, where the number of degrees of freedom of the first fit,
which uses the filtered data, is three less than the fit to the
full data set (the call to
notice() removes the filter since
no arguments were given):
>>> from sherpa.models import Polynom1D
>>> from sherpa.fit import Fit
>>> mdl = Polynom1D()
>>> mdl.c2.thaw()
>>> fit = Fit(d1, mdl, stat=LeastSq(), method=LevMar())
>>> res1 = fit.fit()
>>> d1.notice()
>>> res2 = fit.fit()
>>> print("Degrees of freedom: {} vs {}".format(res1.dof, res2.dof))
Degrees of freedom: 35 vs 38
Visualizing a data set¶
The data objects contain several methods which can be used to
visualize the data, but do not provide any direct plotting
or display capabilities. There are low-level routines which
provide access to the data - these include the
to_plot() and
to_contour() methods - but the
preferred approach is to use the classes defined in the
sherpa.plot module, which are described in the
visualization section:
>>> from sherpa.plot import DataPlot
>>> pdata = DataPlot()
>>> pdata.prepare(d2)
>>> pdata.plot()
Although the data represented by d2 is
a histogram, the values are displayed at the center of the bin.
The plot objects automatically handle any filters applied to the data, as shown below.
>>> d1.ignore(25, 30)
>>> d1.notice(26, 27)
>>> pdata.prepare(d1)
>>> pdata.plot()
Note
The plot object stores the data given in the
prepare() call,
so that changes to the underlying objects will not be reflected
in future calls to
plot()
unless a new call to
prepare() is made.
>>> d1.notice()
At this point, a call to pdata.plot() would re-create the previous
plot, even though the filter has been removed from the underlying
data object.
Evaluating a model¶
The eval_model() and
eval_model_to_fit()
methods can be used
to evaluate a model on the grid defined by the data set. The
first version uses the full grid, whereas the second respects
any filtering applied to the data.
>>> d1.notice(22, 25)
>>> y1 = d1.eval_model(mdl)
>>> y2 = d1.eval_model_to_fit(mdl)
>>> x2 = d1.x[d1.mask]
>>> plt.plot(d1.x, d1.y, 'ko', label='Data')
>>> plt.plot(d1.x, y1, '--', label='Model (all points)')
>>> plt.plot(x2, y2, linewidth=2, label='Model (filtered)')
>>> plt.legend(loc=2)
Reference/API¶
The sherpa.data module¶
Tools for creating, storing, inspecting, and manipulating data sets
Classes
Data(name, indep, y[, staterror, syserror]) |
Data class for generic, N-Dimensional data sets, where N depends on the number of independent axes passed during initialization. |
Data1D(name, x, y[, staterror, syserror]) |
|
Data1DAsymmetricErrs(name, x, y, elo, ehi[, …]) |
1-D data set with asymmetric errors Note: elo and ehi shall be stored as delta values from y |
Data1DInt(name, xlo, xhi, y[, staterror, …]) |
1-D integrated data set |
Data2D(name, x0, x1, y[, shape, staterror, …]) |
|
Data2DInt(name, x0lo, x1lo, x0hi, x1hi, y[, …]) |
2-D integrated data set |
DataSimulFit(name, datasets) |
Store multiple data sets. |
BaseData |
Base class for all data classes. |
DataSpace1D(filter, x) |
Class for representing 1-D Data Space. |
DataSpace2D(filter, x0, x1) |
Class for representing 2-D Data Spaces. |
DataSpaceND(filter, indep) |
Class for representing arbitray N-Dimensional data domains |
Filter() |
A class for representing filters of N-Dimentional datasets. |
IntegratedDataSpace1D(filter, xlo, xhi) |
Same as DataSpace1D, but for supporting integrated data sets. |
IntegratedDataSpace2D(filter, x0lo, x1lo, …) |
Same as DataSpace2D, but for supporting integrated data sets. |
Class Inheritance Diagram¶
The sherpa.astro.data module¶
Classes for storing, inspecting, and manipulating astronomical data sets
Classes
DataPHA(name, channel, counts[, staterror, …]) |
PHA data set, including any associated instrument and background data. |
DataARF(name, energ_lo, energ_hi, specresp) |
ARF data set. |
DataRMF(name, detchans, energ_lo, energ_hi, …) |
RMF data set. |
DataIMG(name, x0, x1, y[, shape, staterror, …]) |
Image data set, including functions for coordinate transformations |
DataIMGInt(name, x0lo, x1lo, x0hi, x1hi, y) |
|
DataRosatRMF(name, detchans, energ_lo, …) |
Class Inheritance Diagram¶
Creating model instances¶
The sherpa.models and sherpa.astro.models
namespaces provides a collection of one- and
two-dimensional models as a convenience; the actual definition of
each particular model depends on its type.
The following modules are assumed to have been imported for this section:
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from sherpa import models
Creating a model instance¶
Models must be created before there parameter values can
be set. In this case a one-dimensional gaussian using the
Gauss1D class:
>>> g = models.Gauss1D()
>>> print(g)
gauss1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
gauss1d.fwhm thawed 10 1.17549e-38 3.40282e+38
gauss1d.pos thawed 0 -3.40282e+38 3.40282e+38
gauss1d.ampl thawed 1 -3.40282e+38 3.40282e+38
A description of the model is provided by help(g).
The parameter values have a current value, a valid range
(as given by the the minimum and maximum columns in the table above),
and a units field. The units field is a string, describing the
expected units for the parameter; there is currently no support for
using astropy.units to set a
parameter value. The “Type” column refers to whether the parameter is
fixed, (frozen) or can be varied during a fit (thawed),
as described below, in the Freezing and Thawing parameters section.
Models can be given a name, to help distinguish multiple versions of the same model type. The default value is the lower-case version of the class name.
>>> g.name
'gauss1d'
>>> h = models.Gauss1D('other')
>>> print(h)
other
Param Type Value Min Max Units
----- ---- ----- --- --- -----
other.fwhm thawed 10 1.17549e-38 3.40282e+38
other.pos thawed 0 -3.40282e+38 3.40282e+38
other.ampl thawed 1 -3.40282e+38 3.40282e+38
>>> h.name
'other'
The model classes are expected to derive from the
ArithmeticModel class.
Combining models¶
Models can be combined and shared by using the standard Python
numerical operators. For instance, a one-dimensional gaussian
plus a flat background - using the
Const1D class - would be
represented by the following model:
>>> src1 = models.Gauss1D('src1')
>>> back = models.Const1D('back')
>>> mdl1 = src1 + back
>>> print(mdl1)
(src1 + back)
Param Type Value Min Max Units
----- ---- ----- --- --- -----
src1.fwhm thawed 10 1.17549e-38 3.40282e+38
src1.pos thawed 0 -3.40282e+38 3.40282e+38
src1.ampl thawed 1 -3.40282e+38 3.40282e+38
back.c0 thawed 1 -3.40282e+38 3.40282e+38
Now consider fitting a second dataset where it is known that the background is two times higher than the first:
>>> src2 = models.Gauss1D('src2')
>>> mdl2 = src2 + 2 * back
>>> print(mdl2)
(src2 + (2 * back))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
src2.fwhm thawed 10 1.17549e-38 3.40282e+38
src2.pos thawed 0 -3.40282e+38 3.40282e+38
src2.ampl thawed 1 -3.40282e+38 3.40282e+38
back.c0 thawed 1 -3.40282e+38 3.40282e+38
The two models can then be fit separately or simultaneously. In this
example the two source models (the Gaussian component) were completely
separate, but they could have been identical - in which case
mdl2 = src1 + 2 * back would have been used instead - or
parameter linking could be used to constrain the
models. An example of the use of linking would be to force the two
FWHM (full-width half-maximum)
parameters to be the same but to let the position and amplitude
values vary independently.
More information is available in the combining models documentation.
Changing a parameter¶
The parameters of a model - those numeric variables that control the
shape of the model, and that can be varied during a fit -
can be accesed as attributes, both to read or change
the current settings. The
val attribute
contains the current value:
>>> print(h.fwhm)
val = 10.0
min = 1.17549435082e-38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 10.0
default_min = 1.17549435082e-38
default_max = 3.40282346639e+38
>>> h.fwhm.val
10.0
>>> h.fwhm.min
1.1754943508222875e-38
>>> h.fwhm.val = 15
>>> print(h.fwhm)
val = 15.0
min = 1.17549435082e-38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 15.0
default_min = 1.17549435082e-38
default_max = 3.40282346639e+38
Assigning a value to a parameter directly (i.e. without using the
val attribute) also works:
>>> h.fwhm = 12
>>> print(h.fwhm)
val = 12.0
min = 1.17549435082e-38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 12.0
default_min = 1.17549435082e-38
default_max = 3.40282346639e+38
The soft and hard limits of a parameter¶
Each parameter has two sets of limits, which are referred to as
“soft” and “hard”. The soft limits are shown when the model
is displayed, and refer to the
min
and
max
attributes for the parameter, whereas the hard limits are
given by the
hard_min
and
hard_max
(which are not displayed, and can not be changed).
>>> print(h)
other
Param Type Value Min Max Units
----- ---- ----- --- --- -----
other.fwhm thawed 12 1.17549e-38 3.40282e+38
other.pos thawed 0 -3.40282e+38 3.40282e+38
other.ampl thawed 1 -3.40282e+38 3.40282e+38
>>> print(h.fwhm)
val = 12.0
min = 1.17549435082e-38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 12.0
default_min = 1.17549435082e-38
default_max = 3.40282346639e+38
These limits act to bound the acceptable parameter range; this is often because certain values are physically impossible, such as having a negative value for the full-width-half-maxium value of a Gaussian, but can also be used to ensure that the fit is restricted to a meaningful part of the search space. The hard limits are set by the model class, and represent the full valid range of the parameter, whereas the soft limits can be changed by the user, although they often default to the same values as the hard limits.
Setting a parameter to a value outside its soft limits will
raise a ParameterErr exception.
During a fit the paramater values are bound by the soft limits, and a screen message will be displayed if an attempt to move outside this range was made. During error analysis the parameter values are allowed outside the soft limits, as long as they remain inside the hard limits.
Guessing a parameter’s value from the data¶
Sherpa models have a
guess()
method which is used to seed the paramters (or
parameter) with values and
soft-limit ranges
which match the data.
The idea is to move the parameters to values appropriate
for the data, which can avoid un-needed computation by
the optimiser.
The existing guess routines are very basic - such as
picking the index of the largest value in the data for
the peak location - and do not always account for the
full complexity of the model expression, so care should
be taken when using this functionality.
The arguments depend on the model type, since both the
independent and dependent axes may be used, but the
to_guess() method of
a data object will return the correct data (assuming the
dimensionality and type match):
>>> mdl.guess(*data.to_guess())
Note that the soft limits can be changed, as in this example which ensures the position of the gaussian falls within the grid of points (since this is the common situation; if the source is meant to lie outside the data range then the limits will need to be increased manually):
>>> yg, xg = np.mgrid[4000:4050:10, 3000:3070:10]
>>> r2 = (xg - 3024.2)**2 + (yg - 4011.7)**2
>>> zg = 2400 * np.exp(-r2 / 1978.2)
>>> d2d = Data2D('example', xg.flatten(), yg.flatten(), zg.flatten(),
shape=zg.shape)
>>> mdl = Gauss2D('mdl')
>>> print(mdl)
mdl
Param Type Value Min Max Units
----- ---- ----- --- --- -----
mdl.fwhm thawed 10 1.17549e-38 3.40282e+38
mdl.xpos thawed 0 -3.40282e+38 3.40282e+38
mdl.ypos thawed 0 -3.40282e+38 3.40282e+38
mdl.ellip frozen 0 0 0.999
mdl.theta frozen 0 -6.28319 6.28319 radians
mdl.ampl thawed 1 -3.40282e+38 3.40282e+38
>>> mdl.guess(*d2d.to_guess())
>>> print(mdl)
mdl
Param Type Value Min Max Units
----- ---- ----- --- --- -----
mdl.fwhm thawed 10 1.17549e-38 3.40282e+38
mdl.xpos thawed 3020 3000 3060
mdl.ypos thawed 4010 4000 4040
mdl.ellip frozen 0 0 0.999
mdl.theta frozen 0 -6.28319 6.28319 radians
mdl.ampl thawed 2375.22 2.37522 2.37522e+06
Freezing and Thawing parameters¶
Not all model parameters should be varied during a fit: perhaps
the data quality is not sufficient to constrain all the parameters,
it is already known, the parameter is highly correlated with
another, or perhaps the parameter value controls a behavior of the
model that should not vary during a fit (such as the interpolation
scheme to use). The frozen
attribute controls whether a fit
should vary that parameter or not; it can be changed directly,
as shown below:
>>> h.fwhm.frozen
False
>>> h.fwhm.frozen = True
or via the freeze()
and thaw()
methods for the parameter.
>>> h.fwhm.thaw()
>>> h.fwhm.frozen
False
There are times when a model parameter should never be varied
during a fit. In this case the
alwaysfrozen
attribute will be set to True (this particular
parameter is read-only).
Linking parameters¶
There are times when it is useful for one parameter to be related to another: this can be equality, such as saying that the width of two model components are the same, or a functional form, such as saying that the position of one component is a certain distance away from another component. This concept is refererred to as linking parameter values. The second case incudes the first - where the functional relationship is equality - but it is treated separately here as it is a common operation. Lnking parameters also reduces the number of free parameters in a fit.
The following examples use the same two model components:
>>> g1 = models.Gauss1D('g1')
>>> g2 = models.Gauss1D('g2')
Linking parameter values requires referring to the parameter, rather
than via the val attribute.
The link attribute
is set to the link value (and is None for parameters that are
not linked).
Equality¶
After the following, the two gaussian components have the same width:
>>> g2.fwhm.val
10.0
>>> g2.fwhm = g1.fwhm
>>> g1.fwhm = 1024
>>> g2.fwhm.val
1024.0
>>> g1.fwhm.link is None
True
>>> g2.fwhm.link
<Parameter 'fwhm' of model 'g1'>
When displaying the model, the value and link expression are included:
>>> print(g2)
g2
Param Type Value Min Max Units
----- ---- ----- --- --- -----
g2.fwhm linked 1024 expr: g1.fwhm
g2.pos thawed 0 -3.40282e+38 3.40282e+38
g2.ampl thawed 1 -3.40282e+38 3.40282e+38
Functional relationship¶
The link can accept anything that evaluates to a value, such as adding a constant.
>>> g2.pos = g1.pos + 8234
>>> g1.pos = 1200
>>> g2.pos.val
9434.0
The CompositeParameter class
controls how parameters are combined. In this case the result
is a BinaryOpParameter object.
Including another parameter¶
It is possible to include other parameters in a link expression,
which can lead to further constraints on the fit. For instance,
rather than using a fixed separation, a range can be used. One
way to do this is to use a Const1D
model, restricting the value its one parameter can vary.
>>> sep = models.Const1D('sep')
>>> print(sep)
sep
Param Type Value Min Max Units
----- ---- ----- --- --- -----
sep.c0 thawed 1 -3.40282e+38 3.40282e+38
>>> g2.fwhm = g1.fwhm + sep.c0
>>> sep.c0 = 1200
>>> sep.c0.min = 800
>>> sep.c0.max = 1600
In this example, the separation of the two components is restricted to lie in the range 800 to 1600.
In order for the optimiser to recognize that it needs to vary the
new parameter (sep.c0), the component must be included in the
model expression. As it does not contribute to the model output
directly, it should be multiplied by zero. So, for this example
the model to be fit would be given by an expression like:
>>> mdl = g1 + g2 + 0 * sep
Resetting parameter values¶
The
reset()
method of a parameter will change the parameter settings (which
includes the status of the thawed flag and allowed ranges,
as well as the value) to the values they had the last time
the parameter was explicitly set. That is, it does not restore
the initial values used when the model was created, but the
last values the user set.
The model class has its own
reset()
method which calls reset on the thawed parameters. This can be used to
change the starting point of a fit
to see how robust the optimiser is by:
- explicitly setting parameter values (or using the default values)
- fit the data
- call reset
- change one or more parameters
- refit
Inspecting models and parameters¶
Models, whether a single component or composite, contain a
pars attribute which is a tuple of all the parameters
for that model. This can be used to programatically query
or change the parameter values.
There are several attributes that return arrays of values
for the thawed parameters of the model expression: the most
useful is thawedpars,
which gives the current values.
Composite models can be queried to find the individual
components using the parts attribute, which contains
a tuple of the components (these components can themselves
be composite objects).
Evaluating a model¶
Binned and Unbinned grids¶
Sherpa supports models for both
unbinned and
binned data
sets. The output of a model depends on how it is called
(is it sent just the grid points or the bin edges), how
the integrate flag of
the model component is set, and whether the model supports
both or just one case.
The Const1D model represents
a constant value, which means that for an unbinned
dataset the model evaluates to a single value (the
c0 parameter):
>>> from sherpa.models.basic import Const1D
>>> mdl = Const1D()
>>> mdl.c0 = 0.1
>>> mdl([1, 2, 3])
array([ 0.1, 0.1, 0.1])
>>> mdl([-4000, 12000])
array([ 0.1, 0.1])
The default value for its
integrate flag is
True:
>>> mdl.integrate
True
which means that this value is multiplied by the bin width when given a binned grid (i.e. when sent in the low and high edges of each bin):
>>> mdl([10, 20, 30], [20, 30, 50])
array([ 1., 1., 2.])
When the integrate flag is unset, the model no longer
multiplies by the bin width, and so acts similarly to the
unbinned case:
>>> mdl.integrate = False
>>> mdl([10, 20, 30], [20, 30, 50])
array([ 0.1, 0.1, 0.1])
The behavior in all these three cases depends on the model - for instance some models may raise an exception, ignore the high-edge values in the binned case, or use the mid-point - and so the model documentation should be reviewed.
The following example uses the
Polynom1D class to model the
linear relation
\(y = mx + c\) with the origin at \(x = 1400\),
an offset of 2, and a gradient of 1:
>>> from sherpa.models.basic import Polynom1D
>>> poly = Polynom1D()
>>> poly.offset = 1400
>>> poly.c0 = 2
>>> poly.c1 = 1
>>> x = [1391, 1396, 1401, 1406, 1411]
>>> poly(x)
array([ -7., -2., 3., 8., 13.])
As the integrate flag is set, the model is integrated across each bin:
>>> poly.integrate
True
>>> xlo, xhi = x[:-1], x[1:]
>>> y = poly(xlo, xhi)
>>> y
array([-22.5, 2.5, 27.5, 52.5])
Thanks to the easy functonal form chosen for this example, it is easy to confirm that these are the values of the integrated model:
>>> (y[:-1] + y[1:]) * 5 / 2.0
array([-22.5, 2.5, 27.5, 52.5])
Turning off the integrate flag for this model shows that it
uses the low-edge of the bin when evaluating the model:
>>> poly.integrate = False
>>> poly(xlo, xhi)
array([-7., -2., 3., 8.])
Combining models and parameters¶
Most of the examples show far have used a single model component, such as a one-dimensional polynomial or a two-dimensional gaussian, but individual components can be combined together by addition, multiplication, subtraction, or even division. Components can also be combined with scalar values or - with great care - NumPy vectors. Parameter values can be “combined” by linking them together using mathematical expressions. The case of one model requiring the results of another model is discussed in the convolution section.
Note
There is currently no restriction on combining models of different types. This means that there is no exception raised when combining a one-dimensional model with a two-dimensional one. It is only when the model is evaluated that an error may be raised.
Model Expressions¶
A model, whether it is required to create a
sherpa.fit.Fit object or the argument to
the sherpa.ui.set_source() call, is expected to
behace like an instance of the
sherpa.models.model.ArithmeticModel class.
Instances can be combined as
numeric types
since the class defines methods for addition, subtraction,
multiplication, division, modulus, and exponentiation.
This means that Sherpa model instances can be combined with
other Python terms, such as the weighted combination of
model components cpt1, cpt2, and cpt3:
cpt1 * (cpt2 + 0.8 * cpt3)
Since the models are evaluated on a grid, it is possible to include a NumPy vector in the expression, but this is only possible in restricted situations, when the grid size is known (i.e. the model expression is not going to be used in a general setting).
Example¶
The following example fits two one-dimensional gaussians to a simulated dataset. It is based on the AstroPy modelling documentation, but has linked the positions of the two gaussians during the fit.
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from sherpa import data, models, stats, fit, plot
Since the example uses many different parts of the Sherpa API, the various modules are imported directly, rather than their contents, to make it easier to work out what each symbol refers to.
Note
Some Sherpa modules re-export symbols from other modules, which
means that a symbol can be found in several modules. An example
is sherpa.models.basic.Gauss1D, which can also be
imported as sherpa.models.Gauss1D.
Creating the simulated data¶
To provide a repeatable example, the NumPy random number generator is set to a fixed value:
>>> np.random.seed(42)
The two components used to create the simulated dataset are called
sim1 and sim2:
>>> s1 = models.Gauss1D('sim1')
>>> s2 = models.Gauss1D('sim2')
The individual components can be displayed, as the __str__
method of the model class creates a display which includes the
model expression and then a list of the paramters:
>>> print(s1)
sim1
Param Type Value Min Max Units
----- ---- ----- --- --- -----
sim1.fwhm thawed 10 1.17549e-38 3.40282e+38
sim1.pos thawed 0 -3.40282e+38 3.40282e+38
sim1.ampl thawed 1 -3.40282e+38 3.40282e+38
The pars attribute contains
a tuple of all the parameters in a model instance. This can be
queried to find the attributes of the parameters (each element
of the tuple is a Parameter
object):
>>> [p.name for p in s1.pars]
['fwhm', 'pos', 'ampl']
These components can be combined using standard mathematical operations; for example addition:
>>> sim_model = s1 + s2
The sim_model object represents the sum of two gaussians, and
contains both the input models (using different names when creating
model components - so here sim1 and sim2 - can make it
easier to follow the logic of more-complicated model combinations):
>>> print(sim_model)
(sim1 + sim2)
Param Type Value Min Max Units
----- ---- ----- --- --- -----
sim1.fwhm thawed 10 1.17549e-38 3.40282e+38
sim1.pos thawed 0 -3.40282e+38 3.40282e+38
sim1.ampl thawed 1 -3.40282e+38 3.40282e+38
sim2.fwhm thawed 10 1.17549e-38 3.40282e+38
sim2.pos thawed 0 -3.40282e+38 3.40282e+38
sim2.ampl thawed 1 -3.40282e+38 3.40282e+38
The pars attribute now includes parameters from both components,
and so
the fullname
attribute is used to discriminate between the two components:
>>> [p.fullname for p in sim_model.pars]
['sim1.fwhm', 'sim1.pos', 'sim1.ampl', 'sim2.fwhm', 'sim2.pos', 'sim2.ampl']
Since the original models are still accessible, they can be used to
change the parameters of the combined model. The following sets the
first component (sim1) to be centered at x = 0 and the
second one at x = 0.5:
>>> s1.ampl = 1.0
>>> s1.pos = 0.0
>>> s1.fwhm = 0.5
>>> s2.ampl = 2.5
>>> s2.pos = 0.5
>>> s2.fwhm = 0.25
The model is evaluated on the grid, and “noise” added to it (using a normal distribution centered on 0 with a standard deviation of 0.2):
>>> x = np.linspace(-1, 1, 200)
>>> y = sim_model(x) + np.random.normal(0., 0.2, x.shape)
These arrays are placed into a Sherpa data object, using the
Data1D class, since it will be fit
below, and then a plot created to show the simulated data:
>>> d = data.Data1D('multiple', x, y)
>>> dplot = plot.DataPlot()
>>> dplot.prepare(d)
>>> dplot.plot()
What is the composite model?¶
The result of the combination is a
BinaryOpModel, which has
op,
lhs,
and rhs
attributes which describe the structure of the combination:
>>> sim_model
<BinaryOpModel model instance '(sim1 + sim2)'>
>>> sim_model.op
<ufunc 'add'>
>>> sim_model.lhs
<Gauss1D model instance 'sim1'>
>>> sim_model.rhs
<Gauss1D model instance 'sim2'>
There is also a
parts attribute
which contains all the elements of the model (in this case the
combination of the lhs and rhs attributes):
>>> sim_model.parts
(<Gauss1D model instance 'sim1'>, <Gauss1D model instance 'sim2'>)
>>> for cpt in sim_model.parts:
...: print(cpt)
sim1
Param Type Value Min Max Units
----- ---- ----- --- --- -----
sim1.fwhm thawed 0.5 1.17549e-38 3.40282e+38
sim1.pos thawed 0 -3.40282e+38 3.40282e+38
sim1.ampl thawed 1 -3.40282e+38 3.40282e+38
sim2
Param Type Value Min Max Units
----- ---- ----- --- --- -----
sim2.fwhm thawed 0.25 1.17549e-38 3.40282e+38
sim2.pos thawed 0.5 -3.40282e+38 3.40282e+38
sim2.ampl thawed 2.5 -3.40282e+38 3.40282e+38
As the BinaryOpModel class is a subclass of the
ArithmeticModel class, the
combined model can be treated as a single model instance; for instance
it can be evaluated on a grid by passing in an array of values:
>>> sim_model([-1.0, 0, 1])
array([ 1.52587891e-05, 1.00003815e+00, 5.34057617e-05])
Setting up the model¶
Rather than use the model components used to simulate the data, new instances are created and combined to create the model:
>>> g1 = models.Gauss1D('g1')
>>> g2 = models.Gauss1D('g2')
>>> mdl = g1 + g2
In this particular fit, the separation of the two models is going
to be assumed to be known, so the two pos parameters can
be linked together, which means that there
is one less free parameter in the fit:
>>> g2.pos = g1.pos + 0.5
The FWHM parameters are changed as the default value of 10 is not appropriate for this data (since the independent axis ranges from -1 to 1):
>>> g1.fwhm = 0.1
>>> g2.fwhm = 0.1
The display of the combined model shows that the g2.pos
parameter is now linked to the g1.pos value:
>>> print(mdl)
(g1 + g2)
Param Type Value Min Max Units
----- ---- ----- --- --- -----
g1.fwhm thawed 0.1 1.17549e-38 3.40282e+38
g1.pos thawed 0 -3.40282e+38 3.40282e+38
g1.ampl thawed 1 -3.40282e+38 3.40282e+38
g2.fwhm thawed 0.1 1.17549e-38 3.40282e+38
g2.pos linked 0.5 expr: (g1.pos + 0.5)
g2.ampl thawed 1 -3.40282e+38 3.40282e+38
Note
It is a good idea to check the parameter ranges - that is their minimum and maximum values - to make sure they are appropriate for the data.
The model is evaluated with its initial parameter values so that it can be compared to the best-fit location later:
>>> ystart = mdl(x)
Fitting the model¶
The initial model can be added to the data plot either directly,
with matplotlib commands, or using the
ModelPlot class to overlay onto the
DataPlot display:
>>> mplot = plot.ModelPlot()
>>> mplot.prepare(d, mdl)
>>> dplot.plot()
>>> mplot.plot(overplot=True)
As can be seen, the initial values for the gaussian positions are close to optimal. This is unlikely to happen in real-world situations!
As there are no errors for the data set, the least-square statistic
(LeastSq) is used (so that
the fit attempts to minimise the separation between the model and
data with no weighting), along with the default optimiser:
>>> f = fit.Fit(d, mdl, stats.LeastSq())
>>> res = f.fit()
>>> res.succeeded
True
When displayig the results, the FitPlot
class is used since it combines both data and model plots (after
updating the mplot object to include the new model parameter
values):
>>> fplot = plot.FitPlot()
>>> mplot.prepare(d, mdl)
>>> fplot.prepare(dplot, mplot)
>>> fplot.plot()
>>> plt.plot(x, ystart, label='Start')
>>> plt.legend(loc=2)
As can be seen below, the position of the g2 gaussian remains
linked to that of g1:
>>> print(mdl)
(g1 + g2)
Param Type Value Min Max Units
----- ---- ----- --- --- -----
g1.fwhm thawed 0.515565 1.17549e-38 3.40282e+38
g1.pos thawed 0.00431538 -3.40282e+38 3.40282e+38
g1.ampl thawed 0.985078 -3.40282e+38 3.40282e+38
g2.fwhm thawed 0.250698 1.17549e-38 3.40282e+38
g2.pos linked 0.504315 expr: (g1.pos + 0.5)
g2.ampl thawed 2.48416 -3.40282e+38 3.40282e+38
Accessing the linked parameter¶
The pars attribute of a model instance provides access to the
individual Parameter objects.
These can be used to query - as shown below - or change the model
values:
>>> for p in mdl.pars:
...: if p.link is None:
...: print("{:10s} -> {:.3f}".format(p.fullname, p.val))
...: else:
...: print("{:10s} -> link to {}".format(p.fullname, p.link.name))
g1.fwhm -> 0.516
g1.pos -> 0.004
g1.ampl -> 0.985
g2.fwhm -> 0.251
g2.pos -> link to (g1.pos + 0.5)
g2.ampl -> 2.484
The linked parameter is actually an instance of the
CompositeParameter
class, which allows parameters to be combined in a similar
manner to models:
>>> g2.pos
<Parameter 'pos' of model 'g2'>
>>> print(g2.pos)
val = 0.504315379302
min = -3.40282346639e+38
max = 3.40282346639e+38
units =
frozen = True
link = (g1.pos + 0.5)
default_val = 0.504315379302
default_min = -3.40282346639e+38
default_max = 3.40282346639e+38
>>> g2.pos.link
<BinaryOpParameter '(g1.pos + 0.5)'>
>>> print(g2.pos.link)
val = 0.504315379302
min = -3.40282346639e+38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 0.504315379302
default_min = -3.40282346639e+38
default_max = 3.40282346639e+38
Convolution¶
A convolution model requires both the evaluation grid and the data to convolve. Examples include using a point-spread function (PSF) model to modify a two-dimensional model to account for blurring due to the instrument (or other sources, such as the atmosphere for ground-based Astronomical data sets), or the redistribution of the counts as a function of energy as modelled by the RMF when analyzing astronomical X-ray spectra.
Two-dimensional convolution with a PSF¶
The sherpa.astro.instrument.PSFModel class augments the
behavior of
sherpa.instrument.PSFModel by supporting images with
a World Coordinate System (WCS).
Including a PSF in a model expression¶
There are two steps to including a PSF:
- create an instance
- apply the instance to the model components
The “kernel” of the PSF is the actual data used to represent the
blurring, and can be given as a numeric array or as a Sherpa model.
In the following example a simple 3 by 3 array is used to represent
the PSF, but it first has to be converted into a
Data object:
>>> from sherpa.data import Data2D
>>> from sherpa.instrument import PSFModel
>>> k = np.asarray([[0, 1, 0], [1, 0, 1], [0, 1, 0]])
>>> yg, xg = np.mgrid[:3, :3]
>>> kernel = Data2D('kdata', xg.flatten(), yg.flatten(), k.flatten(),
shape=k.shape)
>>> psf = PSFModel(kernel=kernel)
>>> print(psf)
psfmodel
Param Type Value Min Max Units
----- ---- ----- --- --- -----
psfmodel.kernel frozen kdata
psfmodel.radial frozen 0 0 1
psfmodel.norm frozen 1 0 1
As shown below, the data in the
PSF is renormalized so that its sum matches the norm parameter,
which here is set to 1.
The following example sets up a simple model expression which represents
the sum of a single pixel and a line of pixels, using
Box2D for both.
>>> from sherpa.models.basic import Box2D
>>> pt = Box2D('pt')
>>> pt.xlow, pt.xhi = 1.5, 2.5
>>> pt.ylow, pt.yhi = 2.5, 3.5
>>> pt.ampl = 8
>>> box = Box2D('box')
>>> box.xlow, box.xhi = 4, 10
>>> box.ylow, box.yhi = 6.5, 7.5
>>> box.ampl = 10
>>> unconvolved_mdl = pt + box
>>> print(unconvolved_mdl)
(pt + box)
Param Type Value Min Max Units
----- ---- ----- --- --- -----
pt.xlow thawed 1.5 -3.40282e+38 3.40282e+38
pt.xhi thawed 2.5 -3.40282e+38 3.40282e+38
pt.ylow thawed 2.5 -3.40282e+38 3.40282e+38
pt.yhi thawed 3.5 -3.40282e+38 3.40282e+38
pt.ampl thawed 10 -3.40282e+38 3.40282e+38
box.xlow thawed 4 -3.40282e+38 3.40282e+38
box.xhi thawed 10 -3.40282e+38 3.40282e+38
box.ylow thawed 6.5 -3.40282e+38 3.40282e+38
box.yhi thawed 7.5 -3.40282e+38 3.40282e+38
box.ampl thawed 10 -3.40282e+38 3.40282e+38
Note
Although Sherpa provides the Delta2D
class, it is suggested that alternatives such as
Box2D be used instead, since a
delta function is very sensitive to the location at which it
is evaluated. However, including a Box2D component in a fit can still
be problematic since the output of the model does not vary smoothly
as any of the bin edges change, which is a challenge for the
optimisers provided with Sherpa.
Rather than being another term in the model expression - that is, an item that is added, subtracted, multiplied, or divided into an existing expression - the PSF model “wraps” the model it is to convolve. This can be a single model or - as in this case - a composite one:
>>> convolved_mdl = psf(unconvolved_mdl)
>>> print(convolved_mdl)
psfmodel((pt + box))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
pt.xlow thawed 1.5 -3.40282e+38 3.40282e+38
pt.xhi thawed 2.5 -3.40282e+38 3.40282e+38
pt.ylow thawed 2.5 -3.40282e+38 3.40282e+38
pt.yhi thawed 3.5 -3.40282e+38 3.40282e+38
pt.ampl thawed 10 -3.40282e+38 3.40282e+38
box.xlow thawed 4 -3.40282e+38 3.40282e+38
box.xhi thawed 10 -3.40282e+38 3.40282e+38
box.ylow thawed 6.5 -3.40282e+38 3.40282e+38
box.yhi thawed 7.5 -3.40282e+38 3.40282e+38
box.ampl thawed 10 -3.40282e+38 3.40282e+38
This new expression can be treated as any other Sherpa model, which means that we can apply extra terms to it, such as adding a background component that is not affected by the PSF:
>>> from sherpa.models.basic import Const2D
>>> bgnd = Const2D('bgnd')
>>> bgnd.c0 = 0.25
>>> print(convolved_mdl + bgnd)
(psfmodel((pt + box)) + bgnd)
Param Type Value Min Max Units
----- ---- ----- --- --- -----
pt.xlow thawed 1.5 -3.40282e+38 3.40282e+38
pt.xhi thawed 2.5 -3.40282e+38 3.40282e+38
pt.ylow thawed 2.5 -3.40282e+38 3.40282e+38
pt.yhi thawed 3.5 -3.40282e+38 3.40282e+38
pt.ampl thawed 10 -3.40282e+38 3.40282e+38
box.xlow thawed 4 -3.40282e+38 3.40282e+38
box.xhi thawed 10 -3.40282e+38 3.40282e+38
box.ylow thawed 6.5 -3.40282e+38 3.40282e+38
box.yhi thawed 7.5 -3.40282e+38 3.40282e+38
box.ampl thawed 10 -3.40282e+38 3.40282e+38
bgnd.c0 thawed 0.25 -3.40282e+38 3.40282e+38
In the following this extra term (bgnd) is not included to simplify
the comparison between the unconvolved and convolved versions.
Evaluating a model including a PSF¶
The PSF-convolved model can be evaluated - in most cases - just as is done for ordinary models. That is by supplying it with the grid coordinates to use. However, the need to convolve the data with a fixed grid does limit this somewhat.
For this example, a grid covering the points 0 to 9 inclusive is used for each axis (with a unit pixel size), which means that the unconvolved model can be evaluated with the following:
>>> yg, xg = np.mgrid[:10, :10]
>>> xg1d, yg1d = xg.flatten(), yg.flatten()
>>> m1 = unconvolved_mdl(xg1d, yg1d).reshape(xg.shape)
An easier alternative, once the PSF is included, is to create an
empty dataset with the given grid (that is, a dataset for which we
do not care about the dependent axis), and use the
eval_model() method to
evaluate the model (the result for m1 is the same whichever
approach is used):
>>> blank = Data2D('blank', xg1d, yg1d, np.ones(xg1d.shape), xg.shape)
>>> m1 = blank.eval_model(unconvolved_mdl).reshape(xg.shape)
The “point source” is located at x = 2, y = 3 and the line
starts at x=5 and extends to the end of the grid (at y=7).
Note
In this example the image coordinates were chosen to be the same
as those drawn by matplotlib. The extent parameter of the
imshow call can be used when this correspondance does not
hold.
The PSF model includes a
fold() method, which is used to
pre-calculate terms needed for the convolution (which is done using
a fourier transform), and so needs the grid over which it is to be
applied. This is done by passing in a Sherpa dataset, such as the
blank example we have just created:
>>> psf.fold(blank)
>>> m2 = blank.eval_model(convolved_mdl).reshape(xg.shape)
The kernel used redistributes flux from the central pixel to its four
immediate neighbors equally, which is what has happened to the point
source at (2, 2). The result for the line is to blur the line
slightly, but note that the convolution has “wrapped around”, so that
the flux that should have been placed into the pixel at (10, 7),
which is off the grid, has been moved to (0, 7).
Note
If the fold method is not called then evaluating the model will raise the following exception:
PSFErr: PSF model has not been folded
Care must be taken to ensure that fold is called whenever the grid has changed. This suggests that the same PSF model should not be used in simultaneous fits, unless it is known that the grid is the same in the multiple datasets.
The PSF Normalization¶
Since the norm parameter of the PSF model was set to 1, the PSF
convolution is flux preserving, even at the edges thanks to the
wrap-around behavior of the fourier transform. This can be seen by
comparing the signal in the unconvolved and convolved images, which
are (to numerical precision) the same:
>>> m1.sum()
60.0
>>> m2.sum()
60.0
The use of a fourier transform means that low-level signal will be
found in many pixels which would expect to be 0. For example,
looking at the row of pixels at y = 7 gives:
>>> m2[7]
array([2.50000000e+00, 1.73472348e-16, 5.20417043e-16, 4.33680869e-16,
2.50000000e+00, 2.50000000e+00, 5.00000000e+00, 5.00000000e+00,
5.00000000e+00, 2.50000000e+00])
Evaluating the model on a different grid¶
Sherpa now provides experimental support for evaluating a model on a different grid to the independent axis. This can be used to better model the underlying distribution (by use of a finer grid), or to include features that lie outside the data range but, due to the use of a convolution model, can affect the values within the data range.
Existing Sherpa models take advantage of this support by inheriting
from the
RegriddableModel1D
or
RegriddableModel2D classes.
At present the only documentation is provided in the
low-level API module.
Simulating data¶
Simulating a data set normally involves:
- evaluate the model
- add in noise
This may need to be repeated several times for complex models, such as when different components have different noise models or the noise needs to be added before evaluation by a component.
The model evaluation would be performed using the techniques
described in this section, and then the noise term can be
handled with sherpa.utils.poisson_noise() or routines from
NumPy or SciPy to evaluate noise, such as numpy.random.standard_normal.
>>> import numpy as np
>>> from sherpa.models.basic import Polynom1D
>>> np.random.seed(235)
>>> x = np.arange(10, 100, 12)
>>> mdl = Polynom1D('mdl')
>>> mdl.offset = 35
>>> mdl.c1 = 0.5
>>> mdl.c2 = 0.12
>>> ymdl = mdl(x)
>>> from sherpa.utils import poisson_noise
>>> ypoisson = poisson_noise(ymdl)
>>> from numpy.random import standard_normal, normal
>>> yconst = ymdl + standard_normal(ymdl.shape) * 10
>>> ydata = ymdl + normal(scale=np.sqrt(ymdl))
Caching model evaluations¶
Sherpa contains a rudimentary system for cacheing the results
of a model evaluation. It is related to the
modelCacher1d()
function decorator.
Examples¶
The following examples show the different ways that a model can be evaluted, for a range of situations. The direct method is often sufficient, but for more complex cases it can be useful to ask a data object to evaluate the model, particularly if you want to include instrumental responses, such as a RMF and ARF.
Evaluating a one-dimensional model directly¶
In the following example a one-dimensional gaussian is evaluated
on a grid of 5 points by
using the model object directly.
The first approch just calls the model with the evaluation
grid (here the array x),
which uses the parameter values as defined in the model itself:
>>> from sherpa.models.basic import Gauss1D
>>> gmdl = Gauss1D()
>>> gmdl.fwhm = 100
>>> gmdl.pos = 5050
>>> gmdl.ampl = 50
>>> x = [4800, 4900, 5000, 5100, 5200]
>>> y1 = gmdl(x)
The second uses the calc()
method, where the parameter values must be specified in the
call along with the grid on which to evaluate the model.
The order matches that of the parameters in the model, which can be
found from the
pars attribute of the model:
>>> [p.name for p in gmdl.pars]
['fwhm', 'pos', 'ampl']
>>> y2 = gmdl.calc([100, 5050, 100], x)
>>> y2 / y1
array([ 2., 2., 2., 2., 2.])
Since in this case the amplitude (the last parameter value) is twice
that used to create y1 the ratio is 2 for each bin.
Evaluating a 2D model to match a Data2D object¶
In the following example the model is evaluated on a grid
specified by a dataset, in this case a set of two-dimensional
points stored in a Data2D object.
First the data is set up (there are only four points
in this example to make things easy to follow).
>>> from sherpa.data import Data2D
>>> x0 = [1.0, 1.9, 2.4, 1.2]
>>> x1 = [-5.0, -7.0, 2.3, 1.2]
>>> y = [12.1, 3.4, 4.8, 5.2]
>>> twod = Data2D('data', x0, x1, y)
For demonstration purposes, the Box2D
model is used, which represents a rectangle (any points within the
xlow
to
xhi
and
ylow
to
yhi
limits are set to the
ampl
value, those outside are zero).
>>> from sherpa.models.basic import Box2D
>>> mdl = Box2D('mdl')
>>> mdl.xlow = 1.5
>>> mdl.xhi = 2.5
>>> mdl.ylow = -9.0
>>> mdl.yhi = 5.0
>>> mdl.ampl = 10.0
The coverage have been set so that some of the points are within the “box”, and so are set to the amplitude value when the model is evaluated.
>>> twod.eval_model(mdl)
array([ 0., 10., 10., 0.])
The eval_model() method evaluates
the model on the grid defined by the data set, so it is the same
as calling the model directly with these values:
>>> twod.eval_model(mdl) == mdl(x0, x1)
array([ True, True, True, True], dtype=bool)
The eval_model_to_fit() method
will apply any filter associated with the data before
evaluating the model. At this time there is no filter
so it returns the same as above.
>>> twod.eval_model_to_fit(mdl)
array([ 0., 10., 10., 0.])
Adding a simple spatial filter - that excludes one of
the points within the box - with
ignore() now results
in a difference in the outputs of
eval_model()
and
eval_model_to_fit(),
as shown below. The call to
get_indep()
is used to show the grid used by
eval_model_to_fit().
>>> twod.ignore(x0lo=2, x0hi=3, x1l0=0, x1hi=10)
>>> twod.eval_model(mdl)
array([ 0., 10., 10., 0.])
>>> twod.get_indep(filter=True)
(array([ 1. , 1.9, 1.2]), array([-5. , -7. , 1.2]))
>>> twod.eval_model_to_fit(mdl)
array([ 0., 10., 0.])
Evaluating a model using a DataPHA object¶
This example is similar to the
two-dimensional case above,
in that it again shows the differences between the
eval_model()
and
eval_model_to_fit()
methods. The added complication in this
case is that the response information provided with a PHA file
is used to convert between the “native” axis of the
PHA file (channels) and that of the model (energy or
wavelength). This conversion is handled automatically
by the two methods (the
following example
shows how this can be done manually).
To start with, the data is loaded from a file, which also loads in the associated ARF and RMF files:
>>> from sherpa.astro.io import read_pha
>>> pha = read_pha('3c273.pi')
WARNING: systematic errors were not found in file '3c273.pi'
statistical errors were found in file '3c273.pi'
but not used; to use them, re-read with use_errors=True
read ARF file 3c273.arf
read RMF file 3c273.rmf
WARNING: systematic errors were not found in file '3c273_bg.pi'
statistical errors were found in file '3c273_bg.pi'
but not used; to use them, re-read with use_errors=True
read background file 3c273_bg.pi
>>> pha
<DataPHA data set instance '3c273.pi'>
>>> pha.get_arf()
<DataARF data set instance '3c273.arf'>
>>> pha.get_rmf()
<DataRMF data set instance '3c273.rmf'>
The returned object - here pha - is an instance of the
sherpa.astro.data.DataPHA class - which has a number
of attributes and methods specialized to handling PHA data.
This particular file has grouping information in it, that it it contains
GROUPING and QUALITY columns, so Sherpa
applies them: that is, the number of bins over which the data is
analysed is smaller than the number of channels in the file because
each bin can consist of multiple channels. For this file,
there are 46 bins after grouping (the filter argument to the
get_dep() call applies both
filtering and grouping steps, but so far no filter has been applied):
>>> pha.channel.size
1024
>>> pha.get_dep().size
1024
>>> pha.grouped
True
>>> pha.get_dep(filter=True).size
46
A filter - in this case to restrict to only bins that cover the
energy range 0.5 to 7.0 keV - is applied with the
notice() call, which
removes four bins for this particular data set:
>>> pha.set_analysis('energy')
>>> pha.notice(0.5, 7.0)
>>> pha.get_dep(filter=True).size
42
A power-law model (PowLaw1D) is
created and evaluated by the data object:
>>> from sherpa.models.basic import PowLaw1D
>>> mdl = PowLaw1D()
>>> y1 = pha.eval_model(mdl)
>>> y2 = pha.eval_model_to_fit(mdl)
>>> y1.size
1024
>>> y2.size
42
The eval_model() call
evaluates the model over the full dataset and does not
apply any grouping, so it returns a vector with 1024 elements.
In contrast, eval_model_to_fit()
applies both filtering and grouping, and returns a vector that
matches the data (i.e. it has 42 elements).
The filtering and grouping information is dynamic, in that it
can be changed without having to re-load the data set. The
ungroup() call removes
the grouping, but leaves the 0.5 to 7.0 keV energy filter:
>>> pha.ungroup()
>>> y3 = pha.eval_model_to_fit(mdl)
>>> y3.size
644
Evaluating a model using PHA responses¶
The sherpa.astro.data.DataPHA class handles the
response information automatically, but it is possible to
directly apply the response information to a model using
the sherpa.astro.instrument module. In the following
example the
RSPModelNoPHA
and
RSPModelPHA
classes are used to wrap a power-law model
(PowLaw1D)
so that the
instrument responses - the ARF and RMF -
are included in the model evaluation.
>>> from sherpa.astro.io import read_arf, read_rmf
>>> arf = read_arf('3c273.arf')
>>> rmf = read_rmf('3c273.rmf')
>>> rmf.detchans
1024
The number of channels in the RMF - that is, the number of bins over which the RMF is defined - is 1024.
>>> from sherpa.models.basic import PowLaw1D
>>> mdl = PowLaw1D()
The RSPModelNoPHA class
models the inclusion of both the ARF and RMF:
>>> from sherpa.astro.instrument import RSPModelNoPHA
>>> inst = RSPModelNoPHA(arf, rmf, mdl)
>>> inst
<RSPModelNoPHA model instance 'apply_rmf(apply_arf(powlaw1d))'>
>>> print(inst)
apply_rmf(apply_arf(powlaw1d))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
powlaw1d.gamma thawed 1 -10 10
powlaw1d.ref frozen 1 -3.40282e+38 3.40282e+38
powlaw1d.ampl thawed 1 0 3.40282e+38
Note
The RMF and ARF are represented as models that “enclose” the
spectrum - that is, they are written apply_rmf(model) and
apply_arf(model) rather than rmf * model - since they
may perform a convolution or rebinning (ARF) of the model
output.
The return value (inst) behaves as a normal Shepra model, for
example:
>>> from sherpa.models.model import ArithmeticModel
>>> isinstance(inst, ArithmeticModel)
True
>>> inst.pars
(<Parameter 'gamma' of model 'powlaw1d'>,
<Parameter 'ref' of model 'powlaw1d'>,
<Parameter 'ampl' of model 'powlaw1d'>)
The model can therefore be evaluated by calling it
with a grid (as used in the first example
above), except that
the input grid is ignored and the “native” grid of the
response information is used. In this case, no matter the
size of the one-dimensional array passed to inst, the
output has 1024 elements (matching the number of channels in
the RMF):
>>> inst(np.arange(1, 1025))
array([ 0., 0., 0., ..., 0., 0., 0.])
>>> inst([0.1, 0.2, 0.3])
array([ 0., 0., 0., ..., 0., 0., 0.])
>>> inst([0.1, 0.2, 0.3]).size
1024
>>> inst([10, 20]) == inst([])
array([ True, True, True, ..., True, True, True], dtype=bool)
The output of this call represents the number of counts expected in each bin:
>>> chans = np.arange(rmf.offset, rmf.offset + rmf.detchans)
>>> ydet = inst(chans)
>>> plt.plot(chans, ydet)
>>> plt.xlabel('Channel')
>>> plt.ylabel('Count / s')
Note
The interpretation of the model output as being in units of “counts” (or a rate) depends on the normalisation (or amplitude) of the model components, and whether any term representing the exposure time has been included.
XSPEC additive models - such as XSapec -
return values that have units of photon/cm^2/s (that is, the spectrum
is integrated across each bin), which when passed through the
ARF and RMF results in count/s (the ARF has units of cm^2 and the
RMF can be thought of as converting photons to counts).
The Sherpa models, such as PowLaw1D,
do not in general have units (so that the models can be applied
to different data sets). This means that the interpretation of
the normalization or amplitude term depends on how the model
is being used.
The data in the EBOUNDS extension of the RMF - which provides
an approximate mapping from channel to energy for visualization
purposes only - is available as the
e_min
and
e_max
attributes of the
DataRMF object returned by
read_rmf().
The ARF object may contain an
exposure time, in its
exposure
attribute:
>>> print(rmf)
name = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp = UInt64[1090]
f_chan = UInt64[2002]
n_chan = UInt64[2002]
matrix = Float64[61834]
offset = 1
e_min = Float64[1024]
e_max = Float64[1024]
ethresh = 1e-10
>>> print(arf)
name = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo = None
bin_hi = None
exposure = 38564.141454905
ethresh = 1e-10
These can be used to create a plot of energy versus counts per energy bin:
>>> # intersperse the low and high edges of each bin
>>> x = np.vstack((rmf.e_min, rmf.e_max)).T.flatten()
>>> # normalize each bin by its width and include the exposure time
>>> y = arf.exposure * ydet / (rmf.e_max - rmf.e_min)
>>> # Repeat for the low and high edges of each bin
>>> y = y.repeat(2)
>>> plt.plot(x, y, '-')
>>> plt.yscale('log')
>>> plt.ylim(1e3, 1e7)
>>> plt.xlim(0, 10)
>>> plt.xlabel('Energy (keV)')
>>> plt.ylabel('Count / keV')
Note
The bin widths are small enough that it is hard to make out each bin on this plot.
The
RSPModelPHA
class adds in a
DataPHA object, which lets the
evaluation grid be determined by any filter applied to the
data object. In the following, the
read_pha() call reads in a PHA
file, along with its associated ARF and RMF (because the
ANCRFILE and RESPFILE keywords are set in the
header of the PHA file), which means that there is no need
to call
read_arf()
and
read_rmf()
to creating the RSPModelPHA instance.
>>> from sherpa.astro.io import read_pha
>>> from sherpa.astro.instrument import RSPModelPHA
>>> pha = read_pha('3c273.pi')
WARNING: systematic errors were not found in file '3c273.pi'
statistical errors were found in file '3c273.pi'
but not used; to use them, re-read with use_errors=True
read ARF file 3c273.arf
read RMF file 3c273.rmf
WARNING: systematic errors were not found in file '3c273_bg.pi'
statistical errors were found in file '3c273_bg.pi'
but not used; to use them, re-read with use_errors=True
read background file 3c273_bg.pi
>>> arf2 = pha2.get_arf()
>>> rmf2 = pha2.get_rmf()
>>> mdl2 = PowLaw1D('mdl2')
>>> inst2 = RSPModelPHA(arf2, rmf2, pha2, mdl2)
>>> print(inst2)
apply_rmf(apply_arf(mdl2))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
mdl2.gamma thawed 1 -10 10
mdl2.ref frozen 1 -3.40282e+38 3.40282e+38
mdl2.ampl thawed 1 0 3.40282e+38
The model again is evaluated on the channel grid defined by the RMF:
>>> inst2([]).size
1024
The DataPHA object can be
adjusted to select a subset of data. The default is to use
the full channel range:
>>> pha2.set_analysis('energy')
>>> pha2.get_filter()
'0.124829999695:12.410000324249'
>>> pha2.get_filter_expr()
'0.1248-12.4100 Energy (keV)'
This can be changed with the
notice()
and
ignore()
methods:
>>> pha2.notice(0.5, 7.0)
>>> pha2.get_filter()
'0.518300011754:8.219800233841'
>>> pha2.get_filter_expr()
'0.5183-8.2198 Energy (keV)'
Note
Since the channels have a finite width, the method of filtering
(in other words, is it notice or ignore)
determines whether a channel that
includes a boundary (in this case 0.5 and 7.0 keV) is included
or excluded from the final range. The dataset used in this example
includes grouping information, which is automatically applied,
which is why the upper limit of the included range is at 8 rather
than 7 keV:
>>> pha2.grouped
True
Ignore a range within the previous range to make the plot more interesting.
>>> pha2.ignore(2.0, 3.0)
>>> pha2.get_filter_expr()
'0.5183-1.9199,3.2339-8.2198 Energy (keV)'
When evaluate, over whole 1-1024 channels, but can take advantage
of the filter if within a pair of calls to
startup()
and
teardown()
(this is performed
automatically by certain routines, such as within a fit):
>>> y1 = inst2([])
>>> inst2.startup()
>>> y2 = inst2([])
>>> inst2.teardown()
>>> y1.size, y2.size
(1024, 1024)
>>> np.all(y1 == y2)
False
>>> plt.plot(pha2.channel, y1, label='all')
>>> plt.plot(pha2.channel, y2, label='filtered')
>>> plt.xscale('log')
>>> plt.yscale('log')
>>> plt.ylim(0.001, 1)
>>> plt.xlim(5, 1000)
>>> plt.legend(loc='center')
Why is the exposure time not being included?
Or maybe this?¶
This could come first, although maybe need a separate section on how to use astro.instruments (since this is geeting quite long now).
>>> from sherpa.astro.io import read_pha
>>> from sherpa.models.basic import PowLaw1D
>>> pha = read_pha('3c273.pi')
>>> pl = PowLaw1D()
>>> from sherpa.astro.instrument import Response1D, RSPModelPHA
>>> rsp = Response1D(pha)
>>> mdl = rsp(pl)
>>> isinstance(mdl, RSPModelPHA)
>>> print(mdl)
apply_rmf(apply_arf((38564.608926889 * powlaw1d)))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
powlaw1d.gamma thawed 1 -10 10
powlaw1d.ref frozen 1 -3.40282e+38 3.40282e+38
powlaw1d.ampl thawed 1 0 3.40282e+38
Note that the exposure time - taken from the PHA or the ARF - is included so that the normalization is correct.
Direct evaluation of the model¶
Normally Sherpa will handle model evaluation automatically, such as
during a fit or displaying the model results. However, the
models can be evalutated directly by passing in the
grid
(the independent axis)
directly. If mdl is an instance of a Sherpa model - that
is it is derived from the
Model
class - then there are two standard ways to perform this
evaluation:
Call the model with the grid directly - e.g. for a one-dimensional grid use one of:
mdl(x) mdl(xlo, xhi)
Use the
calc()method, which requires a sequence of parameter values and then the grid; for the one-dimensional case this would be:mdl.calc(pars, x) mdl.calc(pars, xlo, xhi)
In this case the parameter values do not need to match the values stored in the model itself. This can be useful when a model is to be embedded within another one, as shown in the two-dimensional user model example.
Evaluating a model with a data object¶
It is also possible to pass a model to a data object
and evaluate the model on a grid appropriate for the data,
using the
eval_model() and
eval_model_to_fit() methods.
This can be useful when working in an environment where the mapping
between the “native” grids used to represent data and models is
not a simple one-to-one relation, such as when analyzing
astronomical X-ray spectral data with an associated response
(i.e. a RMF file), or
when using a PSF.
Available Models¶
Note
The models in sherpa.astro.xspec are only available if
Sherpa was built with support for the
XSPEC model library.
This section describes the classes that implement models used to describe and fit data, while the Reference/API section below describes the classes used to create these models.
Writing your own model¶
A model class can be created to fit any function, or interface with external code.
A one-dimensional model¶
An example is the AstroPy trapezoidal model, which has four parameters: the amplitude of the central region, the center and width of this region, and the slope. The following model class, which was not written for efficiancy or robustness, implements this interface:
import numpy as np
from sherpa.models import model
__all__ = ('Trap1D', )
def _trap1d(pars, x):
"""Evaluate the Trapezoid.
Parameters
----------
pars: sequence of 4 numbers
The order is amplitude, center, width, and slope.
These numbers are assumed to be valid (e.g. width
is 0 or greater).
x: sequence of numbers
The grid on which to evaluate the model. It is expected
to be a floating-point type.
Returns
-------
y: sequence of numbers
The model evaluated on the input grid.
Notes
-----
This is based on the interface described at
http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Trapezoid1D.html
but implemented without looking at the code, so any errors
are not due to AstroPy.
"""
(amplitude, center, width, slope) = pars
# There are five segments:
# xlo = center - width/2
# xhi = center + width/2
# x0 = xlo - amplitude/slope
# x1 = xhi + amplitude/slope
#
# flat xlo <= x < xhi
# slope x0 <= x < xlo
# xhi <= x < x1
# zero x < x0
# x >= x1
#
hwidth = width / 2.0
dx = amplitude / slope
xlo = center - hwidth
xhi = center + hwidth
x0 = xlo - dx
x1 = xhi + dx
out = np.zeros(x.size)
out[(x >= xlo) & (x < xhi)] = amplitude
idx = np.where((x >= x0) & (x < xlo))
out[idx] = slope * x[idx] - slope * x0
idx = np.where((x >= xhi) & (x < x1))
out[idx] = - slope * x[idx] + slope * x1
return out
class Trap1D(model.RegriddableModel1D):
"""A one-dimensional trapezoid.
The model parameters are:
ampl
The amplitude of the central (flat) segment (zero or greater).
center
The center of the central segment.
width
The width of the central segment (zero or greater).
slope
The gradient of the slopes (zero or greater).
"""
def __init__(self, name='trap1d'):
self.ampl = model.Parameter(name, 'ampl', 1, min=0, hard_min=0)
self.center = model.Parameter(name, 'center', 1)
self.width = model.Parameter(name, 'width', 1, min=0, hard_min=0)
self.slope = model.Parameter(name, 'slope', 1, min=0, hard_min=0)
model.RegriddableModel1D.__init__(self, name,
(self.ampl, self.center, self.width,
self.slope))
def calc(self, pars, x, *args, **kwargs):
"""Evaluate the model"""
# If given an integrated data set, use the center of the bin
if len(args) == 1:
x = (x + args[0]) / 2
return _trap1d(pars, x)
This can be used in the same manner as the
Gauss1D model
in the quick guide to Sherpa.
First, create the data to fit:
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> np.random.seed(0)
>>> x = np.linspace(-5., 5., 200)
>>> ampl_true = 3
>>> pos_true = 1.3
>>> sigma_true = 0.8
>>> err_true = 0.2
>>> y = ampl_true * np.exp(-0.5 * (x - pos_true)**2 / sigma_true**2)
>>> y += np.random.normal(0., err_true, x.shape)
Now create a Sherpa data object:
>>> from sherpa.data import Data1D
>>> d = Data1D('example', x, y)
Set up the user model:
>>> from trap import Trap1D
>>> t = Trap1D()
>>> print(t)
trap1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
trap1d.ampl thawed 1 0 3.40282e+38
trap1d.center thawed 1 -3.40282e+38 3.40282e+38
trap1d.width thawed 1 0 3.40282e+38
trap1d.slope thawed 1 0 3.40282e+38
Finally, perform the fit:
>>> from sherpa.fit import Fit
>>> from sherpa.stats import LeastSq
>>> from sherpa.optmethods import LevMar
>>> tfit = Fit(d, t, stat=LeastSq(), method=LevMar())
>>> tres = tfit.fit()
>>> if not tres.succeeded: print(tres.message)
Rather than use a ModelPlot object,
the overplot argument can be set to allow multiple values
in the same plot:
>>> from sherpa import plot
>>> dplot = plot.DataPlot()
>>> dplot.prepare(d)
>>> dplot.plot()
>>> mplot = plot.ModelPlot()
>>> mplot.prepare(d, t)
>>> mplot.plot(overplot=True)
A two-dimensional model¶
The two-dimensional case is similar to the one-dimensional case,
with the major difference being the number of independent axes to
deal with. In the following example the model is assumed to only be
applied to non-integrated data sets, as it simplifies the implementation
of the calc method.
It also shows one way of embedding models from a different system, in this case the two-dimemensional polynomial model from the AstroPy package.
from sherpa.models import model
from astropy.modeling.polynomial import Polynomial2D
__all__ = ('WrapPoly2D', )
class WrapPoly2D(model.RegriddableModel2D):
"""A two-dimensional polynomial from AstroPy, restricted to degree=2.
The model parameters (with the same meaning as the underlying
AstroPy model) are:
c0_0
c1_0
c2_0
c0_1
c0_2
c1_1
"""
def __init__(self, name='wrappoly2d'):
self._actual = Polynomial2D(degree=2)
self.c0_0 = model.Parameter(name, 'c0_0', 0)
self.c1_0 = model.Parameter(name, 'c1_0', 0)
self.c2_0 = model.Parameter(name, 'c2_0', 0)
self.c0_1 = model.Parameter(name, 'c0_1', 0)
self.c0_2 = model.Parameter(name, 'c0_2', 0)
self.c1_1 = model.Parameter(name, 'c1_1', 0)
model.RegriddableModel2D.__init__(self, name,
(self.c0_0, self.c1_0, self.c2_0,
self.c0_1, self.c0_2, self.c1_1))
def calc(self, pars, x0, x1, *args, **kwargs):
"""Evaluate the model"""
# This does not support 2D integrated data sets
mdl = self._actual
for n in ['c0_0', 'c1_0', 'c2_0', 'c0_1', 'c0_2', 'c1_1']:
pval = getattr(self, n).val
getattr(mdl, n).value = pval
return mdl(x0, x1)
Repeating the 2D fit by first setting up the data to fit:
>>> np.random.seed(0)
>>> y2, x2 = np.mgrid[:128, :128]
>>> z = 2. * x2 ** 2 - 0.5 * y2 ** 2 + 1.5 * x2 * y2 - 1.
>>> z += np.random.normal(0., 0.1, z.shape) * 50000.
Put this data into a Sherpa data object:
>>> from sherpa.data import Data2D
>>> x0axis = x2.ravel()
>>> x1axis = y2.ravel()
>>> d2 = Data2D('img', x0axis, x1axis, z.ravel(), shape=(128,128))
Create an instance of the user model:
>>> from poly import WrapPoly2D
>>> wp2 = WrapPoly2D('wp2')
>>> wp2.c1_0.frozen = True
>>> wp2.c0_1.frozen = True
Finally, perform the fit:
>>> f2 = Fit(d2, wp2, stat=LeastSq(), method=LevMar())
>>> res2 = f2.fit()
>>> if not res2.succeeded: print(res2.message)
>>> print(res2)
datasets = None
itermethodname = none
methodname = levmar
statname = leastsq
succeeded = True
parnames = ('wp2.c0_0', 'wp2.c2_0', 'wp2.c0_2', 'wp2.c1_1')
parvals = (-80.289475553599914, 1.9894112623565667, -0.4817452191363118, 1.5022711710873158)
statval = 400658883390.6685
istatval = 6571934382318.328
dstatval = 6.17127549893e+12
numpoints = 16384
dof = 16380
qval = None
rstat = None
message = successful termination
nfev = 80
>>> print(wp2)
wp2
Param Type Value Min Max Units
----- ---- ----- --- --- -----
wp2.c0_0 thawed -80.2895 -3.40282e+38 3.40282e+38
wp2.c1_0 frozen 0 -3.40282e+38 3.40282e+38
wp2.c2_0 thawed 1.98941 -3.40282e+38 3.40282e+38
wp2.c0_1 frozen 0 -3.40282e+38 3.40282e+38
wp2.c0_2 thawed -0.481745 -3.40282e+38 3.40282e+38
wp2.c1_1 thawed 1.50227 -3.40282e+38 3.40282e+38
The sherpa.models.basic module¶
Classes
Box1D([name]) |
One-dimensional box function. |
Const1D([name]) |
A constant model for one-dimensional data. |
Cos([name]) |
One-dimensional cosine function. |
Delta1D([name]) |
One-dimensional delta function. |
Erf([name]) |
One-dimensional error function. |
Erfc([name]) |
One-dimensional complementary error function. |
Exp([name]) |
One-dimensional exponential function. |
Exp10([name]) |
One-dimensional exponential function, base 10. |
Gauss1D([name]) |
One-dimensional gaussian function. |
Integrate1D([name]) |
|
Log([name]) |
One-dimensional natural logarithm function. |
Log10([name]) |
One-dimensional logarithm function, base 10. |
LogParabola([name]) |
One-dimensional log-parabolic function. |
NormGauss1D([name]) |
One-dimensional normalised gaussian function. |
Poisson([name]) |
One-dimensional Poisson function. |
Polynom1D([name]) |
One-dimensional polynomial function of order 8. |
PowLaw1D([name]) |
One-dimensional power-law function. |
Scale1D([name]) |
A constant model for one-dimensional data. |
Sin([name]) |
One-dimensional sine function. |
Sqrt([name]) |
One-dimensional square root function. |
StepHi1D([name]) |
One-dimensional step function. |
StepLo1D([name]) |
One-dimensional step function. |
TableModel([name]) |
|
Tan([name]) |
One-dimensional tan function. |
UserModel([name, pars]) |
Support for user-supplied models. |
Box2D([name]) |
Two-dimensional box function. |
Const2D([name]) |
A constant model for two-dimensional data. |
Delta2D([name]) |
Two-dimensional delta function. |
Gauss2D([name]) |
Two-dimensional gaussian function. |
Scale2D([name]) |
A constant model for two-dimensional data. |
SigmaGauss2D([name]) |
Two-dimensional gaussian function (varying sigma). |
NormGauss2D([name]) |
Two-dimensional normalised gaussian function. |
Polynom2D([name]) |
Two-dimensional polynomial function. |
Class Inheritance Diagram¶
The sherpa.astro.models module¶
Classes
Atten([name]) |
Model the attenuation by the Inter-Stellar Medium (ISM). |
BBody([name]) |
A one-dimensional Blackbody model. |
BBodyFreq([name]) |
A one-dimensional Blackbody model (frequency). |
BPL1D([name]) |
One-dimensional broken power-law function. |
Beta1D([name]) |
One-dimensional beta model function. |
Beta2D([name]) |
Two-dimensional beta model function. |
DeVaucouleurs2D([name]) |
Two-dimensional de Vaucouleurs model. |
Dered([name]) |
A de-reddening model. |
Disk2D([name]) |
Two-dimensional uniform disk model. |
Edge([name]) |
Photoabsorption edge model. |
HubbleReynolds([name]) |
Two-dimensional Hubble-Reynolds model. |
JDPileup([name]) |
A CCD pileup model for the ACIS detectors on Chandra. |
LineBroad([name]) |
A one-dimensional line-broadening profile. |
Lorentz1D([name]) |
One-dimensional normalized Lorentz model function. |
Lorentz2D([name]) |
Two-dimensional un-normalised Lorentz function. |
MultiResponseSumModel(name[, pars]) |
|
NormBeta1D([name]) |
One-dimensional normalized beta model function. |
Schechter([name]) |
One-dimensional Schecter model function. |
Sersic2D([name]) |
Two-dimensional Sersic model. |
Shell2D([name]) |
A homogenous spherical 3D shell projected onto 2D. |
Class Inheritance Diagram¶
The sherpa.astro.optical module¶
Optical models intended for SED Analysis
The models match those used by the SpecView application [1], and are intended for un-binned one-dimensional data sets defined on a wavelength grid, with units of Angstroms. When used with a binned data set the lower-edge of each bin is used to evaluate the model. This module does not contain all the spectral components from SpecView ([2]).
References
| [1] | http://www.stsci.edu/institute/software_hardware/specview/ |
| [2] | http://specview.stsci.edu/javahelp/Components.html |
Classes
AbsorptionEdge([name]) |
Optical model of an absorption edge. |
AbsorptionGaussian([name]) |
Gaussian function for modeling absorption (equivalent width). |
AbsorptionLorentz([name]) |
Lorentz function for modeling absorption (equivalent width). |
AbsorptionVoigt([name]) |
Voigt function for modeling absorption (equivalent width). |
AccretionDisk([name]) |
A model of emission due to an accretion disk. |
BlackBody([name]) |
Emission from a black body as a function of wavelength. |
Bremsstrahlung([name]) |
Bremsstrahlung emission. |
BrokenPowerlaw([name]) |
Broken power-law model. |
CCM([name]) |
Galactic extinction: the Cardelli, Clayton, and Mathis model. |
EmissionGaussian([name]) |
Gaussian function for modeling emission. |
EmissionLorentz([name]) |
Lorentz function for modeling emission. |
EmissionVoigt([name]) |
Voigt function for modeling emission. |
FM([name]) |
UV extinction curve: Fitzpatrick and Massa 1988. |
LMC([name]) |
LMC extinction: the Howarth model. |
LogAbsorption([name]) |
Gaussian function for modeling absorption (log of fwhm). |
LogEmission([name]) |
Gaussian function for modeling emission (log of fwhm). |
OpticalGaussian([name]) |
Gaussian function for modeling absorption (optical depth). |
Polynomial([name]) |
Polynomial model of order 5. |
Powerlaw([name]) |
Power-law model. |
Recombination([name]) |
Optically-thin recombination continuum model. |
SM([name]) |
Galactic extinction: the Savage & Mathis model. |
SMC([name]) |
SMC extinction: the Prevot et al. |
Seaton([name]) |
Galactic extinction: the Seaton model from Synphot. |
XGal([name]) |
Extragalactic extinction: Calzetti, Kinney and Storchi-Bergmann |
Class Inheritance Diagram¶
The sherpa.astro.xspec module¶
Support for XSPEC models.
Sherpa supports versions 12.10.1, 12.10.0, 12.9.1, and 12.9.0 of XSPEC [1], and can be built against the model library or the full application. There is no guarantee of support for older or newer versions of XSPEC.
To be able to use most routines from this module, the HEADAS environment variable must be set. The get_xsversion function can be used to return the XSPEC version - including patch level - the module is using:
>>> from sherpa.astro import xspec
>>> xspec.get_xsversion()
'12.10.1b'
Initializing XSPEC¶
The XSPEC model library is initalized so that the cosmology parameters are set to H_0=70, q_0=0.0, and lambda_0=0.73 (they can be changed with set_xscosmo).
The other settings - for example for the abundance and cross-section
tables - follow the standard rules for XSPEC. For XSPEC versions prior
to 12.10.1, this means that the abundance table uses the angr
setting and the cross sections the bcmc setting (see set_xsabund
and set_xsxsect for full details). As of XSPEC 12.10.1, the values
are now taken from the user’s XSPEC configuration file - either
~/.xspec/Xspec.init or $HEADAS/../spectral/manager/Xspec.init -
for these settings. The default value for the photo-ionization table
in this case is now vern rather than bcmc.
References
| [1] | https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/index.html |
Classes
XSagauss([name]) |
The XSPEC agauss model: gaussian line profile in wavelength space. |
XSapec([name]) |
The XSPEC apec model: APEC emission spectrum. |
XSbapec([name]) |
The XSPEC bapec model: velocity broadened APEC thermal plasma model. |
XSbbody([name]) |
The XSPEC bbody model: blackbody spectrum. |
XSbbodyrad([name]) |
The XSPEC bbodyrad model: blackbody spectrum, area normalized. |
XSbexrav([name]) |
The XSPEC bexrav model: reflected e-folded broken power law, neutral medium. |
XSbexriv([name]) |
The XSPEC bexriv model: reflected e-folded broken power law, ionized medium. |
XSbkn2pow([name]) |
The XSPEC bkn2pow model: broken power law with two breaks. |
XSbknpower([name]) |
The XSPEC bknpower model: broken power law. |
XSbmc([name]) |
The XSPEC bmc model: Comptonization by relativistic matter. |
XSbremss([name]) |
The XSPEC bremss model: thermal bremsstrahlung. |
XSbtapec([name]) |
The XSPEC btapec model: velocity broadened APEC emission spectrum with separate continuum and line temperatures. |
XSbvapec([name]) |
The XSPEC bvapec model: velocity broadened APEC thermal plasma model. |
XSbvtapec([name]) |
The XSPEC bvtapec model: velocity broadened APEC emission spectrum with separate continuum and line temperatures. |
XSbvvapec([name]) |
The XSPEC bvvapec model: velocity broadened APEC thermal plasma model. |
XSbvvtapec([name]) |
The XSPEC bvvtapec model: velocity broadened APEC emission spectrum with separate continuum and line temperatures. |
XSc6mekl([name]) |
The XSPEC c6mekl model: differential emission measure using Chebyshev representations with multi-temperature mekal. |
XSc6pmekl([name]) |
The XSPEC c6pmekl model: differential emission measure using Chebyshev representations with multi-temperature mekal. |
XSc6pvmkl([name]) |
The XSPEC c6pvmkl model: differential emission measure using Chebyshev representations with multi-temperature mekal. |
XSc6vmekl([name]) |
The XSPEC c6vmekl model: differential emission measure using Chebyshev representations with multi-temperature mekal. |
XScarbatm([name]) |
The XSPEC carbatm model: Nonmagnetic carbon atmosphere of a neutron star. |
XScemekl([name]) |
The XSPEC cemekl model: plasma emission, multi-temperature using mekal. |
XScevmkl([name]) |
The XSPEC cevmkl model: plasma emission, multi-temperature using mekal. |
XScflow([name]) |
The XSPEC cflow model: cooling flow. |
XScompLS([name]) |
The XSPEC compLS model: Comptonization, Lamb & Sanford. |
XScompPS([name]) |
The XSPEC compPS model: Comptonization, Poutanen & Svenson. |
XScompST([name]) |
The XSPEC compST model: Comptonization, Sunyaev & Titarchuk. |
XScompTT([name]) |
The XSPEC compTT model: Comptonization, Titarchuk. |
XScompbb([name]) |
The XSPEC compbb model: Comptonization, black body. |
XScompmag([name]) |
The XSPEC compmag model: Thermal and bulk Comptonization for cylindrical accretion onto the polar cap of a magnetized neutron star. |
XScomptb([name]) |
The XSPEC comptb model: Thermal and bulk Comptonization of a seed blackbody-like spectrum. |
XScompth([name]) |
The XSPEC compth model: Paolo Coppi’s hybrid (thermal/non-thermal) hot plasma emission models. |
XScplinear([name]) |
The XSPEC cplinear model: a non-physical piecewise-linear model for low count background spectra. |
XScutoffpl([name]) |
The XSPEC cutoffpl model: power law, high energy exponential cutoff. |
XSdisk([name]) |
The XSPEC disk model: accretion disk, black body. |
XSdiskbb([name]) |
The XSPEC diskbb model: accretion disk, multi-black body components. |
XSdiskir([name]) |
The XSPEC diskir model: Irradiated inner and outer disk. |
XSdiskline([name]) |
The XSPEC diskline model: accretion disk line emission, relativistic. |
XSdiskm([name]) |
The XSPEC diskm model: accretion disk with gas pressure viscosity. |
XSdisko([name]) |
The XSPEC disko model: accretion disk, inner, radiation pressure viscosity. |
XSdiskpbb([name]) |
The XSPEC diskpbb model: accretion disk, power-law dependence for T(r). |
XSdiskpn([name]) |
The XSPEC diskpn model: accretion disk, black hole, black body. |
XSeplogpar([name]) |
The XSPEC eplogpar model: log-parabolic blazar model with nu-Fnu normalization. |
XSeqpair([name]) |
The XSPEC eqpair model: Paolo Coppi’s hybrid (thermal/non-thermal) hot plasma emission models. |
XSeqtherm([name]) |
The XSPEC eqtherm model: Paolo Coppi’s hybrid (thermal/non-thermal) hot plasma emission models. |
XSequil([name]) |
The XSPEC equil model: collisional plasma, ionization equilibrium. |
XSexpdec([name]) |
The XSPEC expdec model: exponential decay. |
XSezdiskbb([name]) |
The XSPEC ezdiskbb model: multiple blackbody disk model with zero-torque inner boundary. |
XSgadem([name]) |
The XSPEC gadem model: plasma emission, multi-temperature with gaussian distribution of emission measure. |
XSgaussian([name]) |
The XSPEC gaussian model: gaussian line profile. |
XSgnei([name]) |
The XSPEC gnei model: collisional plasma, non-equilibrium, temperature evolution. |
XSgrad([name]) |
The XSPEC grad model: accretion disk, Schwarzschild black hole. |
XSgrbm([name]) |
The XSPEC grbm model: gamma-ray burst continuum. |
XShatm([name]) |
The XSPEC hatm model: Nonmagnetic hydrogen atmosphere of a neutron star. |
XSkerrbb([name]) |
The XSPEC kerrbb model: multi-temperature blackbody model for thin accretion disk around a Kerr black hole. |
XSkerrd([name]) |
The XSPEC kerrd model: optically thick accretion disk around a Kerr black hole. |
XSkerrdisk([name]) |
The XSPEC kerrdisk model: accretion disk line emission with BH spin as free parameter. |
XSlaor([name]) |
The XSPEC laor model: accretion disk, black hole emission line. |
XSlaor2([name]) |
The XSPEC laor2 model: accretion disk with broken-power law emissivity profile, black hole emission line. |
XSlogpar([name]) |
The XSPEC logpar model: log-parabolic blazar model. |
XSlorentz([name]) |
The XSPEC lorentz model: lorentz line profile. |
XSmeka([name]) |
The XSPEC meka model: emission, hot diffuse gas (Mewe-Gronenschild). |
XSmekal([name]) |
The XSPEC mekal model: emission, hot diffuse gas (Mewe-Kaastra-Liedahl). |
XSmkcflow([name]) |
The XSPEC mkcflow model: cooling flow, mekal. |
XSnei([name]) |
The XSPEC nei model: collisional plasma, non-equilibrium, constant temperature. |
XSnlapec([name]) |
The XSPEC nlapec model: continuum-only APEC emission spectrum. |
XSnpshock([name]) |
The XSPEC npshock model: shocked plasma, plane parallel, separate ion, electron temperatures. |
XSnsa([name]) |
The XSPEC nsa model: neutron star atmosphere. |
XSnsagrav([name]) |
The XSPEC nsagrav model: NS H atmosphere model for different g. |
XSnsatmos([name]) |
The XSPEC nsatmos model: NS Hydrogen Atmosphere model with electron conduction and self-irradiation. |
XSnsmax([name]) |
The XSPEC nsmax model: Neutron Star Magnetic Atmosphere. |
XSnsmaxg([name]) |
The XSPEC nsmaxg model: neutron star with a magnetic atmosphere. |
XSnsx([name]) |
The XSPEC nsx model: neutron star with a non-magnetic atmosphere. |
XSnteea([name]) |
The XSPEC nteea model: non-thermal pair plasma. |
XSnthComp([name]) |
The XSPEC nthComp model: Thermally comptonized continuum. |
XSoptxagn([name]) |
The XSPEC optxagn model: Colour temperature corrected disc and energetically coupled Comptonisation model for AGN. |
XSoptxagnf([name]) |
The XSPEC optxagnf model: Colour temperature corrected disc and energetically coupled Comptonisation model for AGN. |
XSpegpwrlw([name]) |
The XSPEC pegpwrlw model: power law, pegged normalization. |
XSpexmon([name]) |
The XSPEC pexmon model: neutral Compton reflection with self-consistent Fe and Ni lines. |
XSpexrav([name]) |
The XSPEC pexrav model: reflected powerlaw, neutral medium. |
XSpexriv([name]) |
The XSPEC pexriv model: reflected powerlaw, neutral medium. |
XSplcabs([name]) |
The XSPEC plcabs model: powerlaw observed through dense, cold matter. |
XSposm([name]) |
The XSPEC posm model: positronium continuum. |
XSpowerlaw([name]) |
The XSPEC powerlaw model: power law photon spectrum. |
XSpshock([name]) |
The XSPEC pshock model: plane-parallel shocked plasma, constant temperature. |
XSraymond([name]) |
The XSPEC raymond model: emission, hot diffuse gas, Raymond-Smith. |
XSredge([name]) |
The XSPEC redge model: emission, recombination edge. |
XSrefsch([name]) |
The XSPEC refsch model: reflected power law from ionized accretion disk. |
XSrnei([name]) |
The XSPEC rnei model: non-equilibrium recombining collisional plasma. |
XSsedov([name]) |
The XSPEC sedov model: sedov model, separate ion/electron temperature. |
XSsirf([name]) |
The XSPEC sirf model: self-irradiated funnel. |
XSslimbh([name]) |
The XSPEC slimbh model: Stationary slim accretion disk. |
XSsnapec([name]) |
The XSPEC snapec model: galaxy cluster spectrum using SN yields. |
XSsrcut([name]) |
The XSPEC srcut model: synchrotron spectrum, cutoff power law. |
XSsresc([name]) |
The XSPEC sresc model: synchrotron spectrum, cut off by particle escape. |
XSstep([name]) |
The XSPEC step model: step function convolved with gaussian. |
XStapec([name]) |
The XSPEC tapec model: APEC emission spectrum with separate continuum and line temperatures. |
XSvapec([name]) |
The XSPEC vapec model: APEC emission spectrum. |
XSvbremss([name]) |
The XSPEC vbremss model: thermal bremsstrahlung. |
XSvequil([name]) |
The XSPEC vequil model: collisional plasma, ionization equilibrium. |
XSvgadem([name]) |
The XSPEC vgadem model: plasma emission, multi-temperature with gaussian distribution of emission measure. |
XSvgnei([name]) |
The XSPEC vgnei model: collisional plasma, non-equilibrium, temperature evolution. |
XSvmcflow([name]) |
The XSPEC vmcflow model: cooling flow, mekal. |
XSvmeka([name]) |
The XSPEC vmeka model: emission, hot diffuse gas (Mewe-Gronenschild). |
XSvmekal([name]) |
The XSPEC vmekal model: emission, hot diffuse gas (Mewe-Kaastra-Liedahl). |
XSvnei([name]) |
The XSPEC vnei model: collisional plasma, non-equilibrium, constant temperature. |
XSvnpshock([name]) |
The XSPEC vnpshock model: shocked plasma, plane parallel, separate ion, electron temperatures. |
XSvoigt([name]) |
The XSPEC voigt model: Voigt line profile. |
XSvpshock([name]) |
The XSPEC vpshock model: plane-parallel shocked plasma, constant temperature. |
XSvraymond([name]) |
The XSPEC vraymond model: emission, hot diffuse gas, Raymond-Smith. |
XSvrnei([name]) |
The XSPEC vrnei model: non-equilibrium recombining collisional plasma. |
XSvsedov([name]) |
The XSPEC vsedov model: sedov model, separate ion/electron temperature. |
XSvtapec([name]) |
The XSPEC vtapec model: APEC emission spectrum with separate continuum and line temperatures. |
XSvvapec([name]) |
The XSPEC vvapec model: APEC emission spectrum. |
XSvvgnei([name]) |
The XSPEC vvgnei model: collisional plasma, non-equilibrium, temperature evolution. |
XSvvnei([name]) |
The XSPEC vvnei model: collisional plasma, non-equilibrium, constant temperature. |
XSvvnpshock([name]) |
The XSPEC vvnpshock model: shocked plasma, plane parallel, separate ion, electron temperatures. |
XSvvpshock([name]) |
The XSPEC vvpshock model: plane-parallel shocked plasma, constant temperature. |
XSvvrnei([name]) |
The XSPEC vvrnei model: non-equilibrium recombining collisional plasma. |
XSvvsedov([name]) |
The XSPEC vvsedov model: sedov model, separate ion/electron temperature. |
XSvvtapec([name]) |
The XSPEC vvtapec model: APEC emission spectrum with separate continuum and line temperatures. |
XSzagauss([name]) |
The XSPEC zagauss model: gaussian line profile in wavelength space. |
XSzbbody([name]) |
The XSPEC zbbody model: blackbody spectrum. |
XSzbremss([name]) |
The XSPEC zbremss model: thermal bremsstrahlung. |
XSzgauss([name]) |
The XSPEC zgauss model: gaussian line profile. |
XSzpowerlw([name]) |
The XSPEC zpowerlw model: redshifted power law photon spectrum. |
XSSSS_ice([name]) |
The XSPEC sss_ice model: Einstein SSS ice absorption. |
XSTBabs([name]) |
The XSPEC TBabs model: ISM grain absorption. |
XSTBfeo([name]) |
The XSPEC TBfeo model: ISM grain absorption. |
XSTBgas([name]) |
The XSPEC TBgas model: ISM grain absorption. |
XSTBgrain([name]) |
The XSPEC TBgrain model: ISM grain absorption. |
XSTBpcf([name]) |
The XSPEC TBpcf model: ISM grain absorption. |
XSTBrel([name]) |
The XSPEC TBrel model: ISM grain absorption. |
XSTBvarabs([name]) |
The XSPEC TBvarabs model: ISM grain absorption. |
XSabsori([name]) |
The XSPEC absori model: ionized absorber. |
XSacisabs([name]) |
The XSPEC acisabs model: Chandra ACIS q.e. |
XScabs([name]) |
The XSPEC cabs model: Optically-thin Compton scattering. |
XSconstant([name]) |
The XSPEC constant model: energy-independent factor. |
XScyclabs([name]) |
The XSPEC cyclabs model: absorption line, cyclotron. |
XSdust([name]) |
The XSPEC dust model: dust scattering. |
XSedge([name]) |
The XSPEC edge model: absorption edge. |
XSexpabs([name]) |
The XSPEC expabs model: exponential roll-off at low E. |
XSexpfac([name]) |
The XSPEC expfac model: exponential modification. |
XSgabs([name]) |
The XSPEC gabs model: gaussian absorption line. |
XSheilin([name]) |
The XSPEC heilin model: Voigt absorption profiles for He I series. |
XShighecut([name]) |
The XSPEC highecut model: high-energy cutoff. |
XShrefl([name]) |
The XSPEC hrefl model: reflection model. |
XSismabs([name]) |
The XSPEC ismabs model: A high resolution ISM absorption model with variable columns for individual ions. |
XSlyman([name]) |
The XSPEC lyman model: Voigt absorption profiles for H I or He II Lyman series. |
XSnotch([name]) |
The XSPEC notch model: absorption line, notch. |
XSpcfabs([name]) |
The XSPEC pcfabs model: partial covering fraction absorption. |
XSphabs([name]) |
The XSPEC phabs model: photoelectric absorption. |
XSplabs([name]) |
The XSPEC plabs model: power law absorption. |
XSpwab([name]) |
The XSPEC pwab model: power-law distribution of neutral absorbers. |
XSredden([name]) |
The XSPEC redden model: interstellar extinction. |
XSsmedge([name]) |
The XSPEC smedge model: smeared edge. |
XSspexpcut([name]) |
The XSPEC spexpcut model: super-exponential cutoff absorption. |
XSspline([name]) |
The XSPEC spline model: spline modification. |
XSswind1([name]) |
The XSPEC swind1 model: absorption by partially ionized material with large velocity shear. |
XSuvred([name]) |
The XSPEC uvred model: interstellar extinction, Seaton Law. |
XSvarabs([name]) |
The XSPEC varabs model: photoelectric absorption. |
XSvphabs([name]) |
The XSPEC vphabs model: photoelectric absorption. |
XSwabs([name]) |
The XSPEC wabs model: photoelectric absorption, Wisconsin cross-sections. |
XSwndabs([name]) |
The XSPEC wndabs model: photo-electric absorption, warm absorber. |
XSxion([name]) |
The XSPEC xion model: reflected spectrum of photo-ionized accretion disk/ring. |
XSxscat([name]) |
The XSPEC xscat model: dust scattering. |
XSzTBabs([name]) |
The XSPEC zTBabs model: ISM grain absorption. |
XSzbabs([name]) |
The XSPEC zbabs model: EUV ISM attenuation. |
XSzdust([name]) |
The XSPEC zdust model: extinction by dust grains. |
XSzedge([name]) |
The XSPEC zedge model: absorption edge. |
XSzhighect([name]) |
The XSPEC zhighect model: high-energy cutoff. |
XSzigm([name]) |
The XSPEC zigm model: UV/Optical attenuation by the intergalactic medium. |
XSzpcfabs([name]) |
The XSPEC zpcfabs model: partial covering fraction absorption. |
XSzphabs([name]) |
The XSPEC zphabs model: photoelectric absorption. |
XSzredden([name]) |
The XSPEC zredden model: redshifted version of redden. |
XSzsmdust([name]) |
The XSPEC zsmdust model: extinction by dust grains in starburst galaxies. |
XSzvarabs([name]) |
The XSPEC zvarabs model: photoelectric absorption. |
XSzvfeabs([name]) |
The XSPEC zvfeabs model: photoelectric absorption with free Fe edge energy. |
XSzvphabs([name]) |
The XSPEC zvphabs model: photoelectric absorption. |
XSzwabs([name]) |
The XSPEC zwabs model: photoelectric absorption, Wisconsin cross-sections. |
XSzwndabs([name]) |
The XSPEC zwndabs model: photo-electric absorption, warm absorber. |
XSzxipcf([name]) |
The XSPEC zxipcf model: partial covering absorption by partially ionized material. |
Class Inheritance Diagram¶
Reference/API¶
This section describes the classes used to create models and the Available Models section above contains the classes that implement various models.
The sherpa.models.model module¶
Classes
Model(name[, pars]) |
The base class for Sherpa models. |
ArithmeticConstantModel(val[, name]) |
|
ArithmeticFunctionModel(func) |
|
ArithmeticModel(name[, pars]) |
|
CompositeModel(name, parts) |
|
BinaryOpModel(lhs, rhs, op, opstr) |
|
FilterModel(model, filter) |
|
MultigridSumModel(models) |
|
NestedModel(outer, inner, *otherargs, …) |
|
RegriddableModel1D(name[, pars]) |
|
RegriddableModel2D(name[, pars]) |
|
SimulFitModel(name, parts) |
Store multiple models. |
UnaryOpModel(arg, op, opstr) |
Class Inheritance Diagram¶
The sherpa.models.parameter module¶
Support for model parameter values.
Classes
Parameter(modelname, name, val[, min, max, …]) |
Represent a model parameter. |
CompositeParameter(name, parts) |
|
BinaryOpParameter(lhs, rhs, op, opstr) |
|
ConstantParameter(value) |
|
UnaryOpParameter(arg, op, opstr) |
Class Inheritance Diagram¶
The sherpa.models.regrid module¶
Classes
Axis(lo, hi) |
Class for representing N-D axes objects, for both “integrated” and “non-integrated” datasets |
EvaluationSpace1D([x, xhi]) |
Class for 1D Evaluation Spaces. |
EvaluationSpace2D([x, y, xhi, yhi]) |
Class for 2D Evaluation Spaces. |
ModelDomainRegridder1D([evaluation_space, name]) |
Allow 1D models to be evaluated on a different grid. |
ModelDomainRegridder2D([evaluation_space, name]) |
Allow 2D models to be evaluated on a different grid. |
Functions
rebin_2d(y, from_space, to_space) |
|
rebin_int(array, scale_x, scale_y) |
Rebin array by an integer scale on both x and y |
rebin_no_int(array[, dimensions, scale]) |
Rebin the array, conserving flux. |
Class Inheritance Diagram¶
The sherpa.instrument module¶
Classes
ConvolutionKernel(kernel[, name]) |
|
ConvolutionModel(lhs, rhs, psf) |
|
Kernel(dshape, kshape[, norm, frozen, …]) |
Base class for convolution kernels |
PSFKernel(dshape, kshape[, is_model, norm, …]) |
class for PSF convolution kernels |
PSFModel([name, kernel]) |
|
RadialProfileKernel(dshape, kshape[, …]) |
class for 1D radial profile PSF convolution kernels |
Class Inheritance Diagram¶
The sherpa.models.template module¶
Classes
TemplateModel([name, pars, parvals, templates]) |
|
InterpolatingTemplateModel(name, template_model) |
|
KNNInterpolator(name, template_model[, k, order]) |
|
Template(*args, **kwargs) |
Class Inheritance Diagram¶
The sherpa.astro.instrument module¶
Models of common Astronomical models, particularly in X-rays.
The models in this module include support for instrument models that describe how X-ray photons are converted to measurable properties, such as Pulse-Height Amplitudes (PHA) or Pulse-Invariant channels. These ‘responses’ are assumed to follow OGIP standards, such as [1].
References
| [1] | OGIP Calibration Memo CAL/GEN/92-002, “The Calibration Requirements for Spectral Analysis (Definition of RMF and ARF file formats)”, Ian M. George1, Keith A. Arnaud, Bill Pence, Laddawan Ruamsuwan and Michael F. Corcoran, https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html |
Classes
ARF1D(arf[, pha, rmf]) |
|
PileupResponse1D(pha, pileup_model) |
|
RMF1D(rmf[, pha, arf]) |
|
Response1D(pha) |
|
MultipleResponse1D(pha) |
|
PSFModel([name, kernel]) |
|
ARFModel(arf, model) |
Base class for expressing ARF convolution in model expressions. |
ARFModelNoPHA(arf, model) |
ARF convolution model without associated PHA data set. |
ARFModelPHA(arf, pha, model) |
ARF convolution model with associated PHA data set. |
RMFModel(rmf, model) |
Base class for expressing RMF convolution in model expressions. |
RMFModelNoPHA(rmf, model) |
RMF convolution model without an associated PHA data set. |
RMFModelPHA(rmf, pha, model) |
RMF convolution model with associated PHA data set. |
RSPModel(arf, rmf, model) |
Base class for expressing RMF + ARF convolution in model expressions |
RSPModelNoPHA(arf, rmf, model) |
RMF + ARF convolution model without associated PHA data set. |
RSPModelPHA(arf, rmf, pha, model) |
RMF + ARF convolution model with associated PHA. |
MultiResponseSumModel(source, pha) |
|
PileupRMFModel(rmf, model[, pha]) |
Functions
create_arf(elo, ehi[, specresp, exposure, …]) |
Create an ARF. |
create_delta_rmf(rmflo, rmfhi[, offset, …]) |
Create an ideal RMF. |
create_non_delta_rmf(rmflo, rmfhi, fname[, …]) |
Create a RMF using a matrix from a file. |
Class Inheritance Diagram¶
The sherpa.astro.xspec module¶
Support for XSPEC models.
Sherpa supports versions 12.10.1, 12.10.0, 12.9.1, and 12.9.0 of XSPEC [1], and can be built against the model library or the full application. There is no guarantee of support for older or newer versions of XSPEC.
To be able to use most routines from this module, the HEADAS environment variable must be set. The get_xsversion function can be used to return the XSPEC version - including patch level - the module is using:
>>> from sherpa.astro import xspec
>>> xspec.get_xsversion()
'12.10.1b'
Initializing XSPEC¶
The XSPEC model library is initalized so that the cosmology parameters are set to H_0=70, q_0=0.0, and lambda_0=0.73 (they can be changed with set_xscosmo).
The other settings - for example for the abundance and cross-section
tables - follow the standard rules for XSPEC. For XSPEC versions prior
to 12.10.1, this means that the abundance table uses the angr
setting and the cross sections the bcmc setting (see set_xsabund
and set_xsxsect for full details). As of XSPEC 12.10.1, the values
are now taken from the user’s XSPEC configuration file - either
~/.xspec/Xspec.init or $HEADAS/../spectral/manager/Xspec.init -
for these settings. The default value for the photo-ionization table
in this case is now vern rather than bcmc.
References
| [1] | https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/index.html |
This document describes the base classes for XSPEC models, and the utility routines - such as querying and retrieving the abundance table information. The models provided by XSPEC are described in The sherpa.astro.xspec module.
Classes
XSModel(name[, pars])The base class for XSPEC models. XSAdditiveModel(name[, pars])The base class for XSPEC additive models. XSMultiplicativeModel(name[, pars])The base class for XSPEC multiplicative models. XSTableModel(filename[, name, parnames, …])Interface to XSPEC table models. Functions
get_xsabund([element])Return the X-Spec abundance setting or elemental abundance. get_xschatter()Return the chatter level used by X-Spec. get_xscosmo()Return the X-Spec cosmology settings. get_xspath_manager()Return the path to the files describing the XSPEC models. get_xspath_model()Return the path to the model data files. get_xsstate()Return the state of the XSPEC module. get_xsversion()Return the version of the X-Spec model library in use. get_xsxsect()Return the cross sections used by X-Spec models. get_xsxset(name)Return the X-Spec model setting. read_xstable_model(modelname, filename)Create a XSPEC table model. set_xsabund(abundance)Set the elemental abundances used by X-Spec models. set_xschatter(level)Set the chatter level used by X-Spec. set_xscosmo(h0, q0, l0)Set the cosmological parameters used by X-Spec models. set_xspath_manager(path)Set the path to the files describing the XSPEC models. set_xsstate(state)Restore the state of the XSPEC module. set_xsxsect(name)Set the cross sections used by X-Spec models. set_xsxset(name, value)Set a X-Spec model setting.
Class Inheritance Diagram¶
What statistic is to be used?¶
The statistic object defines the numerical quantity which describes the “quality” of the fit, where the definition is that as the statistic value decreases, the model is getting to better describe the data. It is the statistic value that the optimiser uses to determine the “best-fit” parameter settings.
A simple example is the least-squares statistic which, if the data
and model points are \(d_i\) and \(m_i\) respectively,
where the suffix indicates the bin number, then the overall
statistic value is \(s = \sum_i (d_i - m_i)^2\). This is
provided by the LeastSq class, and
an example of its use is shown below. First we import the
classes we need:
>>> import numpy as np
>>> from sherpa.data import Data1D
>>> from sherpa.models.basic import Gauss1D
>>> from sherpa.stats import LeastSq
As the data I use a one-dimensional gaussian with normally-distributed noise:
>>> np.random.seed(0)
>>> x = np.linspace(-5., 5., 200)
>>> gmdl = Gauss1D()
>>> gmdl.fwhm = 1.9
>>> gmdl.pos = 1.3
>>> gmdl.ampl = 3
>>> y = gmdl(x) + np.random.normal(0., 0.2, x.shape)
>>> d = Data1D('stat-example', x, y)
The statistic value, along with per-bin values, is returned by
the calc_stat() method, so we
can find the statistic value of this parameter set with:
>>> stat = LeastSq()
>>> s = stat.calc_stat(d, mdl)
>>> print("Statistic value = {}".format(s[0]))
Statistic value = 8.38666216358492
As the FWHM is varied about the true value we can see that the statistic increases:
>>> fwhm = np.arange(0.5, 3.0, 0.2)
>>> sval = []
>>> for f in fwhm:
... gmdl.fwhm = f
... sval.append(stat.calc_stat(d, gmdl)[0])
...
>>> plt.plot(fwhm, sval)
>>> plt.xlabel('FWHM')
>>> plt.ylabel('Statistic')
The statistic classes provided in Sherpa are given below, and cover a range of possibilities, such as: least-square for when there is no knowledge about the errors on each point, a variety of chi-square statistics for when the errors are assumed to be gaussian, and a variety of maximum-likelihood estimators for Poisson-distributed data. It is also possible to add your own statistic class.
Reference/API¶
The sherpa.stats module¶
Classes
Stat(name) |
The base class for calculating a statistic given data and model. |
Chi2([name]) |
Chi Squared statistic. |
LeastSq([name]) |
Least Squared Statistic. |
Chi2ConstVar([name]) |
Chi Squared with constant variance. |
Chi2DataVar([name]) |
Chi Squared with data variance. |
Chi2Gehrels([name]) |
Chi Squared with Gehrels variance. |
Chi2ModVar([name]) |
Chi Squared with model amplitude variance. |
Chi2XspecVar([name]) |
Chi Squared with data variance (XSPEC style). |
Likelihood([name]) |
Maximum likelihood function |
Cash([name]) |
Maximum likelihood function. |
CStat([name]) |
Maximum likelihood function (XSPEC style). |
WStat([name]) |
Maximum likelihood function including background (XSPEC style). |
UserStat([statfunc, errfunc, name]) |
Support simple user-supplied statistic calculations. |
Class Inheritance Diagram¶
Optimisers: How to improve the current parameter values¶
The optimiser varies the model parameters in an attempt to find the solution which minimises the chosen statistic.
In general it is expected that the optimiser will be used by
a Fit object to
perform the fit, but
it can be used directly using the
fit() method. The optimiser
object allows configuration values to be changed which can
tweak the behavior; for instance the tolerance to determine whether
the fit has converged, the maximum number of iterations to use,
or how much information to display whilst optimising a model.
As an example, the default parameter values for the
Levenberg-Marquardt
optimiser are:
>>> from sherpa.optmethods.LevMar
>>> lm = LevMar()
>>> print(lm)
name = levmar
ftol = 1.19209289551e-07
xtol = 1.19209289551e-07
gtol = 1.19209289551e-07
maxfev = None
epsfcn = 1.19209289551e-07
factor = 100.0
verbose = 0
These settings are available both as fields of the object and via
the config dictionary
field.
Additional optimisers can be built by extending from the
sherpa.optmethods.OptMethod class. This can be used
to provide access to external packages such as
CERN’s MINUIT optimisation library.
Choosing an optimiser¶
Warning
The following may not correctly represent Sherpa’s current capabilities, so please take care when interpreting this section.
The following information is adapted from a memo written by Mark Birkinshaw (1998).
The minimization of mathematical functions is a difficult operation. A general function \(f({\bf x})\) of the vector argument \(\bf x\) may have many isolated local minima, non-isolated minimum hypersurfaces, or even more complicated topologies. No finite minimization routine can guarantee to locate the unique, global, minimum of \(f({\bf x})\) without being fed intimate knowledge about the function by the user.
This does not mean that minimization is a hopeless task. For many problems there are techniques which will locate a local minimum which may be “close enough” to the global minimum, and there are techniques which will find the global minimum a large fraction of the time (in a probabilistic sense). However, the reader should be aware of my philosophy is that there is no “best” algorithm for finding the minimum of a general function. Instead, Sherpa provides tools which will allow the user to look at the overall behavior of the function and find plausible local minima, which will often contain the physically-meaningful minimum in the types of problem with which Sherpa deals.
In general, the best assurance that the correct minimum has been found in a particular calculation is careful examination of the nature of the solution (e.g., by plotting a fitted function over data), and some confidence that the full region that the minimum may lie in has been well searched by the algorithm used. This document seeks to give the reader some information about what the different choices of algorithm will mean in terms of run-time and confidence of locating a good minimum.
Some points to take away from the discussions in the rest of this document.
- Never accept the result of a minimization using a single optimization run; always test the minimum using a different method.
- Check that the result of the minimization does not have parameter values at the edges of the parameter space. If this happens, then the fit must be disregarded since the minimum lies outside the space that has been searched, or the minimization missed the minimum.
- Get a feel for the range of values of the target function (in Sherpa this is the fit statistic), and the stability of the solution, by starting the minimization from several different parameter values.
- Always check that the minimum “looks right” by visualizing the model and the data.
Sherpa contains two types of routine for minimizing a fit statistic. I will call them the “single-shot” routines, which start from a guessed set of parameters, and then try to improve the parameters in a continuous fashion, and the “scatter-shot” routines, which try to look at parameters over the entire permitted hypervolume to see if there are better minima than near the starting guessed set of parameters.
Single-shot techniques¶
As the reader might expect, the single-shot routines are relatively quick, but depend critically on the guessed initial parameter values \({\bf x}_0\) being near (in some sense) to the minimum \({\bf x}_{\rm min}\). All the single-shot routines investigate the local behaviour of the function near \({\bf x}_0\), and then make a guess at the best direction and distance to move to find a better minimum. After testing at the new point, they accept that point as the next guess, \({\bf x}_1\), if the fit statistic is smaller than at the first point, and modify the search procedure if it isn’t smaller. The routines continue to run until one of the following occurs:
- all search directions result in an increased value of the fit statistic;
- an excessive number of steps have been taken; or
- something strange happens to the fit statistic (e.g., it turns out to be discontinuous in some horrible way).
This description indicates that for the single-shot routines, there is a considerable emphasis on the initial search position, \({\bf x}_0\), being reasonable. It may also be apparent that the values of these parameters should be moderate; neither too small (\(10^{-12}\), say), nor too large (\(10^{12}\), say). This is because the initial choice of step size in moving from \({\bf x}_0\) towards the next improved set of parameters, \({\bf x}_1\), is based on the change in the fit statistic, \(f({\bf x})\) as components of \({\bf x}\) are varied by amounts \({\cal O}(1)\). If \(f\) varies little as \({\bf x}\) is varied by this amount, then the calculation of the distance to move to reach the next root may be inaccurate. On the other hand, if \(f\) has a lot of structure (several maxima and minima) as \({\bf x}\) is varied by the initial step size, then these single-shot minimizers may mistakenly jump entirely over the “interesting” region of parameter space.
These considerations suggest that the user should arrange that the search vector is scaled so that the range of parameter space to be searched is neither too large nor too small. To take a concrete example, it would not be a good idea to have \(x_7\) parameterize the Hydrogen column density (\(N_{\rm H}\)) in a spectral fit, with an initial guess of \(10^{20}\ {\rm cm}^{-2}\), and a search range (in units of \({\rm cm}^{-2}\)) of \(10^{16}\) to \(10^{24}\). The minimizers will look for variations in the fit statistic as \(N_{\rm H}\) is varied by \(1\ {\rm cm}^{-2}\), and finding none (to the rounding accuracy likely for the code), will conclude that \(x_7\) is close to being a null parameter and can be ignored in the fitting. It would be much better to have \(x_7 = \log_{10}(N_{\rm H})\), with a search range of 16 to 24. Significant variations in the fit statistic will occur as \(x_7\) is varied by \(\pm 1\), and the code has a reasonable chance of finding a useful solution.
Bearing this in mind, the single-shot minimizers in Sherpa are listed below:
NelderMead- This technique - also known as Simplex - creates a polyhedral search element around the initial position, \({\bf x}_0\), and then grows or shrinks in particular directions while crawling around parameter space, to try to place a minimum within the final search polyhedron. This technique has some hilarious ways of getting stuck in high-dimension parameter spaces (where the polyhedron can become a strange shape), but is very good at finding minima in regions where the fit statistic has a moderately well-defined topology. Since it works in a different way than Levenberg-Marquardt minimization, a good strategy is to combine both minimization to test whether an apparent minimum found by one technique is stable when searched by the other. I regard NelderMead searching as good in smooth and simple parameter spaces, particularly when looking at regions where the fit statistic depends on a parameter in a linear or parabolic fashion, and bad where surfaces of equal value of the fit statistic are complicated. In either case, it is essential that the initial size of the polyhedron (with sides of length 1 unit) is a smallish fraction of the search space.
Levenberg-Marquardt- This can be considered to be a censored maximum-gradients technique which, starting from a first guess, moves towards a minimum by finding a good direction in which to move, and calculating a sensible distance to go. Its principal drawback is that to calculate the distance to move it has to make some assumptions about how large a step size to take, and hence there is an implicit assumption that the search space is reasonably well scaled (to \(\pm 10\) units in each of the search directions, say). It is also important that in finding these gradients, the steps do not miss a lot of important structure; i.e. there should not be too many subsidiary minima. The search directions and distances to move are based on the shape of the target function near the initial guessed minimum, \({\bf x}_0\), with progressive movement towards the dominant local minimum. Since this technique uses information about the local curvature of the fit statistic as well as its local gradients, the approach tends to stabilize the result in somce cases. I regard the techniques implemented in Sherpa as being good minimum-refiners for simple local topologies, since more assumptions about topology are made than in the NelderMead approach, but bad at finding global minima for target functions with complicated topologies.
Scatter-shot techniques¶
Although a bit ad hoc, these techniques attempt to locate a decent minimum over the entire range of the search parameter space. Because they involve searching a lot of the parameter space, they involve many function evaluations, and are somewhere between quite slow and incredibly-tediously slow.
The routines are listed below:
GridSearch- This routine simply searches a grid in each of the search parameters,
where the spacing is uniform between the minimum and maximum
value of each parameter. There is an option to refine the fit
at each point, by setting the
methodattribute to one of the single-shot optimisers, but this is not set by default, as it can significantly increase the time required to fit the data. The coarseness of the grid sets how precise a root will be found, and if the fit statistic has significant structure on a smaller scale, then the grid-searcher will miss it completely. This is a good technique for finding an approximation to the minimum for a slowly-varying function. It is a bad technique for getting accurate estimates of the location of a minimum, or for examining a fit statistic with lots of subsidiary maxima and minima within the search space. It is intended for use withtemplate models. Monte Carlo- This is a simple population based, stochastic function minimizer. At each iteration it combines population vectors - each containing a set of parameter values - using a weighted difference. This optimiser can be used to find solutions to complex search spaces but is not guaranteed to find a global minimum. It is over-kill for relatively simple problems.
Summary and best-buy strategies¶
Overall, the single-shot methods are best regarded as ways of refining minima located in other ways: from good starting guesses, or from the scatter-shot techniques. Using intelligence to come up with a good first-guess solution is the best approach, when the single-shot refiners can be used to get accurate values for the parameters at the minimum. However, I would certainly recommend running at least a second single-shot minimizer after the first, to get some indication that one set of assumptions about the shape of the minimum is not compromising the solution. It is probably best if the code rescales the parameter range between minimizations, so that a completely different sampling of the function near the trial minimum is being made.
| Optimiser | Type | Speed | Commentary |
|---|---|---|---|
| NelderMead | single-shot | fast | OK for refining minima |
| Levenberg-Marquardt | single-shot | fast | OK for refining minima |
| GridSearch | scatter-shot | slow | OK for smooth functions |
| Monte Carlo | scatter-shot | very slow | Good in many cases |
Reference/API¶
The sherpa.optmethods module¶
Classes
OptMethod(name, optfunc) |
Base class for the optimisers. |
LevMar([name]) |
Levenberg-Marquardt optimization method. |
NelderMead([name]) |
Nelder-Mead Simplex optimization method. |
MonCar([name]) |
Monte Carlo optimization method. |
GridSearch([name]) |
Grid Search optimization method. |
Class Inheritance Diagram¶
Fitting the data¶
The Fit class takes the
data and model expression
to be fit, and uses the
optimiser to minimise the
chosen statistic. The basic
approach is to:
- create a
Fitobject;- call its
fit()method one or more times, potentially varying themethodattribute to change the optimiser;- inspect the
FitResultsobject returned byfit()to extract information about the fit;- review the fit quality, perhaps re-fitting with a different set of parameters or using a different optimiser;
- once the “best-fit” solution has been found, calculate error estimates by calling the
est_errors()method, which returns aErrorEstResultsobject detailing the results;- visualize the parameter errors with classes such as
IntervalProjectionandRegionProjection.
The following discussion uses a one-dimensional data set with gaussian errors (it was simulated with gaussian noise with \(\sigma = 5\)):
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from sherpa.data import Data1D
>>> d = Data1D('fit example', [-13, -5, -3, 2, 7, 12],
... [102.3, 16.7, -0.6, -6.7, -9.9, 33.2],
... np.ones(6) * 5)
It is going to be fit with the expression:
which is represented by the Polynom1D
model:
>>> from sherpa.models.basic import Polynom1D
>>> mdl = Polynom1D()
To start with, just the \(c_0\) and \(c_2\) terms are used in the fit:
>>> mdl.c2.thaw()
>>> print(mdl)
polynom1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
polynom1d.c0 thawed 1 -3.40282e+38 3.40282e+38
polynom1d.c1 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c2 thawed 0 -3.40282e+38 3.40282e+38
polynom1d.c3 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c4 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c5 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c6 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c7 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c8 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.offset frozen 0 -3.40282e+38 3.40282e+38
Creating a fit object¶
A Fit object requires both a
data set and a
model object, and will use the
Chi2Gehrels
statistic with the
LevMar
optimiser
unless explicitly over-riden with the stat and
method parameters (when creating the object) or the
stat and
method attributes
(once the object has been created).
>>> from sherpa.fit import Fit
>>> f = Fit(d, mdl)
>>> print(f)
data = fit example
model = polynom1d
stat = Chi2Gehrels
method = LevMar
estmethod = Covariance
>>> print(f.data)
name = fit example
x = [-13, -5, -3, 2, 7, 12]
y = [102.3, 16.7, -0.6, -6.7, -9.9, 33.2]
staterror = Float64[6]
syserror = None
>>> print(f.model)
polynom1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
polynom1d.c0 thawed 1 -3.40282e+38 3.40282e+38
polynom1d.c1 thawed 0 -3.40282e+38 3.40282e+38
polynom1d.c2 thawed 0 -3.40282e+38 3.40282e+38
polynom1d.c3 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c4 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c5 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c6 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c7 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c8 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.offset frozen 0 -3.40282e+38 3.40282e+38
The fit object stores references to objects, such as the model, which means that the fit object reflects the current state, and not just the values when it was created or used. For example, in the following the model is changed directly, and the value stored in the fit object is also changed:
>>> f.model.c2.val
0.0
>>> mdl.c2 = 1
>>> f.model.c2.val
1.0
Using the optimiser and statistic¶
With a Fit object can calculate the statistic value for
the current data and model
(calc_stat()),
summarise how well the current model represents the
data (calc_stat_info()),
calculate the per-bin chi-squared value (for chi-square
statistics; calc_stat_chisqr()),
fit the model to the data
(fit() and
simulfit()), and
calculate the parameter errors
(est_errors()).
>>> print("Starting statistic: {:.3f}".format(f.calc_stat()))
Starting statistic: 840.251
>>> sinfo1 = f.calc_stat_info()
>>> print(sinfo1)
name =
ids = None
bkg_ids = None
statname = chi2
statval = 840.2511999999999
numpoints = 6
dof = 4
qval = 1.46616165292e-180
rstat = 210.06279999999998
The fields in the StatInfoResults depend on
the choice of statistic; in particular,
rstat and
qval are set to None
if the statistic is not based on chi-square. The current data set
has a reduced statistic of \(\sim 200\) which indicates
that the fit is not very good (if the error bars in the
dataset are correct then a good fit is indicated by
a reduced statistic \(\simeq 1\) for
chi-square based statistics).
As with the model and statistic, if the data object is changed then
the results of calculations made on the fit object can be changed.
As shown below, when one data point is
removed, the calculated statistics
are changed (such as the
numpoints value).
>>> d.ignore(0, 5)
>>> sinfo2 = f.calc_stat_info()
>>> d.notice()
>>> sinfo1.numpoints
6
>>> sinfo2.numpoints
5
Note
The objects returned by the fit methods, such as
StatInfoResults and
FitResults, do not contain references
to any of the input objects. This means that the values stored
in these objects are not changed when the input values change.
The fit() method uses the optimiser to
vary the
thawed parameter values
in the model in an attempt to minimize the statistic value.
The method returns a
FitResults object which contains
information on the fit, such as whether it succeeded (found a minimum,
succeeded),
the start and end statistic value
(istatval and
statval),
and the best-fit parameter values
(parnames and
parvals).
>>> res = f.fit()
>>> if res.succeeded: print("Fit succeeded")
Fit succeeded
The return value has a format() method which
provides a summary of the fit:
>>> print(res.format())
Method = levmar
Statistic = chi2
Initial fit statistic = 840.251
Final fit statistic = 96.8191 at function evaluation 6
Data points = 6
Degrees of freedom = 4
Probability [Q-value] = 4.67533e-20
Reduced statistic = 24.2048
Change in statistic = 743.432
polynom1d.c0 -11.0742
polynom1d.c2 0.503612
The information is also available directly as fields of the
FitResults object:
>>> print(res)
datasets = None
itermethodname = none
methodname = levmar
statname = chi2
succeeded = True
parnames = ('polynom1d.c0', 'polynom1d.c2')
parvals = (-11.074165156709268, 0.50361247735062253)
statval = 96.8190901009578
istatval = 840.2511999999999
dstatval = 743.432109899
numpoints = 6
dof = 4
qval = 4.67533320771e-20
rstat = 24.20477252523945
message = successful termination
nfev = 6
The reduced chi square fit is now lower, \(\sim 25\), but still not close to 1.
Visualizing the fit¶
The DataPlot, ModelPlot,
and FitPlot classes can be used to
see how well the model represents the data.
>>> from sherpa.plot import DataPlot, ModelPlot
>>> dplot = DataPlot()
>>> dplot.prepare(f.data)
>>> mplot = ModelPlot()
>>> mplot.prepare(f.data, f.model)
>>> dplot.plot()
>>> mplot.overplot()
As can be seen, the model isn’t able to fully describe the data. One
check to make here is to change the optimiser in case the fit is stuck
in a local minimum. The default optimiser is
LevMar:
>>> f.method.name
'levmar'
>>> original_method = f.method
This can be changed to NelderMead
and the data re-fit.
>>> from sherpa.optmethods import NelderMead
>>> f.method = NelderMead()
>>> resn = f.fit()
In this case the statistic value has not changed, as shown by
dstatval being zero:
>>> print("Change in statistic: {}".format(resn.dstatval))
Change in statistic: 0.0
Note
An alternative approach is to create a new Fit object, with the
new method, and use that instead of changing the
method attribute. For instance:
>>> fit2 = Fit(d, mdl, method=NelderMead())
>>> fit2.fit()
Adjusting the model¶
This suggests that the problem is not that a local minimum has been found, but that the model is not expressive enough to represent the data. One possible approach is to change the set of points used for the comparison, either by removing data points - perhaps because they are not well calibrated or there are known to be issues - or adding extra ones (by removing a previously-applied filter). The approach taken here is to change the model being fit; in this case by allowing the linear term (\(c_1\)) of the polynomial to be fit, but a new model could have been tried.
>>> mdl.c1.thaw()
For the remainder of the analysis the original (Levenberg-Marquardt) optimiser will be used:
>>> f.method = original_method
With \(c_1\) allowed to vary, the fit finds a much better solution, with a reduced chi square value of \(\simeq 1.3\):
>>> res2 = f.fit()
>>> print(res2.format())
Method = levmar
Statistic = chi2
Initial fit statistic = 96.8191
Final fit statistic = 4.01682 at function evaluation 8
Data points = 6
Degrees of freedom = 3
Probability [Q-value] = 0.259653
Reduced statistic = 1.33894
Change in statistic = 92.8023
polynom1d.c0 -9.38489
polynom1d.c1 -2.41692
polynom1d.c2 0.478273
The previous plot objects can be used, but the model plot has to be
updated to reflect the new model values. Three new plot styles are used:
FitPlot shows both the data and model values,
DelchiPlot to show the residuals, and
SplotPlot to control the layout of the plots:
>>> from sherpa.plot import DelchiPlot, FitPlot, SplitPlot
>>> fplot = FitPlot()
>>> rplot = DelchiPlot()
>>> splot = SplitPlot()
>>> mplot.prepare(f.data, f.model)
>>> fplot.prepare(dplot, mplot)
>>> splot.addplot(fplot)
>>> rplot.prepare(f.data, f.model, f.stat)
>>> splot.addplot(rplot)
The residuals plot shows, on the ordinate, \(\sigma = (d - m) / e\) where \(d\), \(m\), and \(e\) are the data, model, and error value for each bin respectively.
Refit from a different location¶
It may be useful to repeat the fit starting the model off at
different parameter locations to check that the fit solution is
robust. This can be done manually, or the
reset() method used to change
the parameters back to
the last values they were explicitly set to:
>>> mdl.reset()
Rather than print out all the components, most of which are fixed at
0, the first three can be looped over using the model’s
pars attribute:
>>> [(p.name, p.val, p.frozen) for p in mdl.pars[:3]]
[('c0', 1.0, False), ('c1', 0.0, False), ('c2', 1.0, False)]
There are many ways to choose the starting location; in this case the value of 10 is picked, as it is “far away” from the original values, but hopefully not so far away that the optimiser will not be able to find the best-fit location.
>>> for p in mdl.pars[:3]:
... p.val = 10
...
Note
Since the parameter object - an instance of the
Parameter class - is being
accessed directly here, the value is changed via the
val attribute,
since it does not contain the same overloaded accessor functionality
that allows the setting of the parameter directly from the model.
The fields of the parameter object are:
>>> print(mdl.pars[1])
val = 10.0
min = -3.40282346639e+38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 10.0
default_min = -3.40282346639e+38
default_max = 3.40282346639e+38
How did the optimiser vary the parameters?¶
It can be instructive to see the “search path” taken by
the optimiser; that is, how the parameter values are changed per
iteration. The fit() method will write
these results to a file if the outfile parameter is set
(clobber is set here
to make sure that any existing file is overwritten):
>>> res3 = f.fit(outfile='fitpath.txt', clobber=True)
>>> print(res3.format())
Method = levmar
Statistic = chi2
Initial fit statistic = 196017
Final fit statistic = 4.01682 at function evaluation 8
Data points = 6
Degrees of freedom = 3
Probability [Q-value] = 0.259653
Reduced statistic = 1.33894
Change in statistic = 196013
polynom1d.c0 -9.38489
polynom1d.c1 -2.41692
polynom1d.c2 0.478273
The output file is a simple ASCII file where the columns give the function evaluation number, the statistic, and the parameter values:
# nfev statistic polynom1d.c0 polynom1d.c1 polynom1d.c2
0.000000e+00 1.960165e+05 1.000000e+01 1.000000e+01 1.000000e+01
1.000000e+00 1.960165e+05 1.000000e+01 1.000000e+01 1.000000e+01
2.000000e+00 1.960176e+05 1.000345e+01 1.000000e+01 1.000000e+01
3.000000e+00 1.960172e+05 1.000000e+01 1.000345e+01 1.000000e+01
4.000000e+00 1.961557e+05 1.000000e+01 1.000000e+01 1.000345e+01
5.000000e+00 4.016822e+00 -9.384890e+00 -2.416915e+00 4.782733e-01
6.000000e+00 4.016824e+00 -9.381649e+00 -2.416915e+00 4.782733e-01
7.000000e+00 4.016833e+00 -9.384890e+00 -2.416081e+00 4.782733e-01
8.000000e+00 4.016879e+00 -9.384890e+00 -2.416915e+00 4.784385e-01
9.000000e+00 4.016822e+00 -9.384890e+00 -2.416915e+00 4.782733e-01
As can be seen, a solution is found quickly in this situation. Is it the same as the previous attempt? Although the final statistic values are not the same, they are very close:
>>> res3.statval == res2.statval
False
>>> res3.statval - res2.statval
1.7763568394002505e-15
as are the parameter values:
>>> res2.parvals
(-9.384889507344186, -2.416915493735619, 0.4782733426100644)
>>> res3.parvals
(-9.384889507268973, -2.4169154937357664, 0.47827334260997567)
>>> for p2, p3 in zip(res2.parvals, res3.parvals):
... print("{:+.2e}".format(p3 - p2))
...
+7.52e-11
-1.47e-13
-8.87e-14
The fact that the parameter values are very similar is good evidence for having found the “best fit” location, although this data set was designed to be easy to fit. Real-world examples often require more in-depth analysis.
Comparing models¶
Sherpa contains the calc_mlr() and
calc_ftest() routines which can be
used to compare the model fits (using the change in the
number of degrees of freedom and chi-square statistic) to
see if freeing up \(c_1\) is statistically meaningful.
>>> from sherpa.utils import calc_mlr
>>> calc_mlr(res.dof - res2.dof, res.statval - res2.statval)
5.7788876325820937e-22
This provides evidence that the three-term model (with \(c_1\) free) is a statistically better fit, but the results of these routines should be reviewed carefully to be sure that they are valid; a suggested reference is Protassov et al. 2002, Astrophysical Journal, 571, 545.
Estimating errors¶
Note
Should I add a paragraph mentioning it can be a good idea to
rerun f.fit() to make sure any changes haven’t messsed up
the location?
There are two methods available to estimate errors from the fit object:
Covariance and
Confidence. The method to use can be set
when creating the fit object - using the estmethod parameter - or
after the object has been created, by changing the
estmethod attribute. The default method
is covariance
>>> print(f.estmethod.name)
covariance
which can be significantly faster faster, but less robust, than the confidence technique shown below.
The est_errors() method is used to calculate
the errors on the parameters and returns a
ErrorEstResults object:
>>> coverrs = f.est_errors()
>>> print(coverrs.format())
Confidence Method = covariance
Iterative Fit Method = None
Fitting Method = levmar
Statistic = chi2gehrels
covariance 1-sigma (68.2689%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
polynom1d.c0 -9.38489 -2.91751 2.91751
polynom1d.c1 -2.41692 -0.250889 0.250889
polynom1d.c2 0.478273 -0.0312677 0.0312677
The EstErrResults object can also be displayed directly
(in a similar manner to the FitResults object returned
by fit):
>>> print(coverrs)
datasets = None
methodname = covariance
iterfitname = none
fitname = levmar
statname = chi2gehrels
sigma = 1
percent = 68.2689492137
parnames = ('polynom1d.c0', 'polynom1d.c1', 'polynom1d.c2')
parvals = (-9.3848895072689729, -2.4169154937357664, 0.47827334260997567)
parmins = (-2.9175079401565722, -0.25088931712955043,-0.031267664298717336)
parmaxes = (2.9175079401565722, 0.25088931712955043, 0.031267664298717336)
nfits = 0
The default is to calculate the one-sigma (68.3%) limits for
each thawed parameter. The
error range can be changed
with the
sigma attribute of the
error estimator, and the
set of parameters used in the
analysis can be changed with the parlist
attribute of the
est_errors()
call.
Warning
The error estimate should only be run when at a “good” fit. The
assumption is that the search statistic is chi-square distributed
so the change in statistic as a statistic varies can be related to
a confidence bound. One requirement is that - for chi-square
statistics - is that the reduced statistic value is small
enough. See the max_rstat field of the
EstMethod object.
Give the error if this does not happen (e.g. if c1 has not been fit)?
Changing the error bounds¶
The default is to calculate one-sigma errors, since:
>>> f.estmethod.sigma
1
This can be changed, e.g. to calculate 90% errors which are approximately \(\sigma = 1.6\):
>>> f.estmethod.sigma = 1.6
>>> coverrs90 = f.est_errors()
>>> print(coverrs90.format())
Confidence Method = covariance
Iterative Fit Method = None
Fitting Method = levmar
Statistic = chi2gehrels
covariance 1.6-sigma (89.0401%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
polynom1d.c0 -9.38489 -4.66801 4.66801
polynom1d.c1 -2.41692 -0.401423 0.401423
polynom1d.c2 0.478273 -0.0500283 0.0500283
Note
As can be seen, 1.6 sigma isn’t quite 90%
>>> coverrs90.percent
89.0401416600884
Accessing the covariance matrix¶
The errors created by Covariance provides
access to the covariance matrix via the
extra_output attribute:
>>> print(coverrs.extra_output)
[[ 8.51185258e+00 -4.39950074e-02 -6.51777887e-02]
[ -4.39950074e-02 6.29454494e-02 6.59925111e-04]
[ -6.51777887e-02 6.59925111e-04 9.77666831e-04]]
The order of these rows and columns matches that of the
parnames field:
>>> print([p.split('.')[1] for p in coverrs.parnames])
['c0', 'c1', 'c2']
Changing the estimator¶
Note
There appears to be a
bug in Sherpa
whereby calling est_errors with the Confidence
method here will fail with the message:
IndexError: tuple index out of range
If you see this message then it appears that the thaw_indices
attribute of the fit object has got out of synchronization and
needs to be reset with:
>>> f.thaw_indices = tuple(i for i, p in enumerate(f.model.pars)
... if not p.frozen)
It should not be needed in most cases (it is only needed here because of the earlier change to thaw c1).
The Confidence method takes each
thawed parameter and varies it until it finds the point where the
statistic has increased enough (the value is determined by the
sigma parameter and assumes the statistic is chi-squared
distributed). This is repeated separately for the upper and
lower bounds for each parameter. Since this can take time for
complicated fits, the default behavior is to display a message
when each limit is found (the
order these messages are displayed can change
from run to run):
>>> from sherpa.estmethods import Confidence
>>> f.estmethod = Confidence()
>>> conferrs = f.est_errors()
polynom1d.c0 lower bound: -2.91751
polynom1d.c1 lower bound: -0.250889
polynom1d.c0 upper bound: 2.91751
polynom1d.c2 lower bound: -0.0312677
polynom1d.c1 upper bound: 0.250889
polynom1d.c2 upper bound: 0.0312677
As with the covariance run, the
return value is a ErrorEstResults object:
>>> print(conferrs.format())
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
----- -------- ----------- -----------
polynom1d.c0 -9.38489 -2.91751 2.91751
polynom1d.c1 -2.41692 -0.250889 0.250889
polynom1d.c2 0.478273 -0.0312677 0.0312677
Unlike the covariance errors, these are not guaranteed to be symmetrical (although in this example they are).
Using multiple cores¶
The example data and model used here does not require many
iterations to fit and calculate errors. However, for real-world
situations the error analysis can often take significantly-longer
than the fit. When using the
Confidence technique, the default
is to use all available CPU cores on the machine (the error range
for each parameter can be thought of as a separate task, and the
Python multiprocessing module is used to evaluate these tasks).
This is why the order of the
screen output from the
est_errors() call can vary.
The
parallel
attribute of the error estimator controls whether the
calculation should be run in parallel (True) or
not, and the
numcores attribute
determines how many CPU cores are used (the default is
to use all available cores).
>>> f.estmethod.name
'confidence'
>>> f.estmethod.parallel
True
The Covariance technique does not
provide a parallel option.
Using a subset of parameters¶
To save time, the error calculation can be restricted to a subset
of parameters by setting the parlist parameter of the
est_errors() call. As an example, the
errors on just \(c_1\) can be evaluated with the following:
>>> c1errs = f.est_errors(parlist=(mdl.c1, ))
polynom1d.c1 lower bound: -0.250889
polynom1d.c1 upper bound: 0.250889
>>> print(c1errs)
datasets = None
methodname = confidence
iterfitname = none
fitname = levmar
statname = chi2gehrels
sigma = 1
percent = 68.2689492137
parnames = ('polynom1d.c1',)
parvals = (-2.4169154937357664,)
parmins = (-0.25088931712955054,)
parmaxes = (0.25088931712955054,)
nfits = 2
Visualizing the errors¶
The IntervalProjection class is used to
show how the statistic varies with a single parameter. It steps
through a range of values for a single parameter, fitting the
other thawed parameters, and displays the statistic value
(this gives an indication if the assumptions used to
calculate the errors
are valid):
>>> from sherpa.plot import IntervalProjection
>>> iproj = IntervalProjection()
>>> iproj.calc(f, mdl.c1)
>>> iproj.plot()
The default options display around one sigma from the best-fit
location (the range is estimated from the covariance error estimate).
The prepare() method is used
to change these defaults - e.g. by explicitly setting the min
and max values to use - or, as shown below, by scaling the
covariance error estimate using the fac argument:
>>> iproj.prepare(fac=5, nloop=51)
The number of points has also been increased (the nloop argument)
to keep the curve smooth. Before re-plotting, the
calc()
method needs to be re-run to calculate the new data.
The one-sigma range is also added as vertical dotted
lines using
vline():
>>> iproj.calc(f, mdl.c1)
>>> iproj.plot()
>>> pmin = c1errs.parvals[0] + c1errs.parmins[0]
>>> pmax = c1errs.parvals[0] + c1errs.parmaxes[0]
>>> iproj.vline(pmin, overplot=True, linestyle='dot')
>>> iproj.vline(pmax, overplot=True, linestyle='dot')
Note
The vline()
method is a simple wrapper routine. Direct calls to
matplotlib can also be used to annotate the
visualization.
The RegionProjection class allows the
correlation between two parameters to be shown as a contour plot.
It follows a similar approach as
IntervalProjection, although the
final visualization is created with a call to
contour() rather than
plot:
>>> from sherpa.plot import RegionProjection
>>> rproj = RegionProjection()
>>> rproj.calc(f, mdl.c0, mdl.c2)
>>> rproj.contour()
The limits can be changed, either using the
fac parameter (as shown in the
interval-projection case), or
with the min and max parameters:
>>> rproj.prepare(min=(-22, 0.35), max=(3, 0.6), nloop=(41, 41))
>>> rproj.calc(f, mdl.c0, mdl.c2)
>>> rproj.contour()
>>> xlo, xhi = plt.xlim()
>>> ylo, yhi = plt.ylim()
>>> def get_limits(i):
... return conferrs.parvals[i] + \
... np.asarray([conferrs.parmins[i],
... conferrs.parmaxes[i]])
...
>>> plt.vlines(get_limits(0), ymin=ylo, ymax=yhi)
>>> plt.hlines(get_limits(2), xmin=xlo, xmax=xhi)
The vertical lines are added to indicate the one-sigma errors calculated by the confidence method earlier.
The calc call
sets the fields of the
RegionProjection object with the
data needed to create the plot. In the following example the
data is used to show the search statistic as an image, with the
contours overplotted. First, the axes of the image
are calculated (the extent argument to matplotlib’s
imshow command requires the pixel edges, not the
centers):
>>> xmin, xmax = rproj.x0.min(), rproj.x0.max()
>>> ymin, ymax = rproj.x1.min(), rproj.x1.max()
>>> nx, ny = rproj.nloop
>>> hx = 0.5 * (xmax - xmin) / (nx - 1)
>>> hy = 0.5 * (ymax - ymin) / (ny - 1)
>>> extent = (xmin - hx, xmax + hx, ymin - hy, ymax + hy)
The statistic data is stored as a one-dimensional array, and needs to be reshaped before it can be displayed:
>>> y = rproj.y.reshape((ny, nx))
>>> plt.clf()
>>> plt.imshow(y, origin='lower', extent=extent, aspect='auto', cmap='viridis_r')
>>> plt.colorbar()
>>> plt.xlabel(rproj.xlabel)
>>> plt.ylabel(rproj.ylabel)
>>> rproj.contour(overplot=True)
Simultaneous fits¶
Sherpa uses the
DataSimulFit and
SimulFitModel
classes to fit multiple data seta and models simultaneously.
This can be done explicitly, by combining individual data sets
and models into instances of these classes, or implicitly
with the simulfit() method.
It only makes sense to perform a simultaneous fit if there is some “link” between the various data sets or models, such as sharing one or model components or linking model parameters.
Poisson stats¶
Should there be something about using Poisson stats?
Reference/API¶
The sherpa.fit module¶
Classes
Fit(data, model[, stat, method, estmethod, …]) |
Fit a model to a data set. |
IterFit(data, model, stat, method[, …]) |
|
FitResults(fit, results, init_stat, …) |
The results of a fit. |
StatInfoResults(statname, statval, …[, …]) |
A summary of the current statistic value for one or more data sets. |
ErrorEstResults(fit, results[, parlist]) |
The results of an error estimation run. |
Class Inheritance Diagram¶
The sherpa.estmethods module¶
Classes
EstMethod(name, estfunc) |
|
Confidence([name]) |
|
Covariance([name]) |
|
Projection([name]) |
|
EstNewMin |
Reached a new minimum fit statistic |
Class Inheritance Diagram¶
Visualisation¶
Overview¶
Sherpa has support for different plot backends, at present limited to matplotlib and ChIPS. Interactive visualizations of images is provided by DS9 - an Astronomical image viewer - if installed, whilst there is limited support for visualizing two-dimensional data sets with matplotlib. The classes described in this document do not need to be used, since the data can be plotted directly, but they do provide some conveniences.
The basic approach to creating a visualization using these classes is:
- create an instance of the relevant class (e.g.
DataPlot);- send it the necessary data with the
prepare()method (optional);- perform any necessary calculation with the
calc()method (optional);- and plot the data with the
plot()orcontour()methods (or theoverplot(), andovercontour()variants).
Note
The sherpa.plot module also includes error-estimation routines, such as the IntervalProjection class. This is mixing analysis with visualization, which may not be ideal.
Example¶
>>> import numpy as np
>>> edges = np.asarray([-10, -5, 5, 12, 17, 20, 30, 56, 60])
>>> y = np.asarray([28, 62, 17, 4, 2, 4, 55, 125])
>>> from sherpa.data import Data1DInt
>>> d = Data1DInt('example histogram', edges[:-1], edges[1:], y)
>>> from sherpa.plot import DataPlot
>>> dplot = DataPlot()
>>> dplot.prepare(d)
>>> dplot.plot()
Some text.
>>> from sherpa.plot import Histogram
>>> hplot = Histogram()
>>> hplot.overplot(d.xlo, d.xhi, d.y)
Some more text.
>>> from sherpa.models.basic import Const1D, Gauss1D
>>> mdl = Const1D('base') - Gauss1D('line')
>>> for p, v in zip(mdl.pars, [10, 25, 22, 10]):
... p.val = v
>>> from sherpa.plot import ModelPlot
>>> mplot = ModelPlot()
>>> mplot.prepare(d, mdl)
>>> mplot.plot()
>>> dplot.overplot()
Blah.
>>> from sherpa.optmethods import NelderMead
>>> from sherpa.stats import Cash
>>> from sherpa.fit import Fit
>>> f = Fit(d, mdl, stat=Cash(), method=NelderMead())
>>> f.fit()
>>> from sherpa.plot import IntervalProjection
>>> iproj = IntervalProjection()
>>> iproj.calc(f, mdl.pars[2])
WARNING: hard minimum hit for parameter base.c0
WARNING: hard maximum hit for parameter base.c0
WARNING: hard minimum hit for parameter line.fwhm
WARNING: hard maximum hit for parameter line.fwhm
WARNING: hard minimum hit for parameter line.ampl
WARNING: hard maximum hit for parameter line.ampl
>>> iproj.plot()
Hmmmm. Not good results. Need to reevaluate the data beng used here.
Reference/API¶
The sherpa.plot module¶
A visualization interface to Sherpa
Classes
Confidence1D() |
|
IntervalProjection() |
|
IntervalUncertainty() |
|
Confidence2D() |
|
RegionProjection() |
|
RegionUncertainty() |
|
Contour() |
Base class for contour plots |
DataContour() |
Create contours of 2D data. |
PSFContour() |
Derived class for creating 2D PSF contours |
FitContour() |
Derived class for creating 2D combination data and model contours |
ModelContour() |
Derived class for creating 2D model contours |
RatioContour() |
Derived class for creating 2D ratio contours (data:model) |
ResidContour() |
Derived class for creating 2D residual contours (data-model) |
SourceContour() |
Derived class for creating 2D model contours |
Histogram() |
Base class for histogram plots |
Plot() |
Base class for line plots |
DataPlot() |
Create 1D plots of data values. |
PSFPlot() |
Derived class for creating 1D PSF kernel data plots |
FitPlot() |
Derived class for creating 1D combination data and model plots |
ModelPlot() |
Create 1D plots of model values. |
ChisqrPlot() |
Create plots of the chi-square value per point. |
ComponentModelPlot() |
|
DelchiPlot() |
Create plots of the delta-chi value per point. |
RatioPlot() |
Create plots of the ratio of data to model per point. |
ResidPlot() |
Create plots of the residuals (data - model) per point. |
SourcePlot() |
Create 1D plots of unconcolved model values. |
ComponentSourcePlot() |
|
SplitPlot([rows, cols]) |
Create multiple plots. |
JointPlot() |
|
Point() |
Base class for point plots |
Class Inheritance Diagram¶
The sherpa.astro.plot module¶
Classes
ChisqrPlot() |
Create plots of the chi-square value per point. |
BkgChisqrPlot() |
Derived class for creating background plots of 1D chi**2 ((data-model)/error)**2 |
DataPlot() |
Create 1D plots of data values. |
BkgDataPlot() |
Derived class for creating plots of background counts |
DelchiPlot() |
Create plots of the delta-chi value per point. |
BkgDelchiPlot() |
Derived class for creating background plots of 1D delchi chi ((data-model)/error) |
FitPlot() |
Derived class for creating 1D combination data and model plots |
BkgFitPlot() |
Derived class for creating plots of background counts with fitted model |
HistogramPlot() |
|
ARFPlot() |
Create plots of the ancillary response file (ARF). |
ModelHistogram() |
Derived class for creating 1D PHA model histogram plots |
BkgModelHistogram() |
Derived class for creating 1D background PHA model histogram plots |
OrderPlot() |
Derived class for creating plots of the convolved source model using selected multiple responses |
SourcePlot() |
Create 1D plots of unconcolved model values. |
BkgSourcePlot() |
Derived class for plotting the background unconvolved source model |
ModelPlot() |
Create 1D plots of model values. |
BkgModelPlot() |
Derived class for creating plots of background model |
RatioPlot() |
Create plots of the ratio of data to model per point. |
BkgRatioPlot() |
Derived class for creating background plots of 1D ratio (data:model) |
ResidPlot() |
Create plots of the residuals (data - model) per point. |
BkgResidPlot() |
Derived class for creating background plots of 1D residual (data-model) |
Class Inheritance Diagram¶
Markov Chain Monte Carlo and Poisson data¶
Sherpa provides a Markov Chain Monte Carlo (MCMC) method designed for Poisson-distributed data. It was originally developed as the Bayesian Low-Count X-ray Spectral (BLoCXS) package, but has since been incorporated into Sherpa. It is developed from the work presented in Analysis of Energy Spectra with Low Photon Counts via Bayesian Posterior Simulation by van Dyk et al.
Unlike many MCMC implementations, idea is that have some idea of the search surface at the optimum - i.e. the covariance matrix - and then use that to explore this region.
Example¶
Note
This example probably needs to be simplified to reduce the run time
Simulate the data¶
Create a simulated data set:
>>> np.random.seed(2)
>>> x0low, x0high = 3000, 4000
>>> x1low, x1high = 4000, 4800
>>> dx = 15
>>> x1, x0 = np.mgrid[x1low:x1high:dx, x0low:x0high:dx]
Convert to 1D arrays:
>>> shape = x0.shape
>>> x0, x1 = x0.flatten(), x1.flatten()
Create the model used to simulate the data:
>>> from sherpa.astro.models import Beta2D
>>> truth = Beta2D()
>>> truth.xpos, truth.ypos = 3512, 4418
>>> truth.r0, truth.alpha = 120, 2.1
>>> truth.ampl = 12
Evaluate the model to calculate the expected values:
>>> mexp = truth(x0, x1).reshape(shape)
Create the simulated data by adding in Poisson-distributed noise:
>>> msim = np.random.poisson(mexp)
What does the data look like?¶
Use an arcsinh transform to view the data, based on the work of Lupton, Gunn & Szalay (1999).
>>> plt.imshow(np.arcsinh(msim), origin='lower', cmap='viridis',
... extent=(x0low, x0high, x1low, x1high),
... interpolation='nearest', aspect='auto')
>>> plt.title('Simulated image')
Find the starting point for the MCMC¶
Set up a model and use the standard Sherpa approach to find a good starting place for the MCMC analysis:
>>> from sherpa import data, stats, optmethods, fit
>>> d = data.Data2D('sim', x0, x1, msim.flatten(), shape=shape)
>>> mdl = Beta2D()
>>> mdl.xpos, mdl.ypos = 3500, 4400
Use a Likelihood statistic and Nelder-Mead algorithm:
>>> f = fit.Fit(d, mdl, stats.Cash(), optmethods.NelderMead())
>>> res = f.fit()
>>> print(res.format())
Method = neldermead
Statistic = cash
Initial fit statistic = 20048.5
Final fit statistic = 607.229 at function evaluation 777
Data points = 3618
Degrees of freedom = 3613
Change in statistic = 19441.3
beta2d.r0 121.945
beta2d.xpos 3511.99
beta2d.ypos 4419.72
beta2d.ampl 12.0598
beta2d.alpha 2.13319
Now calculate the covariance matrix (the default error estimate):
>>> f.estmethod
<Covariance error-estimation method instance 'covariance'>
>>> eres = f.est_errors()
>>> print(eres.format())
Confidence Method = covariance
Iterative Fit Method = None
Fitting Method = neldermead
Statistic = cash
covariance 1-sigma (68.2689%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
beta2d.r0 121.945 -7.12579 7.12579
beta2d.xpos 3511.99 -2.09145 2.09145
beta2d.ypos 4419.72 -2.10775 2.10775
beta2d.ampl 12.0598 -0.610294 0.610294
beta2d.alpha 2.13319 -0.101558 0.101558
The covariance matrix is stored in the extra_output attribute:
>>> cmatrix = eres.extra_output
>>> pnames = [p.split('.')[1] for p in eres.parnames]
>>> plt.imshow(cmatrix, interpolation='nearest', cmap='viridis')
>>> plt.xticks(np.arange(5), pnames)
>>> plt.yticks(np.arange(5), pnames)
>>> plt.colorbar()
Run the chain¶
Finally, run a chain (use a small number to keep the run time low for this example):
>>> from sherpa.sim import MCMC
>>> mcmc = MCMC()
>>> mcmc.get_sampler_name()
>>> draws = mcmc.get_draws(f, cmatrix, niter=1000)
>>> svals, accept, pvals = draws
>>> pvals.shape
(5, 1001)
>>> accept.sum() * 1.0 / 1000
0.48499999999999999
Trace plots¶
>>> plt.plot(pvals[0, :])
>>> plt.xlabel('Iteration')
>>> plt.ylabel('r0')
Or using the sherpa.plot module:
>>> from sherpa import plot
>>> tplot = plot.TracePlot()
>>> tplot.prepare(svals, name='Statistic')
>>> tplot.plot()
PDF of a parameter¶
>>> pdf = plot.PDFPlot()
>>> pdf.prepare(pvals[1, :], 20, False, 'xpos', name='example')
>>> pdf.plot()
Add in the covariance estimate:
>>> xlo, xhi = eres.parmins[1] + eres.parvals[1], eres.parmaxes[1] + eres.parvals[1]
>>> plt.annotate('', (xlo, 90), (xhi, 90), arrowprops={'arrowstyle': '<->'})
>>> plt.plot([eres.parvals[1]], [90], 'ok')
CDF for a parameter¶
Normalise by the actual answer to make it easier to see how well the results match reality:
>>> cdf = plot.CDFPlot()
>>> plt.subplot(2, 1, 1)
>>> cdf.prepare(pvals[1, :] - truth.xpos.val, r'$\Delta x$')
>>> cdf.plot(clearwindow=False)
>>> plt.title('')
>>> plt.subplot(2, 1, 2)
>>> cdf.prepare(pvals[2, :] - truth.ypos.val, r'$\Delta y$')
>>> cdf.plot(clearwindow=False)
>>> plt.title('')
Scatter plot¶
>>> plt.scatter(pvals[0, :] - truth.r0.val,
... pvals[4, :] - truth.alpha.val, alpha=0.3)
>>> plt.xlabel(r'$\Delta r_0$', size=18)
>>> plt.ylabel(r'$\Delta \alpha$', size=18)
This can be compared to the
RegionProjection calculation:
>>> plt.scatter(pvals[0, :], pvals[4, :], alpha=0.3)
>>> from sherpa.plot import RegionProjection
>>> rproj = RegionProjection()
>>> rproj.prepare(min=[95, 1.8], max=[150, 2.6], nloop=[21, 21])
>>> rproj.calc(f, mdl.r0, mdl.alpha)
>>> rproj.contour(overplot=True)
>>> plt.xlabel(r'$r_0$'); plt.ylabel(r'$\alpha$')
Reference/API¶
The sherpa.sim module¶
Monte-Carlo Markov Chain support for low-count data (Poisson statistics).
The sherpa.sim module provides support for exploring the posterior
probability density of parameters in a fit to low-count data, for
which Poisson statistics hold, using a Bayesian algorithm and a
Monte-Carlo Markov Chain (MCMC). It was originally known as the
pyBLoCXS (python Bayesian Low-Count X-ray Spectral) package [1], but
has since been incorporated into Sherpa.
The Sherpa UI modules - e.g. sherpa.ui and sherpa.astro.ui - provide
many of the routines described below (e.g. list_samplers).
Acknowledgements¶
The original version of the code was developed by the CHASC Astro-Statistics collaboration http://hea-www.harvard.edu/AstroStat/, and was called pyBLoCXS. It has since been developed by the Chandra X-ray Center and weas added to Sherpa in version 4.5.1.
Overview¶
The algorithm explores parameter space at a suspected minimum - i.e. after a standard Sherpa fit. It supports a flexible definition of priors and allows for variations in the calibration information. It can be used to compute posterior predictive p-values for the likelihood ratio test [2]. Future versions will allow for the incorporation of calibration uncertainty [3].
MCMC is a complex computational technique that requires some sophistication on the part of its users to ensure that it both converges and explores the posterior distribution properly. The pyBLoCXS code has been tested with a number of simple single-component spectral models. It should be used with great care in more complex settings. The code is based on the methods in [4] but employs a different MCMC sampler than is described in that article. A general description of the techniques employed along with their convergence diagnostics can be found in the Appendices of [4] and in [5].
Jumping Rules¶
The jumping rule determines how each step in the Monte-Carlo Markov
Chain is calculated. The setting can be changed using set_sampler.
The sherpa.sim module provides the following rules, which may
be augmented by other modules:
MHuses a Metropolis-Hastings jumping rule that is a multivariate t-distribution with user-specified degrees of freedom centered on the best-fit parameters, and with multivariate scale determined by thecovarfunction applied at the best-fit location.MetropolisMHmixes this Metropolis Hastings jumping rule with a Metropolis jumping rule centered at the current draw, in both cases drawing from the same t-distribution as used withMH. The probability of using the best-fit location as the start of the jump is given by thep_Mparameter of the rule (useget_samplerorget_sampler_optto view andset_sampler_optto set this value), otherwise the jump is from the previous location in the chain.
Options for the sampler are retrieved and set by get_sampler or
get_sampler_opt, and set_sampler_opt respectively. The list of
available samplers is given by list_samplers.
Choosing the parameter values¶
By default, the prior on each parameter is taken to be flat, varying
from the parameter minima to maxima values. This prior can be changed
using the set_prior function, which can set the prior for a
parameter to a function or Sherpa model. The list of currently set
prior-parameter pairs is returned by the list_priors function, and the
prior function associated with a particular Sherpa model parameter may be
accessed with get_prior.
Running the chain¶
The get_draws function runs a pyBLoCXS chain using fit information
associated with the specified data set(s), and the currently set sampler and
parameter priors, for a specified number of iterations. It returns an array of
statistic values, an array of acceptance Booleans, and a 2-D array of
associated parameter values.
Analyzing the results¶
The module contains several routines to visualize the results of the chain,
including plot_trace, plot_cdf, and plot_pdf, along with
sherpa.utils.get_error_estimates for calculating the limits from a
parameter chain.
References
| [1] | http://hea-www.harvard.edu/AstroStat/pyBLoCXS/ |
| [2] | “Statistics, Handle with Care: Detecting Multiple Model Components with the Likelihood Ratio Test”, Protassov et al., 2002, ApJ, 571, 545 http://adsabs.harvard.edu/abs/2002ApJ…571..545P |
| [3] | “Accounting for Calibration Uncertainties in X-ray Analysis: Effective Areas in Spectral Fitting”, Lee et al., 2011, ApJ, 731, 126 http://adsabs.harvard.edu/abs/2011ApJ…731..126L |
| [4] | (1, 2) “Analysis of Energy Spectra with Low Photon Counts via Bayesian Posterior Simulation”, van Dyk et al. 2001, ApJ, 548, 224 http://adsabs.harvard.edu/abs/2001ApJ…548..224V |
| [5] | Chapter 11 of Gelman, Carlin, Stern, and Rubin (Bayesian Data Analysis, 2nd Edition, 2004, Chapman & Hall/CRC). |
Examples
Analysis proceeds as normal, up to the point that a good fit has
been determined, as shown below (note that a Poisson likelihood,
such as the cash statistic, must be used):
>>> from sherpa.astro import ui
>>> ui.load_pha('src.pi')
>>> ui.notice(0.5, 7)
>>> ui.set_source(ui.xsphabs.gal * ui.xspowerlaw.pl)
>>> ui.set_stat('cash')
>>> ui.set_method('simplex')
>>> ui.fit()
>>> ui.covar()
Once the best-fit location has been determined (which may require
multiple calls to fit), the chain can be run. In this example
the default sampler (MetropolisMH) and default parameter priors
(flat, varying between the minimum and maximum values) are used,
as well as the default number of iterations (1000):
>>> stats, accept, params = ui.get_draws()
The stats array contains the fit statistic for each iteration
(the first element of these arrays is the starting point of the chain,
so there will be 1001 elements in this case). The “trace” - i.e.
statistic versus iteration - can be plotted using:
>>> ui.plot_trace(stats, name='stat')
The accept array indicates whether, at each iteration, the proposed
jump was accepted, (True) or if the previous iterations parameter
values are used. This can be used to look at the acceptance rate for
the chain (dropping the last element and a burn-in period, which
here is arbitrarily taken to be 100):
>>> nburn = 100
>>> arate = accept[nburn:-1].sum() * 1.0 / (len(accept) - nburn - 1)
>>> print("acceptance rate = {}".format(arate))
The trace of the parameter values can also be displayed; in this example a burn-in period has not been removed):
>>> par1 = params[:, 0]
>>> par2 = params[:, 1]
>>> ui.plot_trace(par1, name='par1')
>>> ui.plot_trace(par2, name='par2')
The cumulative distribution can also be viewed:
>>> ui.plot_cdf(par1[nburn:], name='par1')
as well as the probability density:
>>> ui.plot_[pdf(par2[nburn:], name='par2')
The traces can be used to estimate the credible interval for a parameter:
>>> from sherpa.utils import get_error_estimates
>>> pval, plo, phi = get_error_estimates(par1[nburn:])
Classes
MCMC() |
High-level UI to pyBLoCXS that joins the loop in ‘Walk’ with the jumping rule in ‘Sampler’. |
Functions
flat(x) |
The flat prior (returns 1 everywhere). |
inverse(x) |
Returns the inverse of x. |
inverse2(x) |
Returns the invers of x^2. |
Class Inheritance Diagram¶
The sherpa.sim.mh module¶
pyBLoCXS is a sophisticated Markov chain Monte Carlo (MCMC) based algorithm designed to carry out Bayesian Low-Count X-ray Spectral (BLoCXS) analysis in the Sherpa environment. The code is a Python extension to Sherpa that explores parameter space at a suspected minimum using a predefined Sherpa model to high-energy X-ray spectral data. pyBLoCXS includes a flexible definition of priors and allows for variations in the calibration information. It can be used to compute posterior predictive p-values for the likelihood ratio test (see Protassov et al., 2002, ApJ, 571, 545). Future versions will allow for the incorporation of calibration uncertainty (Lee et al., 2011, ApJ, 731, 126).
MCMC is a complex computational technique that requires some sophistication on the part of its users to ensure that it both converges and explores the posterior distribution properly. The pyBLoCXS code has been tested with a number of simple single-component spectral models. It should be used with great care in more complex settings. Readers interested in Bayesian low-count spectral analysis should consult van Dyk et al. (2001, ApJ, 548, 224). pyBLoCXS is based on the methods in van Dyk et al. (2001) but employs a different MCMC sampler than is described in that article. In particular, pyBLoCXS has two sampling modules. The first uses a Metropolis-Hastings jumping rule that is a multivariate t-distribution with user specified degrees of freedom centered on the best spectral fit and with multivariate scale determined by the Sherpa function, covar(), applied to the best fit. The second module mixes this Metropolis Hastings jumping rule with a Metropolis jumping rule centered at the current draw, also sampling according to a t-distribution with user specified degrees of freedom and multivariate scale determined by a user specified scalar multiple of covar() applied to the best fit.
A general description of the MCMC techniques we employ along with their convergence diagnostics can be found in Appendices A.2 - A.4 of van Dyk et al. (2001) and in more detail in Chapter 11 of Gelman, Carlin, Stern, and Rubin (Bayesian Data Analysis, 2nd Edition, 2004, Chapman & Hall/CRC).
http://hea-www.harvard.edu/AstroStat/pyBLoCXS/
Classes
LimitError |
|
MH(fcn, sigma, mu, dof, *args) |
The Metropolis Hastings Sampler |
MetropolisMH(fcn, sigma, mu, dof, *args) |
The Metropolis Metropolis-Hastings Sampler |
Sampler() |
|
Walk([sampler, niter]) |
Functions
dmvnorm(x, mu, sigma[, log]) |
Probability Density of a multi-variate Normal distribution |
dmvt(x, mu, sigma, dof[, log, norm]) |
Probability Density of a multi-variate Student’s t distribution |
Class Inheritance Diagram¶
The sherpa.sim.sample module¶
Classes
Functions
multivariate_t(mean, cov, df[, size]) |
Draw random deviates from a multivariate Student’s T distribution Such a distribution is specified by its mean covariance matrix, and degrees of freedom. |
multivariate_cauchy(mean, cov[, size]) |
This needs to be checked too! A reference to the literature the better |
normal_sample(fit[, num, sigma, correlate, …]) |
Sample the fit statistic by taking the parameter values from a normal distribution. |
uniform_sample(fit[, num, factor, numcores]) |
Sample the fit statistic by taking the parameter values from an uniform distribution. |
t_sample(fit[, num, dof, numcores]) |
Sample the fit statistic by taking the parameter values from a Student’s t-distribution. |
Class Inheritance Diagram¶
The sherpa.sim.simulate module¶
Classes for PPP simulations
Classes
LikelihoodRatioResults(ratios, stats, …) |
The results of a likelihood ratio comparison simulation. |
LikelihoodRatioTest() |
Likelihood Ratio Test. |
Class Inheritance Diagram¶
Utility routines¶
There are a number of utility routines provided by Sherpa that may be useful. Should they be documented here or elsewhere?
Unfortunately it is not always obvious whether a routine is for use with the Object-Oriented API or the Session API.
Reference/API¶
The sherpa.utils.err module¶
Sherpa specific exceptions
Exceptions
SherpaErr(dict, *args) |
Base class for all Sherpa exceptions |
ArgumentErr(key, *args) |
|
ArgumentTypeErr(key, *args) |
|
ConfidenceErr(key, *args) |
|
DS9Err(key, *args) |
|
DataErr(key, *args) |
Error in creating or using a data set |
EstErr(key, *args) |
|
FitErr(key, *args) |
|
IOErr(key, *args) |
|
IdentifierErr(key, *args) |
|
ImportErr(key, *args) |
|
InstrumentErr(key, *args) |
|
ModelErr(key, *args) |
Error in creating or using a model |
NotImplementedErr(key, *args) |
|
PSFErr(key, *args) |
|
ParameterErr(key, *args) |
Error in creating or using a model |
PlotErr(key, *args) |
Error in creating or using a plotting class |
SessionErr(key, *args) |
|
StatErr(key, *args) |
Class Inheritance Diagram¶
The sherpa.utils module¶
Functions
Knuth_close(x, y, tol[, myop]) |
Check whether two floating-point numbers are close together. |
_guess_ampl_scale |
The scaling applied to a value to create its range. |
apache_muller(fcn, xa, xb[, fa, fb, args, …]) |
|
bisection(fcn, xa, xb[, fa, fb, args, …]) |
|
bool_cast(val) |
Convert a string to a boolean. |
calc_ftest(dof1, stat1, dof2, stat2) |
Compare two models using the F test. |
calc_mlr(delta_dof, delta_stat) |
Compare two models using the Maximum Likelihood Ratio test. |
calc_total_error([staterror, syserror]) |
Add statistical and systematic errors in quadrature. |
create_expr(vals[, mask, format, delim]) |
collapse a list of channels into an expression using hyphens and commas to indicate filtered intervals. |
dataspace1d(start, stop[, step, numbins]) |
Populates an integrated grid |
dataspace2d(dim) |
Populates a blank image dataset |
demuller(fcn, xa, xb, xc[, fa, fb, fc, …]) |
A root-finding algorithm using Muller’s method. |
erf(x) |
Calculate the error function. |
export_method(meth[, name, modname]) |
Given a bound instance method, return a simple function that wraps it. |
extract_kernel(kernel, dims_kern, dims_new, …) |
Extract the kernel. |
filter_bins(mins, maxes, axislist) |
What mask represents the given set of filters? |
gamma(z) |
Calculate the Gamma function. |
get_error_estimates(x[, sorted]) |
Compute the median and (-1,+1) sigma values for the data. |
get_fwhm(y, x[, xhi]) |
Estimate the width of the data. |
get_keyword_defaults(func[, skip]) |
Return the keyword arguments and their default values. |
get_keyword_names(func[, skip]) |
Return the names of the keyword arguments. |
get_midpoint(a) |
Estimate the middle of the data. |
get_num_args(func) |
Return the number of arguments for a function. |
get_peak(y, x[, xhi]) |
Estimate the peak position of the data. |
get_position(y, x[, xhi]) |
Get 1D model parameter positions pos (val, min, max) |
get_valley(y, x[, xhi]) |
Estimate the position of the minimum of the data. |
guess_amplitude(y, x[, xhi]) |
Guess model parameter amplitude (val, min, max) |
guess_amplitude2d(y, x0lo, x1lo[, x0hi, x1hi]) |
Guess 2D model parameter amplitude (val, min, max) |
guess_amplitude_at_ref(r, y, x[, xhi]) |
Guess model parameter amplitude (val, min, max) |
guess_bounds(x[, xhi]) |
Guess the bounds of a parameter from the independent axis. |
guess_fwhm(y, x[, xhi, scale]) |
Estimate the value and valid range for the FWHM of the data. |
guess_position(y, x0lo, x1lo[, x0hi, x1hi]) |
Guess 2D model parameter positions xpos, ypos ({val0, min0, max0}, |
guess_radius(x0lo, x1lo[, x0hi, x1hi]) |
Guess the radius parameter of a 2D model. |
guess_reference(pmin, pmax, x[, xhi]) |
Guess model parameter reference (val, min, max) |
histogram1d(x, x_lo, x_hi) |
Create a 1D histogram from a sequence of samples. |
histogram2d(x, y, x_grid, y_grid) |
Create 2D histogram from a sequence of samples. |
igam(a, x) |
Calculate the regularized incomplete Gamma function (lower). |
igamc(a, x) |
Calculate the complement of the regularized incomplete Gamma function (upper). |
incbet(a, b, x) |
Calculate the incomplete Beta function. |
interpolate(xout, xin, yin[, function]) |
One-dimensional interpolation. |
is_binary_file(filename) |
Estimate if a file is a binary file. |
lgam(z) |
Calculate the log (base e) of the Gamma function. |
linear_interp(xout, xin, yin) |
Linear one-dimensional interpolation. |
multinormal_pdf(x, mu, sigma) |
The PDF of a multivariate-normal. |
multit_pdf(x, mu, sigma, dof) |
The PDF of a multivariate student-t. |
nearest_interp(xout, xin, yin) |
Nearest-neighbor one-dimensional interpolation. |
neville(xout, xin, yin) |
Polynomial one-dimensional interpolation using Neville’s method. |
neville2d(xinterp, yinterp, x, y, fval) |
Polynomial two-dimensional interpolation using Neville’s method. |
new_muller(fcn, xa, xb[, fa, fb, args, …]) |
|
normalize(xs) |
Normalize an array. |
numpy_convolve(a, b) |
Convolve two 1D arrays together using NumPy’s FFT. |
pad_bounding_box(kernel, mask) |
Expand the kernel to match the mask. |
parallel_map(function, sequence[, numcores]) |
Run a function on a sequence of inputs in parallel. |
param_apply_limits(param_limits, par[, …]) |
Apply the given limits to a parameter. |
parse_expr(expr) |
parse a filter expression into numerical components for notice/ignore e.g. |
poisson_noise(x) |
Draw samples from a Poisson distribution. |
print_fields(names, vals[, converters]) |
Given a list of strings names and mapping vals, where names is a subset of vals.keys(), return a listing of name/value pairs printed one per line in the format ‘<name> = <value>’. |
quantile(sorted_array, f) |
Return the quantile element from sorted_array, where f is [0,1] using linear interpolation. |
rebin(y0, x0lo, x0hi, x1lo, x1hi) |
Rebin a histogram. |
sao_arange(start, stop[, step]) |
Create a range of values between start and stop. |
sao_fcmp(x, y, tol) |
Compare y to x, using an absolute tolerance. |
set_origin(dims[, maxindex]) |
Return the position of the origin of the kernel. |
sum_intervals(src, indx0, indx1) |
Sum up data within one or more pairs of indexes. |
zeroin(fcn, xa, xb[, fa, fb, args, maxfev, tol]) |
Obtain a zero of a function of one variable using Brent’s root finder. |
Classes
NoNewAttributesAfterInit() |
Prevents attribute deletion and setting of new attributes after __init__ has been called. |
The sherpa.utils.testing module¶
Functions
requires_data(*args, **kwargs) |
|
requires_ds9(*args, **kwargs) |
|
requires_fits(*args, **kwargs) |
|
requires_group(*args, **kwargs) |
|
requires_package(*args) |
|
requires_plotting(*args, **kwargs) |
|
requires_pylab(*args, **kwargs) |
|
requires_stk(*args, **kwargs) |
|
requires_xspec(*args, **kwargs) |
Classes
SherpaTestCase([methodName]) |
Base class for Sherpa unit tests. |
The sherpa.astro.io module¶
Functions
read_table(arg[, ncols, colkeys, dstype]) |
Create a dataset from a tabular file. |
read_image(arg[, coord, dstype]) |
Create an image dataset from a file. |
read_arf(arg) |
Create a DataARF object. |
read_rmf(arg) |
Create a DataRMF object. |
read_arrays(*args) |
Create a dataset from multiple arrays. |
read_pha(arg[, use_errors, use_background]) |
Create a DataPHA object. |
write_image(filename, dataset[, ascii, clobber]) |
Write out an image. |
write_pha(filename, dataset[, ascii, clobber]) |
Write out a PHA dataset. |
write_table(filename, dataset[, ascii, clobber]) |
Write out a table. |
pack_table(dataset) |
Convert a Sherpa data object into an I/O item (tabular). |
pack_image(dataset) |
Convert a Sherpa data object into an I/O item (image). |
pack_pha(dataset) |
Convert a Sherpa PHA data object into an I/O item (tabular). |
read_table_blocks(arg[, make_copy]) |
Return the HDU elements (columns and header) from a FITS table. |
Simple Interpolation¶
Overview¶
Although Sherpa allows you to fit complex models to complex data sets with complex statistics and complex optimisers, it can also be used for simple situations, such as interpolating a function. In this example a one-dimensional set of data is given - i.e. \((x_i, y_i)\) - and a polynomial of order 2 is fit to the data. The model is then used to interpolate (and extrapolate) values. The walk through ends with changing the fit to a linear model (i.e. polynomial of order 1) and a comparison of the two model fits.
Setting up¶
The following sections will load in classes from Sherpa as needed, but it is assumed that the following modules have been loaded:
>>> import numpy as np
>>> import matplotlib.pyplot as plt
Loading the data¶
The data is the following:
| x | y |
|---|---|
| 1 | 1 |
| 1.5 | 1.5 |
| 2 | 1.75 |
| 4 | 3.25 |
| 8 | 6 |
| 17 | 16 |
which can be “loaded” into Sherpa using the
Data1D class:
>>> from sherpa.data import Data1D
>>> x = [1, 1.5, 2, 4, 8, 17]
>>> y = [1, 1.5, 1.75, 3.25, 6, 16]
>>> d = Data1D('interpolation', x, y)
>>> print(d)
name = interpolation
x = [1, 1.5, 2, 4, 8, 17]
y = [1, 1.5, 1.75, 3.25, 6, 16]
staterror = None
syserror = None
None
This can be displayed using the DataPlot class:
>>> from sherpa.plot import DataPlot
>>> dplot = DataPlot()
>>> dplot.prepare(d)
>>> dplot.plot()
Setting up the model¶
For this example, a second-order polynomial is going to be fit to the
data by using the Polynom1D class:
>>> from sherpa.models.basic import Polynom1D
>>> mdl = Polynom1D()
>>> print(mdl)
polynom1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
polynom1d.c0 thawed 1 -3.40282e+38 3.40282e+38
polynom1d.c1 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c2 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c3 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c4 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c5 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c6 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c7 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c8 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.offset frozen 0 -3.40282e+38 3.40282e+38
The help for Polynom1D shows that the model is defined as:
so to get a second-order polynomial we have to
thaw the c2
parameter (the linear term c1 is kept at 0 to show that the
choice of parameter to fit is up to the user):
>>> mdl.c2.thaw()
This model can be compared to the data using the
ModelPlot class (note that, unlike
the data plot, the
prepare() method takes both
the data - needed to know what \(x_i\) to use - and the model):
>>> from sherpa.plot import ModelPlot
>>> mplot = ModelPlot()
>>> mplot.prepare(d, mdl)
>>> dplot.plot()
>>> mplot.overplot()
Since the default parameter values are still being used, the result is not a good description of the data. Let’s fix this!
Fitting the model to the data¶
Since we have no error bars, we are going to use least-squares
minimisation - that is, minimise the square of the distance between
the model and the data using the
LeastSq statisic and the
NelderMead optimiser
(for this case the LevMar optimiser is likely
to produce as good a result but faster, but I have chosen to
select the more robust method):
>>> from sherpa.stats import LeastSq
>>> from sherpa.optmethods import NelderMead
>>> from sherpa.fit import Fit
>>> f = Fit(d, mdl, stat=LeastSq(), method=NelderMead())
>>> print(f)
data = interpolation
model = polynom1d
stat = LeastSq
method = NelderMead
estmethod = Covariance
In this case there is no need to change any of the options for the
optimiser (the least-squares statistic has no options), so the objects
are passed straight to the Fit object.
The fit() method is used to fit the data; as it
returns useful information (in a FitResults
object) we capture this in the res variable, and then check that
the fit was succesfull (i.e. it converged):
>>> res = f.fit()
>>> res.succeeded
True
For this example the time to perform the fit is very short, but for complex data sets and models the call can take a long time!
A quick summary of the fit results is available via the
format() method, while printing the
variable retutrns more details:
>>> print(res.format())
Method = neldermead
Statistic = leastsq
Initial fit statistic = 255.875
Final fit statistic = 2.4374 at function evaluation 264
Data points = 6
Degrees of freedom = 4
Change in statistic = 253.438
polynom1d.c0 1.77498
polynom1d.c2 0.0500999
>>> print(res)
datasets = None
itermethodname = none
methodname = neldermead
statname = leastsq
succeeded = True
parnames = ('polynom1d.c0', 'polynom1d.c2')
parvals = (1.7749826216226083, 0.050099944904353017)
statval = 2.4374045728256455
istatval = 255.875
dstatval = 253.437595427
numpoints = 6
dof = 4
qval = None
rstat = None
message = Optimization terminated successfully
nfev = 264
The best-fit parameter values can also be retrieved from the model itself:
>>> print(mdl)
polynom1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
polynom1d.c0 thawed 1.77498 -3.40282e+38 3.40282e+38
polynom1d.c1 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c2 thawed 0.0500999 -3.40282e+38 3.40282e+38
polynom1d.c3 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c4 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c5 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c6 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c7 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c8 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.offset frozen 0 -3.40282e+38 3.40282e+38
as can the current fit statistic (as this is for fitting a second-order polynomial I’ve chosen to label the variable with a suffix of 2, which will make more sense below):
>>> stat2 = f.calc_stat()
>>> print("Statistic = {:.4f}".format(stat2))
Statistic = 2.4374
Note
In an actual analysis session the fit would probably be repeated, perhaps with a different optimiser, and starting from a different set of parameter values, to give more confidence that the fit has not been caught in a local minimum. This example is simple enough that this is not needed here.
To compare the new model to the data I am going to use a
FitPlot - which combines a DataPlot
and ModelPlot - and a ResidPlot - to look
at the residuals, defined as \({\rm data}_i - {\rm model}_i\),
using the SplitPlot class to orchestrate
the display (note that mplot needs to be re-created since the
model has changed since the last time its prepare method
was called):
>>> from sherpa.plot import FitPlot, ResidPlot, SplitPlot
>>> fplot = FitPlot()
>>> mplot.prepare(d, mdl)
>>> fplot.prepare(dplot, mplot)
>>> splot = SplitPlot()
>>> splot.addplot(fplot)
>>> rplot = ResidPlot()
>>> rplot.prepare(d, mdl, stat=LeastSq())
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with leastsq
>>> rplot.plot_prefs['yerrorbars'] = False
>>> splot.addplot(rplot)
The default behavior for the residual plot is to include error bars,
here calculated using the Chi2XspecVar class,
but they have been turned off - by setting the
yerrorbars option to False - since they are not meaningful here.
Interpolating values¶
The model can be evaluated directly by supplying it with the independent-axis values; for instance for \(x\) equal to 2, 5, and 10:
>>> print(mdl([2, 5, 10]))
[ 1.9753824 3.02748124 6.78497711]
It can also be used to extrapolate the model outside the range of the data (as long as the model is defined for these values):
>>> print(mdl([-100]))
[ 502.77443167]
>>> print(mdl([234.56]))
[ 2758.19347071]
Changing the fit¶
Let’s see how the fit looks if we use a linear model instead. This
means thawing out the c1 parameter and clearing c2:
>>> mdl.c1.thaw()
>>> mdl.c2 = 0
>>> mdl.c2.freeze()
>>> f.fit()
<Fit results instance>
As this is a simple case, I am ignoring the return value from the
fit() method, but in an actual analysis
session it should be checked to ensure the fit converged.
The new model parameters are:
>>> print(mdl)
polynom1d
Param Type Value Min Max Units
----- ---- ----- --- --- -----
polynom1d.c0 thawed -0.248624 -3.40282e+38 3.40282e+38
polynom1d.c1 thawed 0.925127 -3.40282e+38 3.40282e+38
polynom1d.c2 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c3 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c4 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c5 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c6 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c7 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.c8 frozen 0 -3.40282e+38 3.40282e+38
polynom1d.offset frozen 0 -3.40282e+38 3.40282e+38
and the best-fit statistic value can be compared to the earlier version:
>>> stat1 = f.calc_stat()
>>> print("Statistic: order 1 = {:.3f} order 2 = {:.3f}".format(stat1, stat2))
Statistic: order 1 = 1.898 order 2 = 2.437
Note
Sherpa provides several routines for comparing statistic values,
such as sherpa.utils.calc_ftest() and
sherpa.utils.calc_mlr(), to see if one can be preferred
over the other, but these are not relevant here, as the statistic
being used is just the least-squared difference.
The two models can be visually compared by taking advantage of the previous plot objects retaining the values from the previous fit:
>>> mplot2 = ModelPlot()
>>> mplot2.prepare(d, mdl)
>>> mplot.plot()
>>> mplot2.overplot()
An alternative would be to create the plots directly (the order=2 parameter values are restored from the res object created from the first fit to the data), in which case we are not limited to calculating the model on the independent axis of the input data (the order is chosen to match the colors of the previous plot):
>>> xgrid = np.linspace(0, 20, 21)
>>> y1 = mdl(xgrid)
>>> mdl.c0 = res.parvals[0]
>>> mdl.c1 = 0
>>> mdl.c2 = res.parvals[1]
>>> y2 = mdl(xgrid)
>>> plt.clf()
>>> plt.plot(xgrid, y2, label='order=2');
>>> plt.plot(xgrid, y1, label='order=1');
>>> plt.legend();
>>> plt.title("Manual evaluation of the models");
Simple user model¶
This example works through a fit to a small one-dimensional dataset which includes errors. This means that, unlike the Simple Interpolation example, an analysis of the parameter errors can be made. The fit begins with the use of the basic Sherpa models, but this turns out to be sub-optimal - since the model parameters do not match the required parameters - so a user model is created, which recasts the Sherpa model parameters into the desired form. It also has the advantage of simplifying the model, which avoids the need for manual intervention required with the Sherpa version.
Introduction¶
For this example, a data set from the literature was chosen, looking at non-Astronomy papers to show that Sherpa can be used in a variety of fields. There is no attempt made here to interpret the results, and the model parameters, and their errors, derived here should not be assumed to have any meaning compared to the results of the paper.
The data used in this example is taken from
Zhao Y, Tang X, Zhao X, Wang Y (2017) Effect of various nitrogen
conditions on population growth, temporary cysts and cellular biochemical
compositions of Karenia mikimotoi. PLoS ONE 12(2): e0171996.
doi:10.1371/journal.pone.0171996. The
Supporting information section of the paper includes a
spreadsheet containing the data for the figures, and this was
downloaded and stored as the file pone.0171996.s001.xlsx.
The aim is to fit a similar model to that described in Table 5, that is
where \(t\) and \(y\) are the abscissa (independent axis) and ordinate (dependent axis), respectively. The idea is to see if we can get a similar result rather than to make any inferences based on the data. For this example only the “NaNO3” dataset is going to be used.
Setting up¶
Both NumPy and Matplotlib are required:
>>> import numpy as np
>>> import matplotlib.pyplot as plt
Reading in the data¶
The openpyxl package (version 2.5.3) is used to read in the data from the Excel spreadsheet. This is not guaranteed to be the optimal means of reading in the data (and relies on hard-coded knowledge of the column numbers):
>>> from openpyxl import load_workbook
>>> wb = load_workbook('pone.0171996.s001.xlsx')
>>> fig4 = wb['Fig4data']
>>> t = []; y = []; dy = []
>>> for r in list(fig4.values)[2:]:
... t.append(r[0])
... y.append(r[3])
... dy.append(r[4])
...
With these arrays, a data object can be created:
>>> from sherpa.data import Data1D
>>> d = Data1D('NaNO_3', t, y, dy)
Unlike the first worked example,
this data set includes an error column, so the data plot
created by DataPlot contains
error bars (although not obvious for the first point,
which has an error of 0):
>>> from sherpa.plot import DataPlot
>>> dplot = DataPlot()
>>> dplot.prepare(d)
>>> dplot.plot()
The data can also be inspected directly (as there aren’t many data points):
>>> print(d)
name = NaNO_3
x = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
y = [1, 4.125, 7.5, 10.0417, 12.25, 13.75, 14.0833, 14.7917, 15.5, 14.6667, 14.875]
staterror = [0, 0.9214, 1.1273, 1.9441, 2.3363, 0.9289, 1.6615, 1.1726, 1.8066, 2.149, 1.983]
syserror = None
Restricting the data¶
Trying to fit the whole data set will fail because the first data
point has an error of 0, so it is necessary to
restrict, or filter out,
this data point. The simplest way is to select a data range to ignore using
ignore(), in this
case everything where \(x < 1\):
>>> d.get_filter()
'0.0000:20.0000'
>>> d.ignore(None, 1)
>>> d.get_filter()
'2.0000:20.0000'
The get_filter() routine returns a
text description of the filters applied to the data; it starts
with all the data being included (0 to 20) and then after
excluding all points less than 1 the filter is now 2 to 20.
The format can be changed to something more appropriate for
this data set:
>>> d.get_filter(format='%d')
'2:20'
Since the data has been changed, the data plot object is updated so that the following plots reflect the new filter:
>>> dplot.prepare(d)
Creating the model¶
Table 5 lists the model fit to this dataset as
which can be constructed from components using the
Const1D
and Exp models, as shown below:
>>> from sherpa.models.basic import Const1D, Exp
>>> plateau = Const1D('plateau')
>>> rise = Exp('rise')
>>> mdl = plateau / (1 + rise)
>>> print(mdl)
(plateau / (1 + rise))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
plateau.c0 thawed 1 -3.40282e+38 3.40282e+38
rise.offset thawed 0 -3.40282e+38 3.40282e+38
rise.coeff thawed -1 -3.40282e+38 3.40282e+38
rise.ampl thawed 1 0 3.40282e+38
The amplitude of the exponential is fixed at 1, but the other
terms will remain free in the fit, with plateau.c0 representing
the normalization, and the rise.offset and rise.coeff terms
the exponent term. The offset and coeff terms do not
match the form used in the paper, namely \(a + b t\),
which has some interesting consequences for the fit, as will
be discussed below in the
user-model section.
>>> rise.ampl.freeze()
>>> print(mdl)
(plateau / (1 + rise))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
plateau.c0 thawed 1 -3.40282e+38 3.40282e+38
rise.offset thawed 0 -3.40282e+38 3.40282e+38
rise.coeff thawed -1 -3.40282e+38 3.40282e+38
rise.ampl frozen 1 0 3.40282e+38
The funtional form of the exponential model provided by Sherpa, assuming an amplitude of unity, is
which means that I expect the final values to be \({\rm coeff} \simeq -0.5\) and, as \(- {\rm coeff} * {\rm offset} \simeq 1.9\), then \({\rm offset} \simeq 4\). The plateau value should be close to 15.
The model and data can be shown together, but as the fit has not yet been made then showing on the same plot is not very instructive, so here’s two plots one above the other, created by mixing the Sherpa and Matplotlib APIs:
>>> from sherpa.plot import ModelPlot
>>> mplot = ModelPlot()
>>> mplot.prepare(d, mdl)
>>> plt.subplot(2, 1, 1)
>>> mplot.plot(clearwindow=False)
>>> plt.subplot(2, 1, 2)
>>> dplot.plot(clearwindow=False)
>>> plt.title('')
The title of the data plot was removed since it overlaped the X axis of the model plot above it.
Fitting the data¶
The main difference to fitting the first example is that the
Chi2 statistic is used,
since the data contains error values.
>>> from sherpa.stats import Chi2
>>> from sherpa.fit import Fit
>>> f = Fit(d, mdl, stat=Chi2())
>>> print(f)
data = NaNO_3
model = (plateau / (1 + rise))
stat = Chi2
method = LevMar
estmethod = Covariance
>>> print("Starting statistic: {}".format(f.calc_stat()))
Starting statistic: 633.2233812020354
The use of a Chi-square statistic means that the fit also calculates the reduced statistic (the final statistic value divided by the degrees of freedom), which should be \(\sim 1\) for a “good” fit, and an estimate of the probability (Q value) that the fit is good (this is also based on the statistic and number of degrees of freedom).
>>> fitres = f.fit()
>>> print(fitres.format())
Method = levmar
Statistic = chi2
Initial fit statistic = 633.223
Final fit statistic = 101.362 at function evaluation 17
Data points = 10
Degrees of freedom = 7
Probability [Q-value] = 5.64518e-19
Reduced statistic = 14.4802
Change in statistic = 531.862
plateau.c0 10.8792 +/- 0.428815
rise.offset 457.221 +/- 0
rise.coeff 24.3662 +/- 0
Changed in version 4.10.1: The implementation of the LevMar
class has been changed from Fortran to C++ in the 4.10.1 release.
The results of the optimiser are expected not to change
significantly, but one of the more-noticeable changes is that
the covariance matrix is now returned directly from a fit,
which results in an error estimate provided as part of the
fit output (the values after the +/- terms above).
The reduced chi-square value is large, as shown in the screen output above and the explicit access below, the probability value is essentially 0, and the parameters are nowhere near the expected values.
>>> print("Reduced chi square = {:.2f}".format(fitres.rstat))
Reduced chi square = 14.48
Visually comparing the model and data values highlights how poor
this fit is (the data plot does not need regenerating in this
case, but prepare() is called
just to make sure that the correct data is being displayed):
>>> dplot.prepare(d)
>>> mplot.prepare(d, mdl)
>>> dplot.plot()
>>> mplot.overplot()
Either the model has got caught in a local minimum, or it is not
a good description of the data. To investigate further, a useful technique
is to switch the optimiser and re-fit; the hope is that the different
optimiser will be able to escape the local minima in the search
space. The default optimiser used by
Fit is
LevMar, which is often a good
choice for data with errors. The other standard optimiser
provided by Sherpa is
NelderMead, which is often slower
than LevMar - as it requires more model evaluations - but
less-likely to get stuck:
>>> from sherpa.optmethods import NelderMead
>>> f.method = NelderMead()
>>> fitres2 = f.fit()
>>> print(mdl)
(plateau / (1 + rise))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
plateau.c0 thawed 10.8792 -3.40282e+38 3.40282e+38
rise.offset thawed 457.221 -3.40282e+38 3.40282e+38
rise.coeff thawed 24.3662 -3.40282e+38 3.40282e+38
rise.ampl frozen 1 0 3.40282e+38
An alternative to replacing the
method attribute, as done above, would be
to create a new Fit object - changing the
method using the method attribute of the initializer, and use
that to fit the model and data.
As can be seen, the parameter values have not changed; the
dstatval attribute contains the
change in the statsitic value, and as shown below, it has
not improved:
>>> fitres2.dstatval
0.0
The failure of this fit is actually down to the coupling of
the offset and coeff parameters of the
Exp model, as will be
discussed below,
but a good solution can be found by tweaking the starting
parameter values.
Restarting the fit¶
The reset() will change the
parameter values back to the
last values you set them to,
which may not be the same as their
default settings
(in this case the difference is in the state of the rise.ampl
parameter, which has remained frozen):
>>> mdl.reset()
>>> print(mdl)
(plateau / (1 + rise))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
plateau.c0 thawed 1 -3.40282e+38 3.40282e+38
rise.offset thawed 0 -3.40282e+38 3.40282e+38
rise.coeff thawed -1 -3.40282e+38 3.40282e+38
rise.ampl frozen 1 0 3.40282e+38
Note
It is not always necessary to reset the parameter values when trying to get out of a local minimum, but it can be a useful strategy to avoid getting trapped in the same area.
One of the simplest changes to make here is to set the plateau term to the maximum data value, as the intention is for this term to represent the asymptote of the curve.
>>> plateau.c0 = np.max(d.y)
>>> mplot.prepare(d, mdl)
>>> dplot.plot()
>>> mplot.overplot()
A new fit object could be created, but it is also possible
to re-use the existing object. This leaves the optimiser set to
NelderMead, although in this
case the same parameter values are found if the method
attribute had been changed back to
LevMar:
>>> fitres3 = f.fit()
>>> print(fitres3.format())
Method = neldermead
Statistic = chi2
Initial fit statistic = 168.42
Final fit statistic = 0.299738 at function evaluation 42
Data points = 10
Degrees of freedom = 7
Probability [Q-value] = 0.9999
Reduced statistic = 0.0428198
Change in statistic = 168.12
plateau.c0 14.9694 +/- 0.859633
rise.offset 4.17729 +/- 0.630148
rise.coeff -0.420696 +/- 0.118487
These results already look a lot better than the previous attempt; the reduced statistic is much smaller, and the values are similar to the reported values. As shown in the plot below, the model also well describes the data:
>>> mplot.prepare(d, mdl)
>>> dplot.plot()
>>> mplot.overplot()
The residuals can also be displayed, in this case normalizing by
the error values by using a
DelchiPlot plot:
>>> from sherpa.plot import DelchiPlot
>>> residplot = DelchiPlot()
>>> residplot.prepare(d, mdl, f.stat)
>>> residplot.plot()
Unlike the data and model plots, the
prepare() method of the
residual plot requires a statistic object, so the value
in the fit object (using the stat
attribute) is used.
Given that the reduced statistic for the fit is a lot smaller than 1 (\(\sim 0.04\)), the residuals are all close to 0: the ordinate axis shows \((d - m) / e\) where \(d\), \(m\), and \(e\) are data, model, and error value respectively.
What happens at \(t = 0\)?¶
The filtering applied earlier
can be removed, to see how the model behaves at low times. Calling
the notice() without any arguments
removes any previous filter:
>>> d.notice()
>>> d.get_filter(format='%d')
'0:20'
For this plot, the FitPlot class is going
to be used to show both the data and model rather than doing it
manually as above:
>>> from sherpa.plot import FitPlot
>>> fitplot = FitPlot()
>>> dplot.prepare(d)
>>> mplot.prepare(d, mdl)
>>> fitplot.prepare(dplot, mplot)
>>> fitplot.plot()
Note
The prepare method on the
components of the Fit plot (in this case dplot and
mplot) must be called with their appropriate arguments
to ensure that the latest changes - such as filters and
parameter values - are picked up.
Warning
Trying to create a residual plot for this new data range,
will end up with a division-by-zero warning from the
prepare call, as the first data point has an error
of 0 and the residual plot shows \((d - m) / e\).
For the rest of this example the first data point has been removed:
>>> d.ignore(None, 1)
Estimating parameter errors¶
The calc_stat_info() method returns
an overview of the current fit:
>>> statinfo = f.calc_stat_info()
>>> print(statinfo)
name =
ids = None
bkg_ids = None
statname = chi2
statval = 0.2997382864907501
numpoints = 10
dof = 7
qval = 0.999900257643
rstat = 0.04281975521296431
It is another way of getting at some of the information in the
FitResults object; for instance
>>> statinfo.rstat == fitres3.rstat
True
Note
The FitResults object refers to the model at the time
the fit was made, whereas calc_stat_info is calculated
based on the current values, and so the results can be
different.
The est_errors() method is used to
estimate error ranges for the parameter values. It does this by
varying the parameters around the best-fit location
until the statistic value has increased by a set amount.
The default method for estimating errors is
Covariance
>>> f.estmethod.name
'covariance'
which has the benefit of being fast, but may not be as robust as other techniques.
>>> coverrs = f.est_errors()
>>> print(coverrs.format())
Confidence Method = covariance
Iterative Fit Method = None
Fitting Method = levmar
Statistic = chi2
covariance 1-sigma (68.2689%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
plateau.c0 14.9694 -0.880442 0.880442
rise.offset 4.17729 -0.646012 0.646012
rise.coeff -0.420696 -0.12247 0.12247
These errors are similar to those reported during the fit.
As shown below,
the error values can be extracted from the output of
est_errors().
The default is to calculate “one sigma” error bounds
(i.e. those that cover 68.3% of the expected parameter range),
but this can be changed by altering the
sigma attribute
of the error estimator.
>>> f.estmethod.sigma
1
Changing this value to 1.6 means that the errors are close to the 90% bounds (for a single parameter):
>>> f.estmethod.sigma = 1.6
>>> coverrs90 = f.est_errors()
>>> print(coverrs90.format())
Confidence Method = covariance
Iterative Fit Method = None
Fitting Method = neldermead
Statistic = chi2
covariance 1.6-sigma (89.0401%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
plateau.c0 14.9694 -1.42193 1.42193
rise.offset 4.17729 -1.04216 1.04216
rise.coeff -0.420696 -0.19679 0.19679
The covariance method uses the covariance matrix to estimate
the error surface, and so the parameter errors are symmetric.
A more-robust, but often significantly-slower, approach is to
use the Confidence approach:
>>> from sherpa.estmethods import Confidence
>>> f.estmethod = Confidence()
>>> conferrs = f.est_errors()
plateau.c0 lower bound: -0.804259
rise.offset lower bound: -0.590258
rise.coeff lower bound: -0.148887
rise.offset upper bound: 0.714407
plateau.c0 upper bound: 0.989664
rise.coeff upper bound: 0.103391
The error estimation for the confidence technique is run in parallel - if the machine has multiple cores usable by the Python multiprocessing module - which can mean that the screen output above is not always in the same order. As shown below, the confidence-derived error bounds are similar to the covariance bounds, but are not symmetric.
>>> print(conferrs.format())
Confidence Method = confidence
Iterative Fit Method = None
Fitting Method = neldermead
Statistic = chi2
confidence 1-sigma (68.2689%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
plateau.c0 14.9694 -0.804259 0.989664
rise.offset 4.17729 -0.590258 0.714407
rise.coeff -0.420696 -0.148887 0.103391
The default is to use all
thawed parameters
in the error analysis, but the est_errors()
method has a parlist attribute which can be used to restrict
the parameters used, for example to just the offset term:
>>> offseterrs = f.est_errors(parlist=(mdl.pars[1], ))
rise.offset lower bound: -0.590258
rise.offset upper bound: 0.714407
>>> print(offseterrs)
datasets = None
methodname = confidence
iterfitname = none
fitname = neldermead
statname = chi2
sigma = 1
percent = 68.2689492137
parnames = ('rise.offset',)
parvals = (4.177287700807689,)
parmins = (-0.59025803525842369,)
parmaxes = (0.71440700826435144,)
nfits = 8
The covariance and confidence limits can be compared by
accessing the fields of the
ErrorEstResults object:
>>> fmt = "{:13s} covar=±{:4.2f} conf={:+5.2f} {:+5.2f}"
>>> for i in range(len(conferrs.parnames)):
... print(fmt.format(conferrs.parnames[i], coverrs.parmaxes[i],
... conferrs.parmins[i], conferrs.parmaxes[i]))
...
plateau.c0 covar=±0.88 conf=-0.80 +0.99
rise.offset covar=±0.65 conf=-0.59 +0.71
rise.coeff covar=±0.12 conf=-0.15 +0.10
The est_errors() method returns a
range, but often it is important to visualize the error
surface, which can be done using the interval projection
(for one parameter) and region projection (for two parameter)
routines. The one-dimensional version is created with the
IntervalProjection
class, as shown in the following, which shows how the statistic
varies with the plateau term (the vertical dashed line indicates
the best-fit location for the parameter, and the horizontal
line the statistic value for the best-fit location):
>>> from sherpa.plot import IntervalProjection
>>> intproj = IntervalProjection()
>>> intproj.calc(f, plateau.c0)
>>> intproj.plot()
Unlike the previous plots, this requires calling the
calc() method
before plot(). As
the prepare()
method was not called, it used the default options to
calculate the plot range (i.e. the range over which
plateau.c0 would be varied), which turns out in this
case to be close to the one-sigma limits.
The range, and number of points, can also be set explicitly:
>>> intproj.prepare(min=12.5, max=20, nloop=51)
>>> intproj.calc(f, plateau.c0)
>>> intproj.plot()
>>> s0 = f.calc_stat()
>>> for ds in [1, 4, 9]:
... intproj.hline(s0 + ds, overplot=True, linestyle='dot', linecolor='gray')
...
The horizontal lines indicate the statistic value for one, two, and three sigma limits for a single parameter value (and assuming a Chi-square statistic). The plot shows how, as the parameter moves away from its best-fit location, the search space becomes less symmetric.
Following the same approach, the RegionProjection
class calculates the statistic value as two parameters are varied,
displaying the results as a contour plot. It requires two parameters
and the visualization is
created with the contour()
method:
>>> from sherpa.plot import RegionProjection
>>> regproj = RegionProjection()
>>> regproj.calc(f, rise.offset, rise.coeff)
>>> regproj.contour()
The contours show the one, two, and three sigma contours, with the
cross indicating the best-fit value. As with the interval-projection plot,
the prepare() method can
be used to define the grid of points to use; the values below are
chosen to try and cover the full three-sigma range as well as improve
the smoothness of the contours by increasing the number of points
that are looped over:
>>> regproj.prepare(min=(2, -1.2), max=(8, -0.1), nloop=(21, 21))
>>> regproj.calc(f, rise.offset, rise.coeff)
>>> regproj.contour()
Writing your own model¶
The approach above has provided fit results, but they do not match those of the paper and, since
it is hard to transform the values from above to get
accurate results. An alternative approach is to
create a model with the parameters
in the required form, which requires a small amount
of code (by using the
Exp class to do the actual
model evaluation).
The following class (MyExp) creates a model that has
two parameters (a and b) that represents
\(f(x) = e^{a + b x}\). The starting values for these
parameters are chosen to match the default values of the
Exp parameters,
where \({\rm coeff} = -1\) and \({\rm offset} = 0\):
from sherpa.models.basic import RegriddableModel1D
from sherpa.models.parameter import Parameter
class MyExp(RegriddableModel1D):
"""A simpler form of the Exp model.
The model is f(x) = exp(a + b * x).
"""
def __init__(self, name='myexp'):
self.a = Parameter(name, 'a', 0)
self.b = Parameter(name, 'b', -1)
# The _exp instance is used to perform the model calculation,
# as shown in the calc method.
self._exp = Exp('hidden')
return RegriddableModel1D.__init__(self, name, (self.a, self.b))
def calc(self, pars, *args, **kwargs):
"""Calculate the model"""
# Tell the exp model to evaluate the model, after converting
# the parameter values to the required form, and order, of:
# offset, coeff, ampl.
#
coeff = pars[1]
offset = -1 * pars[0] / coeff
ampl = 1.0
return self._exp.calc([offset, coeff, ampl], *args, **kwargs)
This can be used as any other Sherpa model:
>>> plateau2 = Const1D('plateau2')
>>> rise2 = MyExp('rise2')
>>> mdl2 = plateau2 / (1 + rise2)
>>> print(mdl2)
(plateau2 / (1 + rise2))
Param Type Value Min Max Units
----- ---- ----- --- --- -----
plateau2.c0 thawed 1 -3.40282e+38 3.40282e+38
rise2.a thawed 0 -3.40282e+38 3.40282e+38
rise2.b thawed -1 -3.40282e+38 3.40282e+38
>>> fit2 = Fit(d, mdl2, stat=Chi2())
>>> res2 = fit2.fit()
>>> print(res2.format())
Method = levmar
Statistic = chi2
Initial fit statistic = 633.223
Final fit statistic = 0.299738 at function evaluation 52
Data points = 10
Degrees of freedom = 7
Probability [Q-value] = 0.9999
Reduced statistic = 0.0428198
Change in statistic = 632.924
plateau2.c0 14.9694 +/- 0.859768
rise2.a 1.75734 +/- 0.419169
rise2.b -0.420685 +/- 0.118473
>>> dplot.prepare(d)
>>> mplot2 = ModelPlot()
>>> mplot2.prepare(d, mdl2)
>>> dplot.plot()
>>> mplot2.overplot()
Unlike the initial attempt,
this version did not require any manual intervention to find the
best-fit solution. This is because the degeneracy between the two
terms of the exponential in the
Exp model have been broken in
this version, and so the optimiser work better.
It also has the advantage that the parameters match the
problem, and so the parameter limits determined below can be
used directly, without having to transform them.
>>> fit2.estmethod = Confidence()
>>> conferrs2 = fit2.est_errors()
plateau2.c0 lower bound: -0.804444
rise2.b lower bound: -0.148899
rise2.a lower bound: -0.38086
rise2.b upper bound: 0.10338
plateau2.c0 upper bound: 0.989623
rise2.a upper bound: 0.489919
>>> print(conferrs2.format())
Confidence Method = confidence
Iterative Fit Method = None
Fitting Method = levmar
Statistic = chi2
confidence 1-sigma (68.2689%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
plateau2.c0 14.9694 -0.804444 0.989623
rise2.a 1.75734 -0.38086 0.489919
rise2.b -0.420685 -0.148899 0.10338
The difference in the model parameterisation can also be seen in the various error-analysis plots, such as the region-projection contour plot (where the limits have been chosen to cover the three-sigma contour), and a marker has been added to show the result listed in Table 5 of Zhao et al:
>>> regproj2 = RegionProjection()
>>> regproj2.prepare(min=(0.5, -1.2), max=(5, -0.1), nloop=(21, 21))
>>> regproj2.calc(fit2, rise2.a, rise2.b)
>>> regproj2.contour()
>>> plt.plot(1.941, -0.453, 'ko', label='NaNO$_3$ Table 5')
>>> plt.legend(loc=1)
Using Sessions to manage models and data¶
So far we have discussed the object-based API of Sherpa -
where it is up to the user to manage the creation
and handling of
data,
model,
fit and related objects. Sherpa
also provides a “Session” class that handles much of this,
and it can be used
directly - via the sherpa.ui.utils.Session or
sherpa.astro.ui.utils.Session classes - or
indirectly using the routines in the
sherpa.ui and sherpa.astro.ui modules.
The session API is intended to be used in an interactive setting, and so deals with object management. Rather than deal with objects, the API uses labels (numeric or string) to identify data sets and model components. The Astronomy-specific version adds domain-specific functionality; in this case support for Astronomical data analysis, with a strong focus on high-energy (X-ray) data. It is is currently documented on the http://cxc.harvard.edu/sherpa/ web site.
The Session object provides methods
that allow you to:
- load data
- set the model
- change the statistic and optimiser
- fit
- calculate errors
- visualize the results
These are the same stages as described in the getting started section, but the syntax is different, since the Session object handles the creation of, and passing around, the underlying Sherpa objects.
The sherpa.ui module provides an interface where
the Session object is hidden from the user, which makes it
more appropriate for an interactive analysis session.
Examples¶
The following examples are very basic, since they are intended to highlight how the Sesssion API is used. The CIAO documentation for Sherpa at http://cxc.harvard.edu/sherpa/ provides more documentation and examples.
There are two examples which show the same process -
finding out what value best represents a small dataset -
using the
Session object directly and then via the
sherpa.ui module.
The data to be fit is the four element array:
>>> x = [100, 200, 300, 400]
>>> y = [10, 12, 9, 13]
For this example the Cash statistic will
be used, along with the
NelderMead optimiser.
Note
Importing the Session object - whether directly or via the ui module - causes several checks to be run, to see what parts of the system may not be available. This can lead to warning messages such as the following to be displayed:
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'
Other checks are to see if the chosen I/O and plotting backends are present, and if support for the XSPEC model library is available.
Using the Session object¶
By default the Session object has no available models associated
with it. The
_add_model_types()
method is used to register the models from
sherpa.models.basic with the session (by default it will
add any class in the module that is derived from the
ArithmeticModel
class):
>>> from sherpa.ui.utils import Session
>>> import sherpa.models.basic
>>> s = Session()
>>> s._add_model_types(sherpa.models.basic)
The load_arrays() is used to
create a Data1D object, which is managed
by the Session class and referenced by the identifier 1
(this is in fact the default identifier, which can be
manipulated by the
get_default_id()
and
set_default_id()
methods, and can be a string or an integer).
Many methods will default to using the default identifier,
but load_arrays requires it:
>>> s.load_arrays(1, x, y)
Note
The session object is not just limited to handling
Data1D data sets. The
load_arrays takes an optional argument which defines
the class of the data (e.g. Data2D),
and there are several other methods which can be used to
create a data object, such as
load_data
and
set_data.
The list_data_ids() method
returns the list of available data sets (i.e. those that have
been loaded into the session):
>>> s.list_data_ids()
[1]
The get_data() method lets a user
access the underlying data object. This method uses the default
identifier if not specified:
>>> s.get_data()
<Data1D data set instance ''>
>>> print(s.get_data())
name =
x = Int64[4]
y = Int64[4]
staterror = None
syserror = None
The default statistic and optimiser are set to values useful for data with Gaussian errors:
>>> s.get_stat_name()
'chi2gehrels'
>>> s.get_method_name()
'levmar'
As the data here is counts based, and is to be fit with Poisson
statitics, the
set_stat()
and
set_method()
methods are used to change the statistic and optimiser.
Note that they take a string as an argument
(rather than an instance of a
Stat
or OptMethod
class):
>>> s.set_stat('cash')
>>> s.set_method('simplex')
The set_source() method is
used to define the model expression that is to be fit to the
data. It can be sent a model expression created using the
model classes directly, as described in the
Creating Model Instances section above.
However, in this case a string is used to define the model, and
references each model component using the form
modelname.instancename. The modelname defines the
type of model - in this case the
Const1D model - and it must
have been registered with the session object using
_add_model_types. The
list_models() method
can be used to find out what models are available.
The instancename is used as an identifier for the
component, and can be used with other methods,
such as set_par().
>>> s.set_source('const1d.mdl')
The instancename value is also used to create a Python variable
which provides direct access to the model component (it can
also be retrieved with
get_model_component()):
>>> print(mdl)
const1d.mdl
Param Type Value Min Max Units
----- ---- ----- --- --- -----
mdl.c0 thawed 1 -3.40282e+38 3.40282e+38
The source model can be retrievd with
get_source(), which in this
example is just the single model component mdl:
>>> s.get_source()
<Const1D model instance 'const1d.mdl'>
With the data, model, statistic, and optimiser set, it
is now possible to perform a fit. The
fit() method defaults to
a simultaneous fit of all the loaded data sets; in this case
there is only one:
>>> s.fit()
Dataset = 1
Method = neldermead
Statistic = cash
Initial fit statistic = 8
Final fit statistic = -123.015 at function evaluation 90
Data points = 4
Degrees of freedom = 3
Change in statistic = 131.015
mdl.c0 11
The fit results are displayed to the screen, but can also be accessed
with methods such as
calc_stat(),
calc_stat_info(),
and
get_fit_results().
>>> r = s.get_fit_results()
>>> print(r)
datasets = (1,)
itermethodname = none
methodname = neldermead
statname = cash
succeeded = True
parnames = ('mdl.c0',)
parvals = (11.0,)
statval = -123.01478400625663
istatval = 8.0
dstatval = 131.014784006
numpoints = 4
dof = 3
qval = None
rstat = None
message = Optimization terminated successfully
nfev = 90
There are also methods which allow you to plot the data, model,
fit, and residuals (amongst others):
plot_data(),
plot_model(),
plot_fit(),
plot_resid().
The following
hides the automatically-created error bars on the data points
(but unfortunately not the warning message)
by changing a setting in dictionary returned by
get_data_plot_prefs(),
and then displays the data along with the model:
>>> s.get_data_plot_prefs()['yerrorbars'] = False
>>> s.plot_fit()
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with leastsq
Using the UI module¶
Using the UI module is very similar to the Session object, since it automatically creates a global Session object, and registers the available models, when imported. This means that the preceeding example can be replicated but without the need for the Session object.
Since the module is intended for an interactive environment, in this example the symbols are loaded into the default namespace to avoid having to qualify each function with the module name. For commentary, please refer to the preceeding example:
>>> from sherpa.ui import *
>>> load_arrays(1, x, y)
>>> list_data_ids()
[1]
>>> get_data()
<Data1D data set instance ''>
>>> print(get_data())
name =
x = Int64[4]
y = Int64[4]
staterror = None
syserror = None
>>> get_stat_name()
'chi2gehrels'
>>> get_method_name()
'levmar'
>>> set_stat('cash')
>>> set_method('simplex')
>>> set_source('const1d.mdl')
>>> print(mdl)
const1d.mdl
Param Type Value Min Max Units
----- ---- ----- --- --- -----
mdl.c0 thawed 1 -3.40282e+38 3.40282e+38
>>> get_source()
<Const1D model instance 'const1d.mdl'>
>>> fit()
Dataset = 1
Method = neldermead
Statistic = cash
Initial fit statistic = 8
Final fit statistic = -123.015 at function evaluation 90
Data points = 4
Degrees of freedom = 3
Change in statistic = 131.015
mdl.c0 11
>>> r = get_fit_results()
>>> print(r)
datasets = (1,)
itermethodname = none
methodname = neldermead
statname = cash
succeeded = True
parnames = ('mdl.c0',)
parvals = (11.0,)
statval = -123.01478400625663
istatval = 8.0
dstatval = 131.014784006
numpoints = 4
dof = 3
qval = None
rstat = None
message = Optimization terminated successfully
nfev = 90
>>> get_data_plot_prefs()['yerrorbars'] = False
>>> plot_fit()
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with leastsq
The plot created by this function is the same as shown in the previous example.
Reference/API¶
The sherpa.ui module¶
The sherpa.ui module provides an interface to the
sherpa.ui.utils.Session
object, where a singleton class is used to provide the access but
hidden away. This needs better explanation…
Functions
add_model(modelclass[, args, kwargs])Create a user-defined model class. add_user_pars(modelname, parnames[, …])Add parameter information to a user model. calc_chisqr([id])Calculate the per-bin chi-squared statistic. calc_stat([id])Calculate the fit statistic for a data set. calc_stat_info()Display the statistic values for the current models. clean()Clear out the current Sherpa session. conf(*args)Estimate parameter confidence intervals using the confidence method. confidence(*args)Estimate parameter confidence intervals using the confidence method. contour(*args)Create a contour plot for an image data set. contour_data([id])Contour the values of an image data set. contour_fit([id])Contour the fit to a data set. contour_fit_resid([id, replot, overcontour])Contour the fit and the residuals to a data set. contour_kernel([id])Contour the kernel applied to the model of an image data set. contour_model([id])Create a contour plot of the model. contour_psf([id])Contour the PSF applied to the model of an image data set. contour_ratio([id])Contour the ratio of data to model. contour_resid([id])Contour the residuals of the fit. contour_source([id])Create a contour plot of the unconvolved spatial model. copy_data(fromid, toid)Copy a data set, creating a new identifier. covar(*args)Estimate parameter confidence intervals using the covariance method. covariance(*args)Estimate parameter confidence intervals using the covariance method. create_model_component([typename, name])Create a model component. dataspace1d(start, stop[, step, numbins, …])Create the independent axis for a 1D data set. dataspace2d(dims[, id, dstype])Create the independent axis for a 2D data set. delete_data([id])Delete a data set by identifier. delete_model([id])Delete the model expression for a data set. delete_model_component(name)Delete a model component. delete_psf([id])Delete the PSF model for a data set. fake([id, method])Simulate a data set. fit([id])Fit a model to one or more data sets. freeze(*args)Fix model parameters so they are not changed by a fit. get_cdf_plot()Return the data used to plot the last CDF. get_chisqr_plot([id])Return the data used by plot_chisqr. get_conf()Return the confidence-interval estimation object. get_conf_opt([name])Return one or all of the options for the confidence interval method. get_conf_results()Return the results of the last conf run. get_confidence_results()Return the results of the last conf run. get_covar()Return the covariance estimation object. get_covar_opt([name])Return one or all of the options for the covariance method. get_covar_results()Return the results of the last covar run. get_covariance_results()Return the results of the last covar run. get_data([id])Return the data set by identifier. get_data_contour([id])Return the data used by contour_data. get_data_contour_prefs()Return the preferences for contour_data. get_data_image([id])Return the data used by image_data. get_data_plot([id])Return the data used by plot_data. get_data_plot_prefs()Return the preferences for plot_data. get_default_id()Return the default data set identifier. get_delchi_plot([id])Return the data used by plot_delchi. get_dep([id, filter])Return the dependent axis of a data set. get_dims([id, filter])Return the dimensions of the data set. get_draws([id, otherids, niter, covar_matrix])Run the pyBLoCXS MCMC algorithm. get_error([id, filter])Return the errors on the dependent axis of a data set. get_filter([id])Return the filter expression for a data set. get_fit_contour([id])Return the data used by contour_fit. get_fit_plot([id])Return the data used to create the fit plot. get_fit_results()Return the results of the last fit. get_functions()Return the functions provided by Sherpa. get_indep([id])Return the independent axes of a data set. get_int_proj([par, id, otherids, recalc, …])Return the interval-projection object. get_int_unc([par, id, otherids, recalc, …])Return the interval-uncertainty object. get_iter_method_name()Return the name of the iterative fitting scheme. get_iter_method_opt([optname])Return one or all options for the iterative-fitting scheme. get_kernel_contour([id])Return the data used by contour_kernel. get_kernel_image([id])Return the data used by image_kernel. get_kernel_plot([id])Return the data used by plot_kernel. get_method([name])Return an optimization method. get_method_name()Return the name of current Sherpa optimization method. get_method_opt([optname])Return one or all of the options for the current optimization method. get_model([id])Return the model expression for a data set. get_model_autoassign_func()Return the method used to create model component identifiers. get_model_component(name)Returns a model component given its name. get_model_component_image(id[, model])Return the data used by image_model_component. get_model_component_plot(id[, model])Return the data used to create the model-component plot. get_model_contour([id])Return the data used by contour_model. get_model_contour_prefs()Return the preferences for contour_model. get_model_image([id])Return the data used by image_model. get_model_pars(model)Return the names of the parameters of a model. get_model_plot([id])Return the data used to create the model plot. get_model_plot_prefs()Return the preferences for plot_model. get_model_type(model)Describe a model expression. get_num_par([id])Return the number of parameters in a model expression. get_num_par_frozen([id])Return the number of frozen parameters in a model expression. get_num_par_thawed([id])Return the number of thawed parameters in a model expression. get_par(par)Return a parameter of a model component. get_pdf_plot()Return the data used to plot the last PDF. get_prior(par)Return the prior function for a parameter (MCMC). get_proj()Return the confidence-interval estimation object. get_proj_opt([name])Return one or all of the options for the confidence interval method. get_proj_results()Return the results of the last proj run. get_projection_results()Return the results of the last proj run. get_psf([id])Return the PSF model defined for a data set. get_psf_contour([id])Return the data used by contour_psf. get_psf_image([id])Return the data used by image_psf. get_psf_plot([id])Return the data used by plot_psf. get_pvalue_plot([null_model, alt_model, …])Return the data used by plot_pvalue. get_pvalue_results()Return the data calculated by the last plot_pvalue call. get_ratio_contour([id])Return the data used by contour_ratio. get_ratio_image([id])Return the data used by image_ratio. get_ratio_plot([id])Return the data used by plot_ratio. get_reg_proj([par0, par1, id, otherids, …])Return the region-projection object. get_reg_unc([par0, par1, id, otherids, …])Return the region-uncertainty object. get_resid_contour([id])Return the data used by contour_resid. get_resid_image([id])Return the data used by image_resid. get_resid_plot([id])Return the data used by plot_resid. get_sampler()Return the current MCMC sampler options. get_sampler_name()Return the name of the current MCMC sampler. get_sampler_opt(opt)Return an option of the current MCMC sampler. get_scatter_plot()Return the data used to plot the last scatter plot. get_source([id])Return the source model expression for a data set. get_source_component_image(id[, model])Return the data used by image_source_component. get_source_component_plot(id[, model])Return the data used by plot_source_component. get_source_contour([id])Return the data used by contour_source. get_source_image([id])Return the data used by image_source. get_source_plot([id])Return the data used to create the source plot. get_split_plot()Return the plot attributes for displays with multiple plots. get_stat([name])Return the fit statisic. get_stat_info()Return the statistic values for the current models. get_stat_name()Return the name of the current fit statistic. get_staterror([id, filter])Return the statistical error on the dependent axis of a data set. get_syserror([id, filter])Return the systematic error on the dependent axis of a data set. get_trace_plot()Return the data used to plot the last trace. guess([id, model, limits, values])Estimate the parameter values and ranges given the loaded data. ignore([lo, hi])Exclude data from the fit. ignore_id(ids[, lo, hi])Exclude data from the fit for a data set. image_close()Close the image viewer. image_data([id, newframe, tile])Display a data set in the image viewer. image_deleteframes()Delete all the frames open in the image viewer. image_fit([id, newframe, tile, deleteframes])Display the data, model, and residuals for a data set in the image viewer. image_getregion([coord])Return the region defined in the image viewer. image_kernel([id, newframe, tile])Display the 2D kernel for a data set in the image viewer. image_model([id, newframe, tile])Display the model for a data set in the image viewer. image_model_component(id[, model, newframe, …])Display a component of the model in the image viewer. image_open()Start the image viewer. image_psf([id, newframe, tile])Display the 2D PSF model for a data set in the image viewer. image_ratio([id, newframe, tile])Display the ratio (data/model) for a data set in the image viewer. image_resid([id, newframe, tile])Display the residuals (data - model) for a data set in the image viewer. image_setregion(reg[, coord])Set the region to display in the image viewer. image_source([id, newframe, tile])Display the source expression for a data set in the image viewer. image_source_component(id[, model, …])Display a component of the source expression in the image viewer. image_xpaget(arg)Return the result of an XPA call to the image viewer. image_xpaset(arg[, data])Return the result of an XPA call to the image viewer. int_proj(par[, id, otherids, replot, fast, …])Calculate and plot the fit statistic versus fit parameter value. int_unc(par[, id, otherids, replot, min, …])Calculate and plot the fit statistic versus fit parameter value. link(par, val)Link a parameter to a value. list_data_ids()List the identifiers for the loaded data sets. list_functions([outfile, clobber])Display the functions provided by Sherpa. list_iter_methods()List the iterative fitting schemes. list_methods()List the optimization methods. list_model_components()List the names of all the model components. list_model_ids()List of all the data sets with a source expression. list_models([show])List the available model types. list_priors()Return the priors set for model parameters, if any. list_samplers()List the MCMC samplers. list_stats()List the fit statistics. load_arrays(id, *args)Create a data set from array values. load_conv(modelname, filename_or_model, …)Load a 1D convolution model. load_data(id[, filename, ncols, colkeys, …])Load a data set from an ASCII file. load_filter(id[, filename, ignore, ncols])Load the filter array from an ASCII file and add to a data set. load_psf(modelname, filename_or_model, …)Create a PSF model. load_staterror(id[, filename, ncols])Load the statistical errors from an ASCII file. load_syserror(id[, filename, ncols])Load the systematic errors from an ASCII file. load_table_model(modelname, filename[, …])Load ASCII tabular data and use it as a model component. load_template_interpolator(name, …)Set the template interpolation scheme. load_template_model(modelname, templatefile)Load a set of templates and use it as a model component. load_user_model(func, modelname[, filename, …])Create a user-defined model. load_user_stat(statname, calc_stat_func[, …])Create a user-defined statistic. normal_sample([num, sigma, correlate, id, …])Sample the fit statistic by taking the parameter values from a normal distribution. notice([lo, hi])Include data in the fit. notice_id(ids[, lo, hi])Include data from the fit for a data set. paramprompt([val])Should the user be asked for the parameter values when creating a model? plot(*args)Create one or more plot types. plot_cdf(points[, name, xlabel, replot, …])Plot the cumulative density function of an array of values. plot_chisqr([id])Plot the chi-squared value for each point in a data set. plot_data([id])Plot the data values. plot_delchi([id])Plot the ratio of residuals to error for a data set. plot_fit([id])Plot the fit results (data, model) for a data set. plot_fit_delchi([id, replot, overplot, …])Plot the fit results, and the residuals, for a data set. plot_fit_resid([id, replot, overplot, …])Plot the fit results, and the residuals, for a data set. plot_kernel([id])Plot the 1D kernel applied to a data set. plot_model([id])Plot the model for a data set. plot_model_component(id[, model])Plot a component of the model for a data set. plot_pdf(points[, name, xlabel, bins, …])Plot the probability density function of an array of values. plot_psf([id])Plot the 1D PSF model applied to a data set. plot_pvalue(null_model, alt_model[, …])Compute and plot a histogram of likelihood ratios by simulating data. plot_ratio([id])Plot the ratio of data to model for a data set. plot_resid([id])Plot the residuals (data - model) for a data set. plot_scatter(x, y[, name, xlabel, ylabel, …])Create a scatter plot. plot_source([id])Plot the source expression for a data set. plot_source_component(id[, model])Plot a component of the source expression for a data set. plot_trace(points[, name, xlabel, replot, …])Create a trace plot of row number versus value. proj(*args)Estimate parameter confidence intervals using the projection method. projection(*args)Estimate parameter confidence intervals using the projection method. reg_proj(par0, par1[, id, otherids, replot, …])Plot the statistic value as two parameters are varied. reg_unc(par0, par1[, id, otherids, replot, …])Plot the statistic value as two parameters are varied. reset([model, id])Reset the model parameters to their default settings. restore([filename])Load in a Sherpa session from a file. save([filename, clobber])Save the current Sherpa session to a file. save_arrays(filename, args[, fields, …])Write a list of arrays to an ASCII file. save_data(id[, filename, fields, sep, …])Save the data to a file. save_delchi(id[, filename, clobber, sep, …])Save the ratio of residuals (data-model) to error to a file. save_error(id[, filename, clobber, sep, …])Save the errors to a file. save_filter(id[, filename, clobber, sep, …])Save the filter array to a file. save_model(id[, filename, clobber, sep, …])Save the model values to a file. save_resid(id[, filename, clobber, sep, …])Save the residuals (data-model) to a file. save_source(id[, filename, clobber, sep, …])Save the model values to a file. save_staterror(id[, filename, clobber, sep, …])Save the statistical errors to a file. save_syserror(id[, filename, clobber, sep, …])Save the statistical errors to a file. set_conf_opt(name, val)Set an option for the confidence interval method. set_covar_opt(name, val)Set an option for the covariance method. set_data(id[, data])Set a data set. set_default_id(id)Set the default data set identifier. set_dep(id[, val])Set the dependent axis of a data set. set_filter(id[, val, ignore])Set the filter array of a data set. set_full_model(id[, model])Define the convolved model expression for a data set. set_iter_method(meth)Set the iterative-fitting scheme used in the fit. set_iter_method_opt(optname, val)Set an option for the iterative-fitting scheme. set_method(meth)Set the optimization method. set_method_opt(optname, val)Set an option for the current optimization method. set_model(id[, model])Set the source model expression for a data set. set_model_autoassign_func([func])Set the method used to create model component identifiers. set_par(par[, val, min, max, frozen])Set the value, limits, or behavior of a model parameter. set_prior(par, prior)Set the prior function to use with a parameter. set_proj_opt(name, val)Set an option for the projection method. set_psf(id[, psf])Add a PSF model to a data set. set_sampler(sampler)Set the MCMC sampler. set_sampler_opt(opt, value)Set an option for the current MCMC sampler. set_source(id[, model])Set the source model expression for a data set. set_stat(stat)Set the statistical method. set_staterror(id[, val, fractional])Set the statistical errors on the dependent axis of a data set. set_syserror(id[, val, fractional])Set the systematic errors on the dependent axis of a data set. set_xlinear([plottype])New plots will display a linear X axis. set_xlog([plottype])New plots will display a logarithmically-scaled X axis. set_ylinear([plottype])New plots will display a linear Y axis. set_ylog([plottype])New plots will display a logarithmically-scaled Y axis. show_all([id, outfile, clobber])Report the current state of the Sherpa session. show_conf([outfile, clobber])Display the results of the last conf evaluation. show_covar([outfile, clobber])Display the results of the last covar evaluation. show_data([id, outfile, clobber])Summarize the available data sets. show_filter([id, outfile, clobber])Show any filters applied to a data set. show_fit([outfile, clobber])Summarize the fit results. show_kernel([id, outfile, clobber])Display any kernel applied to a data set. show_method([outfile, clobber])Display the current optimization method and options. show_model([id, outfile, clobber])Display the model expression used to fit a data set. show_proj([outfile, clobber])Display the results of the last proj evaluation. show_psf([id, outfile, clobber])Display any PSF model applied to a data set. show_source([id, outfile, clobber])Display the source model expression for a data set. show_stat([outfile, clobber])Display the current fit statistic. simulfit([id])Fit a model to one or more data sets. t_sample([num, dof, id, otherids, numcores])Sample the fit statistic by taking the parameter values from a Student’s t-distribution. thaw(*args)Allow model parameters to be varied during a fit. uniform_sample([num, factor, id, otherids, …])Sample the fit statistic by taking the parameter values from an uniform distribution. unlink(par)Unlink a parameter value. unpack_arrays(*args)Create a sherpa data object from arrays of data. unpack_data(filename[, ncols, colkeys, …])Create a sherpa data object from an ASCII file.
The sherpa.astro.ui module¶
Since some of these are re-exports of the
sherpa.ui version, should
we only document the new/different ones here? That would be a maintenance
burden…
Functions
add_model(modelclass[, args, kwargs])Create a user-defined model class. add_user_pars(modelname, parnames[, …])Add parameter information to a user model. calc_chisqr([id])Calculate the per-bin chi-squared statistic. calc_data_sum([lo, hi, id, bkg_id])Sum up the data values over a pass band. calc_data_sum2d([reg, id])Sum up the data values of a 2D data set. calc_energy_flux([lo, hi, id, bkg_id])Integrate the unconvolved source model over a pass band. calc_kcorr(z, obslo, obshi[, restlo, …])Calculate the K correction for a model. calc_model_sum([lo, hi, id, bkg_id])Sum up the fitted model over a pass band. calc_model_sum2d([reg, id])Sum up the convolved model for a 2D data set. calc_photon_flux([lo, hi, id, bkg_id])Integrate the unconvolved source model over a pass band. calc_source_sum([lo, hi, id, bkg_id])Sum up the source model over a pass band. calc_source_sum2d([reg, id])Sum up the unconvolved model for a 2D data set. calc_stat([id])Calculate the fit statistic for a data set. calc_stat_info()Display the statistic values for the current models. clean()Clear out the current Sherpa session. conf(*args)Estimate parameter confidence intervals using the confidence method. confidence(*args)Estimate parameter confidence intervals using the confidence method. contour(*args)Create a contour plot for an image data set. contour_data([id])Contour the values of an image data set. contour_fit([id])Contour the fit to a data set. contour_fit_resid([id, replot, overcontour])Contour the fit and the residuals to a data set. contour_kernel([id])Contour the kernel applied to the model of an image data set. contour_model([id])Create a contour plot of the model. contour_psf([id])Contour the PSF applied to the model of an image data set. contour_ratio([id])Contour the ratio of data to model. contour_resid([id])Contour the residuals of the fit. contour_source([id])Create a contour plot of the unconvolved spatial model. copy_data(fromid, toid)Copy a data set, creating a new identifier. covar(*args)Estimate parameter confidence intervals using the covariance method. covariance(*args)Estimate parameter confidence intervals using the covariance method. create_arf(elo, ehi[, specresp, exposure, …])Create an ARF. create_model_component([typename, name])Create a model component. create_rmf(rmflo, rmfhi[, startchan, e_min, …])Create an RMF. dataspace1d(start, stop[, step, numbins, …])Create the independent axis for a 1D data set. dataspace2d(dims[, id, dstype])Create the independent axis for a 2D data set. delete_bkg_model([id, bkg_id])Delete the background model expression for a data set. delete_data([id])Delete a data set by identifier. delete_model([id])Delete the model expression for a data set. delete_model_component(name)Delete a model component. delete_psf([id])Delete the PSF model for a data set. eqwidth(src, combo[, id, lo, hi, bkg_id, …])Calculate the equivalent width of an emission or absorption line. fake([id, method])Simulate a data set. fake_pha(id, arf, rmf, exposure[, backscal, …])Simulate a PHA data set from a model. fit([id])Fit a model to one or more data sets. fit_bkg([id])Fit a model to one or more background PHA data sets. freeze(*args)Fix model parameters so they are not changed by a fit. get_analysis([id])Return the units used when fitting spectral data. get_areascal([id, bkg_id])Return the fractional area factor of a PHA data set. get_arf([id, resp_id, bkg_id])Return the ARF associated with a PHA data set. get_arf_plot([id, resp_id])Return the data used by plot_arf. get_axes([id, bkg_id])Return information about the independent axes of a data set. get_backscal([id, bkg_id])Return the area scaling of a PHA data set. get_bkg([id, bkg_id])Return the background for a PHA data set. get_bkg_arf([id])Return the background ARF associated with a PHA data set. get_bkg_chisqr_plot([id, bkg_id])Return the data used by plot_bkg_chisqr. get_bkg_delchi_plot([id, bkg_id])Return the data used by plot_bkg_delchi. get_bkg_fit_plot([id, bkg_id])Return the data used by plot_bkg_fit. get_bkg_model([id, bkg_id])Return the model expression for the background of a PHA data set. get_bkg_model_plot([id, bkg_id])Return the data used by plot_bkg_model. get_bkg_plot([id, bkg_id])Return the data used by plot_bkg. get_bkg_ratio_plot([id, bkg_id])Return the data used by plot_bkg_ratio. get_bkg_resid_plot([id, bkg_id])Return the data used by plot_bkg_resid. get_bkg_rmf([id])Return the background RMF associated with a PHA data set. get_bkg_scale([id])Return the background scaling factor for a PHA data set. get_bkg_source([id, bkg_id])Return the model expression for the background of a PHA data set. get_bkg_source_plot([id, lo, hi, bkg_id])Return the data used by plot_bkg_source. get_cdf_plot()Return the data used to plot the last CDF. get_chisqr_plot([id])Return the data used by plot_chisqr. get_conf()Return the confidence-interval estimation object. get_conf_opt([name])Return one or all of the options for the confidence interval method. get_conf_results()Return the results of the last conf run. get_confidence_results()Return the results of the last conf run. get_coord([id])Get the coordinate system used for image analysis. get_counts([id, filter, bkg_id])Return the dependent axis of a data set. get_covar()Return the covariance estimation object. get_covar_opt([name])Return one or all of the options for the covariance method. get_covar_results()Return the results of the last covar run. get_covariance_results()Return the results of the last covar run. get_data([id])Return the data set by identifier. get_data_contour([id])Return the data used by contour_data. get_data_contour_prefs()Return the preferences for contour_data. get_data_image([id])Return the data used by image_data. get_data_plot([id])Return the data used by plot_data. get_data_plot_prefs()Return the preferences for plot_data. get_default_id()Return the default data set identifier. get_delchi_plot([id])Return the data used by plot_delchi. get_dep([id, filter, bkg_id])Return the dependent axis of a data set. get_dims([id, filter])Return the dimensions of the data set. get_draws([id, otherids, niter, covar_matrix])Run the pyBLoCXS MCMC algorithm. get_energy_flux_hist([lo, hi, id, num, …])Return the data displayed by plot_energy_flux. get_error([id, filter, bkg_id])Return the errors on the dependent axis of a data set. get_exposure([id, bkg_id])Return the exposure time of a PHA data set. get_filter([id])Return the filter expression for a data set. get_fit_contour([id])Return the data used by contour_fit. get_fit_plot([id])Return the data used to create the fit plot. get_fit_results()Return the results of the last fit. get_functions()Return the functions provided by Sherpa. get_grouping([id, bkg_id])Return the grouping array for a PHA data set. get_indep([id, filter, bkg_id])Return the independent axes of a data set. get_int_proj([par, id, otherids, recalc, …])Return the interval-projection object. get_int_unc([par, id, otherids, recalc, …])Return the interval-uncertainty object. get_iter_method_name()Return the name of the iterative fitting scheme. get_iter_method_opt([optname])Return one or all options for the iterative-fitting scheme. get_kernel_contour([id])Return the data used by contour_kernel. get_kernel_image([id])Return the data used by image_kernel. get_kernel_plot([id])Return the data used by plot_kernel. get_method([name])Return an optimization method. get_method_name()Return the name of current Sherpa optimization method. get_method_opt([optname])Return one or all of the options for the current optimization method. get_model([id])Return the model expression for a data set. get_model_autoassign_func()Return the method used to create model component identifiers. get_model_component(name)Returns a model component given its name. get_model_component_image(id[, model])Return the data used by image_model_component. get_model_component_plot(id[, model])Return the data used to create the model-component plot. get_model_contour([id])Return the data used by contour_model. get_model_contour_prefs()Return the preferences for contour_model. get_model_image([id])Return the data used by image_model. get_model_pars(model)Return the names of the parameters of a model. get_model_plot([id])Return the data used to create the model plot. get_model_plot_prefs()Return the preferences for plot_model. get_model_type(model)Describe a model expression. get_num_par([id])Return the number of parameters in a model expression. get_num_par_frozen([id])Return the number of frozen parameters in a model expression. get_num_par_thawed([id])Return the number of thawed parameters in a model expression. get_order_plot([id, orders])Return the data used by plot_order. get_par(par)Return a parameter of a model component. get_pdf_plot()Return the data used to plot the last PDF. get_photon_flux_hist([lo, hi, id, num, …])Return the data displayed by plot_photon_flux. get_pileup_model([id])Return the pile up model for a data set. get_prior(par)Return the prior function for a parameter (MCMC). get_proj()Return the confidence-interval estimation object. get_proj_opt([name])Return one or all of the options for the confidence interval method. get_proj_results()Return the results of the last proj run. get_projection_results()Return the results of the last proj run. get_psf([id])Return the PSF model defined for a data set. get_psf_contour([id])Return the data used by contour_psf. get_psf_image([id])Return the data used by image_psf. get_psf_plot([id])Return the data used by plot_psf. get_pvalue_plot([null_model, alt_model, …])Return the data used by plot_pvalue. get_pvalue_results()Return the data calculated by the last plot_pvalue call. get_quality([id, bkg_id])Return the quality flags for a PHA data set. get_rate([id, filter, bkg_id])Return the count rate of a PHA data set. get_ratio_contour([id])Return the data used by contour_ratio. get_ratio_image([id])Return the data used by image_ratio. get_ratio_plot([id])Return the data used by plot_ratio. get_reg_proj([par0, par1, id, otherids, …])Return the region-projection object. get_reg_unc([par0, par1, id, otherids, …])Return the region-uncertainty object. get_resid_contour([id])Return the data used by contour_resid. get_resid_image([id])Return the data used by image_resid. get_resid_plot([id])Return the data used by plot_resid. get_response([id, bkg_id])Return the response information applied to a PHA data set. get_rmf([id, resp_id, bkg_id])Return the RMF associated with a PHA data set. get_sampler()Return the current MCMC sampler options. get_sampler_name()Return the name of the current MCMC sampler. get_sampler_opt(opt)Return an option of the current MCMC sampler. get_scatter_plot()Return the data used to plot the last scatter plot. get_source([id])Return the source model expression for a data set. get_source_component_image(id[, model])Return the data used by image_source_component. get_source_component_plot(id[, model])Return the data used by plot_source_component. get_source_contour([id])Return the data used by contour_source. get_source_image([id])Return the data used by image_source. get_source_plot([id, lo, hi])Return the data used by plot_source. get_specresp([id, filter, bkg_id])Return the effective area values for a PHA data set. get_split_plot()Return the plot attributes for displays with multiple plots. get_stat([name])Return the fit statisic. get_stat_info()Return the statistic values for the current models. get_stat_name()Return the name of the current fit statistic. get_staterror([id, filter, bkg_id])Return the statistical error on the dependent axis of a data set. get_syserror([id, filter, bkg_id])Return the systematic error on the dependent axis of a data set. get_trace_plot()Return the data used to plot the last trace. group([id, bkg_id])Turn on the grouping for a PHA data set. group_adapt(id[, min, bkg_id, maxLength, …])Adaptively group to a minimum number of counts. group_adapt_snr(id[, min, bkg_id, …])Adaptively group to a minimum signal-to-noise ratio. group_bins(id[, num, bkg_id, tabStops])Group into a fixed number of bins. group_counts(id[, num, bkg_id, maxLength, …])Group into a minimum number of counts per bin. group_snr(id[, snr, bkg_id, maxLength, …])Group into a minimum signal-to-noise ratio. group_width(id[, num, bkg_id, tabStops])Group into a fixed bin width. guess([id, model, limits, values])Estimate the parameter values and ranges given the loaded data. ignore([lo, hi])Exclude data from the fit. ignore2d([val])Exclude a spatial region from all data sets. ignore2d_id(ids[, val])Exclude a spatial region from a data set. ignore2d_image([ids])Exclude pixels using the region defined in the image viewer. ignore_bad([id, bkg_id])Exclude channels marked as bad in a PHA data set. ignore_id(ids[, lo, hi])Exclude data from the fit for a data set. image_close()Close the image viewer. image_data([id, newframe, tile])Display a data set in the image viewer. image_deleteframes()Delete all the frames open in the image viewer. image_fit([id, newframe, tile, deleteframes])Display the data, model, and residuals for a data set in the image viewer. image_getregion([coord])Return the region defined in the image viewer. image_kernel([id, newframe, tile])Display the 2D kernel for a data set in the image viewer. image_model([id, newframe, tile])Display the model for a data set in the image viewer. image_model_component(id[, model, newframe, …])Display a component of the model in the image viewer. image_open()Start the image viewer. image_psf([id, newframe, tile])Display the 2D PSF model for a data set in the image viewer. image_ratio([id, newframe, tile])Display the ratio (data/model) for a data set in the image viewer. image_resid([id, newframe, tile])Display the residuals (data - model) for a data set in the image viewer. image_setregion(reg[, coord])Set the region to display in the image viewer. image_source([id, newframe, tile])Display the source expression for a data set in the image viewer. image_source_component(id[, model, …])Display a component of the source expression in the image viewer. image_xpaget(arg)Return the result of an XPA call to the image viewer. image_xpaset(arg[, data])Return the result of an XPA call to the image viewer. int_proj(par[, id, otherids, replot, fast, …])Calculate and plot the fit statistic versus fit parameter value. int_unc(par[, id, otherids, replot, min, …])Calculate and plot the fit statistic versus fit parameter value. link(par, val)Link a parameter to a value. list_bkg_ids([id])List all the background identifiers for a data set. list_data_ids()List the identifiers for the loaded data sets. list_functions([outfile, clobber])Display the functions provided by Sherpa. list_iter_methods()List the iterative fitting schemes. list_methods()List the optimization methods. list_model_components()List the names of all the model components. list_model_ids()List of all the data sets with a source expression. list_models([show])List the available model types. list_priors()Return the priors set for model parameters, if any. list_response_ids([id, bkg_id])List all the response identifiers of a data set. list_samplers()List the MCMC samplers. list_stats()List the fit statistics. load_arf(id[, arg, resp_id, bkg_id])Load an ARF from a file and add it to a PHA data set. load_arrays(id, *args)Create a data set from array values. load_ascii(id[, filename, ncols, colkeys, …])Load an ASCII file as a data set. load_ascii_with_errors(id[, filename, …])Load an ASCII file with asymmetric errors as a data set. load_bkg(id[, arg, use_errors, bkg_id])Load the background from a file and add it to a PHA data set. load_bkg_arf(id[, arg])Load an ARF from a file and add it to the background of a PHA data set. load_bkg_rmf(id[, arg])Load a RMF from a file and add it to the background of a PHA data set. load_conv(modelname, filename_or_model, …)Load a 1D convolution model. load_data(id[, filename])Load a data set from a file. load_filter(id[, filename, bkg_id, ignore, …])Load the filter array from a file and add to a data set. load_grouping(id[, filename, bkg_id])Load the grouping scheme from a file and add to a PHA data set. load_image(id[, arg, coord, dstype])Load an image as a data set. load_multi_arfs(id, filenames[, resp_ids])Load multiple ARFs for a PHA data set. load_multi_rmfs(id, filenames[, resp_ids])Load multiple RMFs for a PHA data set. load_pha(id[, arg, use_errors])Load a PHA data set. load_psf(modelname, filename_or_model, …)Create a PSF model. load_quality(id[, filename, bkg_id])Load the quality array from a file and add to a PHA data set. load_rmf(id[, arg, resp_id, bkg_id])Load a RMF from a file and add it to a PHA data set. load_staterror(id[, filename, bkg_id])Load the statistical errors from a file. load_syserror(id[, filename, bkg_id])Load the systematic errors from a file. load_table(id[, filename, ncols, colkeys, …])Load a FITS binary file as a data set. load_table_model(modelname, filename[, method])Load tabular or image data and use it as a model component. load_template_interpolator(name, …)Set the template interpolation scheme. load_template_model(modelname, templatefile)Load a set of templates and use it as a model component. load_user_model(func, modelname[, filename])Create a user-defined model. load_user_stat(statname, calc_stat_func[, …])Create a user-defined statistic. load_xstable_model(modelname, filename)Load a XSPEC table model. normal_sample([num, sigma, correlate, id, …])Sample the fit statistic by taking the parameter values from a normal distribution. notice([lo, hi])Include data in the fit. notice2d([val])Include a spatial region of all data sets. notice2d_id(ids[, val])Include a spatial region of a data set. notice2d_image([ids])Include pixels using the region defined in the image viewer. notice_id(ids[, lo, hi])Include data from the fit for a data set. pack_image([id])Convert a data set into an image structure. pack_pha([id])Convert a PHA data set into a file structure. pack_table([id])Convert a data set into a table structure. paramprompt([val])Should the user be asked for the parameter values when creating a model? plot(*args)Create one or more plot types. plot_arf([id, resp_id])Plot the ARF associated with a data set. plot_bkg([id, bkg_id])Plot the background values for a PHA data set. plot_bkg_chisqr([id, bkg_id])Plot the chi-squared value for each point of the background of a PHA data set. plot_bkg_delchi([id, bkg_id])Plot the ratio of residuals to error for the background of a PHA data set. plot_bkg_fit([id, bkg_id])Plot the fit results (data, model) for the background of a PHA data set. plot_bkg_fit_delchi([id, bkg_id, replot, …])Plot the fit results, and the residuals, for the background of a PHA data set. plot_bkg_fit_resid([id, bkg_id, replot, …])Plot the fit results, and the residuals, for the background of a PHA data set. plot_bkg_model([id, bkg_id])Plot the model for the background of a PHA data set. plot_bkg_ratio([id, bkg_id])Plot the ratio of data to model values for the background of a PHA data set. plot_bkg_resid([id, bkg_id])Plot the residual (data-model) values for the background of a PHA data set. plot_bkg_source([id, lo, hi, bkg_id])Plot the model expression for the background of a PHA data set. plot_cdf(points[, name, xlabel, replot, …])Plot the cumulative density function of an array of values. plot_chisqr([id])Plot the chi-squared value for each point in a data set. plot_data([id])Plot the data values. plot_delchi([id])Plot the ratio of residuals to error for a data set. plot_energy_flux([lo, hi, id, num, bins, …])Display the energy flux distribution. plot_fit([id])Plot the fit results (data, model) for a data set. plot_fit_delchi([id, replot, overplot, …])Plot the fit results, and the residuals, for a data set. plot_fit_resid([id, replot, overplot, …])Plot the fit results, and the residuals, for a data set. plot_kernel([id])Plot the 1D kernel applied to a data set. plot_model([id])Plot the model for a data set. plot_model_component(id[, model])Plot a component of the model for a data set. plot_order([id, orders])Plot the model for a data set convolved by the given response. plot_pdf(points[, name, xlabel, bins, …])Plot the probability density function of an array of values. plot_photon_flux([lo, hi, id, num, bins, …])Display the photon flux distribution. plot_psf([id])Plot the 1D PSF model applied to a data set. plot_pvalue(null_model, alt_model[, …])Compute and plot a histogram of likelihood ratios by simulating data. plot_ratio([id])Plot the ratio of data to model for a data set. plot_resid([id])Plot the residuals (data - model) for a data set. plot_scatter(x, y[, name, xlabel, ylabel, …])Create a scatter plot. plot_source([id, lo, hi])Plot the source expression for a data set. plot_source_component(id[, model])Plot a component of the source expression for a data set. plot_trace(points[, name, xlabel, replot, …])Create a trace plot of row number versus value. proj(*args)Estimate parameter confidence intervals using the projection method. projection(*args)Estimate parameter confidence intervals using the projection method. reg_proj(par0, par1[, id, otherids, replot, …])Plot the statistic value as two parameters are varied. reg_unc(par0, par1[, id, otherids, replot, …])Plot the statistic value as two parameters are varied. resample_data([id, niter, seed])Resample data with asymmetric error bars. reset([model, id])Reset the model parameters to their default settings. restore([filename])Load in a Sherpa session from a file. sample_energy_flux([lo, hi, id, num, …])Return the energy flux distribution of a model. sample_flux([modelcomponent, lo, hi, id, …])Return the flux distribution of a model. sample_photon_flux([lo, hi, id, num, …])Return the photon flux distribution of a model. save([filename, clobber])Save the current Sherpa session to a file. save_all([outfile, clobber])Save the information about the current session to a text file. save_arrays(filename, args[, fields, ascii, …])Write a list of arrays to a file. save_data(id[, filename, bkg_id, ascii, clobber])Save the data to a file. save_delchi(id[, filename, bkg_id, ascii, …])Save the ratio of residuals (data-model) to error to a file. save_error(id[, filename, bkg_id, ascii, …])Save the errors to a file. save_filter(id[, filename, bkg_id, ascii, …])Save the filter array to a file. save_grouping(id[, filename, bkg_id, ascii, …])Save the grouping scheme to a file. save_image(id[, filename, ascii, clobber])Save the pixel values of a 2D data set to a file. save_model(id[, filename, bkg_id, ascii, …])Save the model values to a file. save_pha(id[, filename, bkg_id, ascii, clobber])Save a PHA data set to a file. save_quality(id[, filename, bkg_id, ascii, …])Save the quality array to a file. save_resid(id[, filename, bkg_id, ascii, …])Save the residuals (data-model) to a file. save_source(id[, filename, bkg_id, ascii, …])Save the model values to a file. save_staterror(id[, filename, bkg_id, …])Save the statistical errors to a file. save_syserror(id[, filename, bkg_id, ascii, …])Save the systematic errors to a file. save_table(id[, filename, ascii, clobber])Save a data set to a file as a table. set_analysis(id[, quantity, type, factor])Set the units used when fitting and displaying spectral data. set_areascal(id[, area, bkg_id])Change the fractional area factor of a PHA data set. set_arf(id[, arf, resp_id, bkg_id])Set the ARF for use by a PHA data set. set_backscal(id[, backscale, bkg_id])Change the area scaling of a PHA data set. set_bkg(id[, bkg, bkg_id])Set the background for a PHA data set. set_bkg_full_model(id[, model, bkg_id])Define the convolved background model expression for a PHA data set. set_bkg_model(id[, model, bkg_id])Set the background model expression for a PHA data set. set_bkg_source(id[, model, bkg_id])Set the background model expression for a PHA data set. set_conf_opt(name, val)Set an option for the confidence interval method. set_coord(id[, coord])Set the coordinate system to use for image analysis. set_counts(id[, val, bkg_id])Set the dependent axis of a data set. set_covar_opt(name, val)Set an option for the covariance method. set_data(id[, data])Set a data set. set_default_id(id)Set the default data set identifier. set_dep(id[, val, bkg_id])Set the dependent axis of a data set. set_exposure(id[, exptime, bkg_id])Change the exposure time of a PHA data set. set_filter(id[, val, bkg_id, ignore])Set the filter array of a data set. set_full_model(id[, model])Define the convolved model expression for a data set. set_grouping(id[, val, bkg_id])Apply a set of grouping flags to a PHA data set. set_iter_method(meth)Set the iterative-fitting scheme used in the fit. set_iter_method_opt(optname, val)Set an option for the iterative-fitting scheme. set_method(meth)Set the optimization method. set_method_opt(optname, val)Set an option for the current optimization method. set_model(id[, model])Set the source model expression for a data set. set_model_autoassign_func([func])Set the method used to create model component identifiers. set_par(par[, val, min, max, frozen])Set the value, limits, or behavior of a model parameter. set_pileup_model(id[, model])Include a model of the Chandra ACIS pile up when fitting PHA data. set_prior(par, prior)Set the prior function to use with a parameter. set_proj_opt(name, val)Set an option for the projection method. set_psf(id[, psf])Add a PSF model to a data set. set_quality(id[, val, bkg_id])Apply a set of quality flags to a PHA data set. set_rmf(id[, rmf, resp_id, bkg_id])Set the RMF for use by a PHA data set. set_sampler(sampler)Set the MCMC sampler. set_sampler_opt(opt, value)Set an option for the current MCMC sampler. set_source(id[, model])Set the source model expression for a data set. set_stat(stat)Set the statistical method. set_staterror(id[, val, fractional, bkg_id])Set the statistical errors on the dependent axis of a data set. set_syserror(id[, val, fractional, bkg_id])Set the systematic errors on the dependent axis of a data set. set_xlinear([plottype])New plots will display a linear X axis. set_xlog([plottype])New plots will display a logarithmically-scaled X axis. set_ylinear([plottype])New plots will display a linear Y axis. set_ylog([plottype])New plots will display a logarithmically-scaled Y axis. show_all([id, outfile, clobber])Report the current state of the Sherpa session. show_bkg([id, bkg_id, outfile, clobber])Show the details of the PHA background data sets. show_bkg_model([id, bkg_id, outfile, clobber])Display the background model expression used to fit a data set. show_bkg_source([id, bkg_id, outfile, clobber])Display the background model expression for a data set. show_conf([outfile, clobber])Display the results of the last conf evaluation. show_covar([outfile, clobber])Display the results of the last covar evaluation. show_data([id, outfile, clobber])Summarize the available data sets. show_filter([id, outfile, clobber])Show any filters applied to a data set. show_fit([outfile, clobber])Summarize the fit results. show_kernel([id, outfile, clobber])Display any kernel applied to a data set. show_method([outfile, clobber])Display the current optimization method and options. show_model([id, outfile, clobber])Display the model expression used to fit a data set. show_proj([outfile, clobber])Display the results of the last proj evaluation. show_psf([id, outfile, clobber])Display any PSF model applied to a data set. show_source([id, outfile, clobber])Display the source model expression for a data set. show_stat([outfile, clobber])Display the current fit statistic. simulfit([id])Fit a model to one or more data sets. subtract([id])Subtract the background estimate from a data set. t_sample([num, dof, id, otherids, numcores])Sample the fit statistic by taking the parameter values from a Student’s t-distribution. thaw(*args)Allow model parameters to be varied during a fit. ungroup([id, bkg_id])Turn off the grouping for a PHA data set. uniform_sample([num, factor, id, otherids, …])Sample the fit statistic by taking the parameter values from an uniform distribution. unlink(par)Unlink a parameter value. unpack_arf(arg)Create an ARF data structure. unpack_arrays(*args)Create a sherpa data object from arrays of data. unpack_ascii(filename[, ncols, colkeys, …])Unpack an ASCII file into a data structure. unpack_bkg(arg[, use_errors])Create a PHA data structure for a background data set. unpack_data(filename, *args, **kwargs)Create a sherpa data object from a file. unpack_image(arg[, coord, dstype])Create an image data structure. unpack_pha(arg[, use_errors])Create a PHA data structure. unpack_rmf(arg)Create a RMF data structure. unpack_table(filename[, ncols, colkeys, dstype])Unpack a FITS binary file into a data structure. unsubtract([id])Undo any background subtraction for the data set.
The sherpa.ui.utils module¶
The sherpa.ui.utils module contains the
Session class that
provides the data-management code used by the
sherpa.ui module.
Classes
ModelWrapper(session, modeltype[, args, kwargs])Session()
Class Inheritance Diagram¶
The sherpa.astro.ui.utils module¶
The sherpa.astro.ui.utils module contains the Astronomy-specific
Session class that provides the
data-management code used by the sherpa.astro.ui module.
Classes
Session()
Class Inheritance Diagram¶
Bug Reports¶
If you have found a bug in Sherpa please report it. The preferred way is to create a new issue on the Sherpa GitHub issue page; that requires creating a free account on GitHub if you do not have one. For those using Sherpa as part of CIAO, please use the CXC HelpDesk system.
Please include an example that demonstrates the issue that will allow the developers to reproduce and fix the problem. You may be asked to also provide information about your operating system and a full Python stack trace; the Sherpa developers will walk you through obtaining a stack trace if it is necessary.
Contributing to Sherpa development¶
Contributions to Sherpa - whether it be bug reports, documentation updates, or new code - are highly encouraged.
At present we do not have any explicit documentation on how to contribute to Sherpa, but it is similar to other open-source packages such as AstroPy.
The developer documentation is also currently lacking.
Indices and tables¶
Glossary¶
- ARF
- The Ancillary Response Function used to describe the effective area curve of an X-ray telescope; that is, the area of the telescope and detector tabulated as a function of energy. The FITS format used to represent ARFs is defined in the OGIP Calibration Memo CAL/GEN/02-002.
- Astropy
- A community Python library for Astronomy: http://www.astropy.org/.
- ChIPS
- The Chandra Imaging and Plotting System, ChIPS, is one of the plotting backends supported by Sherpa. It is the default choice when using Sherpa as part of CIAO.
- CIAO
- The data reduction and analysis provided by the CXC for users of the Chandra X-ray telescope. Sherpa is provided as part of CIAO, but can also be used separately. The CIAO system is available from http://cxc.harvard.edu/ciao/.
- Crates
- The Input/Output library provided as part of CIAO. It provides read and write access to FITS data files, with speciality support for X-ray data formats that follow OGIP standards (such as ARF and RMF files).
- CXC
- The Chandra X-ray Center.
- DS9
- An external image viewer designed to allow users to interact with gridded data sets (2D and 3D). Is is used by Sherpa to display image data, and is available from http://ds9.si.edu/. It uses the XPA messaging system to communicate with external processes.
- FITS
- The Flexible Image Transport System is a common data format in Astronomy, originally defined to represent imaging data from radio telescopes, but has since been extended to contain a mixture of imaging and tabular data. Information on the various standards related to the format are available from the FITS documentation page at HEASARC.
- HEASARC
- NASA’s High Energy Astrophysics Science Archive Research Center at Goddard Space Flight Center: https://heasarc.gsfc.nasa.gov/.
- matplotlib
- The matplotlib plotting package is is one of the plotting backends supported by Sherpa. More information can be found at http://matplotlib.org/.
- OGIP
- The Office for Guest Investigator Programs (OGIP) was a division of the Laboratory for High Energy Astrophysics at Goddard Space Flight Center. The activities of that group have now become the responsibility of the HEASARC FITS Working Group (HFWG), and supports the use of high-energy astrophysics data through multimission standards and archive access. Of particular note for users of Sherpa are the standard documents produced by this group that define the data formats and standards used by high-energy Astrophysics missions.
- PHA
- The standard file format used to store astronomical X-ray spectral data. The format is defined as part of the OGIP set of standards, in particular OGIP memos OGIP/92-007 and OGIP/92-007a. Confusingly, PHA can also refer to the Pulse Height Amplitude (the amount of charge detected) of an event, which is one of the two channel types that can be found in a PHA format file.
- PSF
- The Point Spread Function. This represents the response of an imaging system to a delta function: e.g. what is the shape that a point source would produce when observed by a system. It is dependent on the optical design of the system but can also be influenced by other factors (e.g. for ground-based observatories the atmosphere can add additional blurring).
- RMF
- The Redistribution Matrix Function used to describe the response of an Astronomical X-ray detector. It is a matrix containing the probability of detecting a photon of a given energy at a given detector channel. The FITS format used to represent RMFs is defined in the OGIP Calibration Memo CAL/GEN/02-002.
- WCS
- The phrase World Coordinate System for an Astronomical data set represents the mapping between the measured position on the detector and a “celestial” coordinate. The most common case is in providing a location on the sky (e.g. in Equatorial or Galactic coordinates) for a given image pixel, but it can also be used to map between row on a spectrograph and the corresponding wavelength of light.
- XPA
- The XPA messaging system is used by DS9 to communicate with external programs. Sherpa uses this functionality to control DS9 - by sending it images to display and retriving any regions a used may have created on the image data. The command-line tools used for this commiunication may be available via the package manager for a particular operating system, such as xpa-tools for Ubuntu, or they can be built from source.
- XSPEC
This can refer to either the X-ray Spectral fitting package, or the models from this package. XSPEC is distributed by HEASARC and its home page is https://heasarc.gsfc.nasa.gov/xanadu/xspec/. Sherpa can be built with support for the models from XSPEC.
Sherpa can be built to use XSPEC versions 12.10.1 (patch level a or later), 12.10.0, 12.9.1, or 12.9.0.
At present there is no developer mailing list for Sherpa.