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 Xray Center for use in analysing Xray data (both spectral and imaging) from the Chandra Xray telescope, but it is designed to be a generalpurpose package, which can be enhanced with domainspecific tasks (such as Xray 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 OSX (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 onedimensional data or models, one or two dimensional error analysis, and the results of MonteCarlo 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 twodimensional 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 precompiled version of Sherpa¶
Additional useful Python packages include astropy
, matplotlib
,
and ipythonnotebook
.
Using the Anaconda python distribution¶
The Chandra Xray 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 dryrun
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 pytestxvfb. 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 nolonger 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
fftwinclude_dirs=/usr/local/include
fftwlibdirs=/use/local/lib
fftwlibraries=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:
withxspec = 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
withxspec
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:
withxspec = 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 = wcs5.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:
withxspec = 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 = wcs5.16
If the full XSPEC 12.9.x system has been built then use:
withxspec = 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 = wcs5.16
changing
12.9.1
to12.9.0
as appropriate.If the modelonly 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
, andwcs
libraries were installed.Note that XSPEC 12.10.0 introduces a new
enablexsmodelsonly
flag 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 enablexsmodelsonly
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
pytestargs
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 reinstall 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 Xwindows is set up). This can be avoided by installing the X virtualframe 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 verybrief “smoke” test can be run from the commandline 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
commandline
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 onedimensional 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 onedimensional 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 perfectlydescribed 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 precanned 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 backend 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 onedimensional
gaussian provided by the Gauss1D
class  but more complex examples are possible: these
include multiple components,
sharing models between data sets, and
adding userdefined 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.17549e38 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 justaseasy 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 bestfit 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 bestfit 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 leastsquare 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
(LevenbergMarquardt). The latter is often quicker, but less robust,
so we start with it (the optimiser can be changed and the data refit):
>>> from sherpa.optmethods import LevMar
>>> opt = LevMar()
>>> print(opt)
name = levmar
ftol = 1.19209289551e07
xtol = 1.19209289551e07
gtol = 1.19209289551e07
maxfev = None
epsfcn = 1.19209289551e07
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 bestfit 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 bestfit 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 reused, 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 fullwidth
halfmaximum 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.17549e38 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('witherrors', x, y, staterror=dy)
>>> print(de)
name = witherrors
x = Float64[200]
y = Float64[200]
staterror = Float64[200]
syserror = None
The statistic is changed from least squares to
chisquare (Chi2
), to take advantage
of this extra knowledge (i.e. the Chisquare 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 [Qvalue] = 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.17549e38 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.17549e38 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 1sigma (68.2689%) bounds:
Param BestFit 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 threesigma limit from the confidence analysis above,
ahd the dotted line is added to indicate the threesigma
limit above the bestfit 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 singleparameter 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 twodimensional data¶
Sherpa has support for twodimensional 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 irregularlygridded data, the multidimensional
data classes require
that the coordinate arrays and data values are onedimensional.
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 onedimensional case; in this
case the Polynom2D
class is used
to create a loworder 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 reused 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 = 1e2
>>> 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 bestfit
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 Xray Center (CXC) as a general purpose fitting and modeling tool, with specializations for handling Xray 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
sessionbased API
provided by the
sherpa.astro.ui
and sherpa.ui
modules.
These are wrappers around the ObjectOriented
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
multipledimensional data sets, the current classes only
support one and twodimensional 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 multidimensional data sets can be stored as arrays with dimensionality greater than one, the internal representation used by Sherpa is often a flattened  i.e. onedimensional  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 lowlevel 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 recreate 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
BaseData () 
Base class for all data set types 
Data (name, indep, dep[, staterror, syserror]) 
Generic data set 
DataND (name, indep, dep[, staterror, syserror]) 
Base class for Data1D, Data2D, etc. 
Data1D (name, x, y[, staterror, syserror]) 
1D data set 
Data1DInt (name, xlo, xhi, y[, staterror, …]) 
1D integrated data set 
Data2D (name, x0, x1, y[, shape, staterror, …]) 
2D data set 
Data2DInt (name, x0lo, x1lo, x0hi, x1hi, y[, …]) 
2D integrated data set 
DataSimulFit (name, datasets) 
Store multiple 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) 
Class Inheritance Diagram¶
Creating model instances¶
The sherpa.models
and sherpa.astro.models
namespaces provides a collection of one and
twodimensional 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 onedimensional gaussian using the
Gauss1D
class:
>>> g = models.Gauss1D()
>>> print(g)
gauss1d
Param Type Value Min Max Units
     
gauss1d.fwhm thawed 10 1.17549e38 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 lowercase 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.17549e38 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 onedimensional 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.17549e38 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.17549e38 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 (fullwidth halfmaximum)
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.17549435082e38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 10.0
default_min = 1.17549435082e38
default_max = 3.40282346639e+38
>>> h.fwhm.val
10.0
>>> h.fwhm.min
1.1754943508222875e38
>>> h.fwhm.val = 15
>>> print(h.fwhm)
val = 15.0
min = 1.17549435082e38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 15.0
default_min = 1.17549435082e38
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.17549435082e38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 12.0
default_min = 1.17549435082e38
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.17549e38 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.17549435082e38
max = 3.40282346639e+38
units =
frozen = False
link = None
default_val = 12.0
default_min = 1.17549435082e38
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 fullwidthhalfmaxium 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
softlimit ranges
which match the data.
The idea is to move the parameters to values appropriate
for the data, which can avoid unneeded 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.17549e38 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.17549e38 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 readonly).
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 highedge values in the binned case, or use the midpoint  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 lowedge 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 onedimensional polynomial or a twodimensional 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 onedimensional model with a twodimensional 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 onedimensional 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 reexport 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.17549e38 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 morecomplicated model combinations):
>>> print(sim_model)
(sim1 + sim2)
Param Type Value Min Max Units
     
sim1.fwhm thawed 10 1.17549e38 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.17549e38 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.17549e38 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.17549e38 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.52587891e05, 1.00003815e+00, 5.34057617e05])
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.17549e38 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.17549e38 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 bestfit 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 realworld situations!
As there are no errors for the data set, the leastsquare 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.17549e38 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.17549e38 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 pointspread function (PSF) model to modify a twodimensional model to account for blurring due to the instrument (or other sources, such as the atmosphere for groundbased Astronomical data sets), or the redistribution of the counts as a function of energy as modelled by the RMF when analyzing astronomical Xray spectra.
Twodimensional 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 PSFconvolved 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
precalculate 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
wraparound 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 lowlevel 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.73472348e16, 5.20417043e16, 4.33680869e16,
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
lowlevel 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 onedimensional model directly¶
In the following example a onedimensional 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 twodimensional
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
twodimensional 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, reread 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, reread 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 powerlaw 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 reload 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 powerlaw 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 onedimensional 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 = 1e10
>>> 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 = 1e10
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, reread 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, reread 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.124812.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.51838.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.51831.9199,3.23398.2198 Energy (keV)'
When evaluate, over whole 11024 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 onedimensional 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 onedimensional 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 twodimensional 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 onetoone relation, such as when analyzing
astronomical Xray 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 onedimensional 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 floatingpoint 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 onedimensional 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 twodimensional model¶
The twodimensional case is similar to the onedimensional 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 nonintegrated 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 twodimemensional polynomial model from the AstroPy package.
from sherpa.models import model
from astropy.modeling.polynomial import Polynomial2D
__all__ = ('WrapPoly2D', )
class WrapPoly2D(model.RegriddableModel2D):
"""A twodimensional 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]) 
Onedimensional box function. 
Const1D ([name]) 
A constant model for onedimensional data. 
Cos ([name]) 
Onedimensional cosine function. 
Delta1D ([name]) 
Onedimensional delta function. 
Erf ([name]) 
Onedimensional error function. 
Erfc ([name]) 
Onedimensional complementary error function. 
Exp ([name]) 
Onedimensional exponential function. 
Exp10 ([name]) 
Onedimensional exponential function, base 10. 
Gauss1D ([name]) 
Onedimensional gaussian function. 
Integrate1D ([name]) 

Log ([name]) 
Onedimensional natural logarithm function. 
Log10 ([name]) 
Onedimensional logarithm function, base 10. 
LogParabola ([name]) 
Onedimensional logparabolic function. 
NormGauss1D ([name]) 
Onedimensional normalised gaussian function. 
Poisson ([name]) 
Onedimensional Poisson function. 
Polynom1D ([name]) 
Onedimensional polynomial function of order 8. 
PowLaw1D ([name]) 
Onedimensional powerlaw function. 
Scale1D ([name]) 
A constant model for onedimensional data. 
Sin ([name]) 
Onedimensional sine function. 
Sqrt ([name]) 
Onedimensional square root function. 
StepHi1D ([name]) 
Onedimensional step function. 
StepLo1D ([name]) 
Onedimensional step function. 
TableModel ([name]) 

Tan ([name]) 
Onedimensional tan function. 
UserModel ([name, pars]) 
Support for usersupplied models. 
Box2D ([name]) 
Twodimensional box function. 
Const2D ([name]) 
A constant model for twodimensional data. 
Delta2D ([name]) 
Twodimensional delta function. 
Gauss2D ([name]) 
Twodimensional gaussian function. 
Scale2D ([name]) 
A constant model for twodimensional data. 
SigmaGauss2D ([name]) 
Twodimensional gaussian function (varying sigma). 
NormGauss2D ([name]) 
Twodimensional normalised gaussian function. 
Polynom2D ([name]) 
Twodimensional polynomial function. 
Class Inheritance Diagram¶
The sherpa.astro.models module¶
Classes
Atten ([name]) 
Model the attenuation by the InterStellar Medium (ISM). 
BBody ([name]) 
A onedimensional Blackbody model. 
BBodyFreq ([name]) 
A onedimensional Blackbody model (frequency). 
BPL1D ([name]) 
Onedimensional broken powerlaw function. 
Beta1D ([name]) 
Onedimensional beta model function. 
Beta2D ([name]) 
Twodimensional beta model function. 
DeVaucouleurs2D ([name]) 
Twodimensional de Vaucouleurs model. 
Dered ([name]) 
A dereddening model. 
Disk2D ([name]) 
Twodimensional uniform disk model. 
Edge ([name]) 
Photoabsorption edge model. 
HubbleReynolds ([name]) 
Twodimensional HubbleReynolds model. 
JDPileup ([name]) 
A CCD pileup model for the ACIS detectors on Chandra. 
LineBroad ([name]) 
A onedimensional linebroadening profile. 
Lorentz1D ([name]) 
Onedimensional normalized Lorentz model function. 
Lorentz2D ([name]) 
Twodimensional unnormalised Lorentz function. 
MultiResponseSumModel (name[, pars]) 

NormBeta1D ([name]) 
Onedimensional normalized beta model function. 
Schechter ([name]) 
Onedimensional Schecter model function. 
Sersic2D ([name]) 
Twodimensional 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 unbinned onedimensional data sets defined on a wavelength grid, with units of Angstroms. When used with a binned data set the loweredge 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 powerlaw 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]) 
Powerlaw model. 
Recombination ([name]) 
Opticallythin 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 StorchiBergmann 
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 crosssection
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 photoionization 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 efolded broken power law, neutral medium. 
XSbexriv ([name]) 
The XSPEC bexriv model: reflected efolded 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 multitemperature mekal. 
XSc6pmekl ([name]) 
The XSPEC c6pmekl model: differential emission measure using Chebyshev representations with multitemperature mekal. 
XSc6pvmkl ([name]) 
The XSPEC c6pvmkl model: differential emission measure using Chebyshev representations with multitemperature mekal. 
XSc6vmekl ([name]) 
The XSPEC c6vmekl model: differential emission measure using Chebyshev representations with multitemperature mekal. 
XScarbatm ([name]) 
The XSPEC carbatm model: Nonmagnetic carbon atmosphere of a neutron star. 
XScemekl ([name]) 
The XSPEC cemekl model: plasma emission, multitemperature using mekal. 
XScevmkl ([name]) 
The XSPEC cevmkl model: plasma emission, multitemperature 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 blackbodylike spectrum. 
XScompth ([name]) 
The XSPEC compth model: Paolo Coppi’s hybrid (thermal/nonthermal) hot plasma emission models. 
XScplinear ([name]) 
The XSPEC cplinear model: a nonphysical piecewiselinear 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, multiblack 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, powerlaw dependence for T(r). 
XSdiskpn ([name]) 
The XSPEC diskpn model: accretion disk, black hole, black body. 
XSeplogpar ([name]) 
The XSPEC eplogpar model: logparabolic blazar model with nuFnu normalization. 
XSeqpair ([name]) 
The XSPEC eqpair model: Paolo Coppi’s hybrid (thermal/nonthermal) hot plasma emission models. 
XSeqtherm ([name]) 
The XSPEC eqtherm model: Paolo Coppi’s hybrid (thermal/nonthermal) 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 zerotorque inner boundary. 
XSgadem ([name]) 
The XSPEC gadem model: plasma emission, multitemperature with gaussian distribution of emission measure. 
XSgaussian ([name]) 
The XSPEC gaussian model: gaussian line profile. 
XSgnei ([name]) 
The XSPEC gnei model: collisional plasma, nonequilibrium, temperature evolution. 
XSgrad ([name]) 
The XSPEC grad model: accretion disk, Schwarzschild black hole. 
XSgrbm ([name]) 
The XSPEC grbm model: gammaray burst continuum. 
XShatm ([name]) 
The XSPEC hatm model: Nonmagnetic hydrogen atmosphere of a neutron star. 
XSkerrbb ([name]) 
The XSPEC kerrbb model: multitemperature 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 brokenpower law emissivity profile, black hole emission line. 
XSlogpar ([name]) 
The XSPEC logpar model: logparabolic blazar model. 
XSlorentz ([name]) 
The XSPEC lorentz model: lorentz line profile. 
XSmeka ([name]) 
The XSPEC meka model: emission, hot diffuse gas (MeweGronenschild). 
XSmekal ([name]) 
The XSPEC mekal model: emission, hot diffuse gas (MeweKaastraLiedahl). 
XSmkcflow ([name]) 
The XSPEC mkcflow model: cooling flow, mekal. 
XSnei ([name]) 
The XSPEC nei model: collisional plasma, nonequilibrium, constant temperature. 
XSnlapec ([name]) 
The XSPEC nlapec model: continuumonly 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 selfirradiation. 
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 nonmagnetic atmosphere. 
XSnteea ([name]) 
The XSPEC nteea model: nonthermal 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 selfconsistent 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: planeparallel shocked plasma, constant temperature. 
XSraymond ([name]) 
The XSPEC raymond model: emission, hot diffuse gas, RaymondSmith. 
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: nonequilibrium recombining collisional plasma. 
XSsedov ([name]) 
The XSPEC sedov model: sedov model, separate ion/electron temperature. 
XSsirf ([name]) 
The XSPEC sirf model: selfirradiated 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, multitemperature with gaussian distribution of emission measure. 
XSvgnei ([name]) 
The XSPEC vgnei model: collisional plasma, nonequilibrium, temperature evolution. 
XSvmcflow ([name]) 
The XSPEC vmcflow model: cooling flow, mekal. 
XSvmeka ([name]) 
The XSPEC vmeka model: emission, hot diffuse gas (MeweGronenschild). 
XSvmekal ([name]) 
The XSPEC vmekal model: emission, hot diffuse gas (MeweKaastraLiedahl). 
XSvnei ([name]) 
The XSPEC vnei model: collisional plasma, nonequilibrium, 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: planeparallel shocked plasma, constant temperature. 
XSvraymond ([name]) 
The XSPEC vraymond model: emission, hot diffuse gas, RaymondSmith. 
XSvrnei ([name]) 
The XSPEC vrnei model: nonequilibrium 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, nonequilibrium, temperature evolution. 
XSvvnei ([name]) 
The XSPEC vvnei model: collisional plasma, nonequilibrium, constant temperature. 
XSvvnpshock ([name]) 
The XSPEC vvnpshock model: shocked plasma, plane parallel, separate ion, electron temperatures. 
XSvvpshock ([name]) 
The XSPEC vvpshock model: planeparallel shocked plasma, constant temperature. 
XSvvrnei ([name]) 
The XSPEC vvrnei model: nonequilibrium 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: Opticallythin Compton scattering. 
XSconstant ([name]) 
The XSPEC constant model: energyindependent 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 rolloff 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: highenergy 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: powerlaw 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: superexponential 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 crosssections. 
XSwndabs ([name]) 
The XSPEC wndabs model: photoelectric absorption, warm absorber. 
XSxion ([name]) 
The XSPEC xion model: reflected spectrum of photoionized 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: highenergy 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 crosssections. 
XSzwndabs ([name]) 
The XSPEC zwndabs model: photoelectric 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¶
Evaluate a model on a different grid to the requested one.
This is intended to support convolutionstyle models, where the convolved model should be evaluated on a different grid to the data  e.g. a larger grid, since the convolution will account for signal outside the data range  and then be regridded to match the desired grid.
Classes
Axis (lo, hi) 

EvaluationSpace1D ([x, xhi]) 

EvaluationSpace2D ([x, y, xhi, yhi]) 

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, custom_space, requested_space) 

rebin_flux (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 Xrays.
The models in this module include support for instrument models that describe how Xray photons are converted to measurable properties, such as PulseHeight Amplitudes (PHA) or PulseInvariant channels. These ‘responses’ are assumed to follow OGIP standards, such as [1].
References
[1]  OGIP Calibration Memo CAL/GEN/92002, “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 crosssection
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 photoionization 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 XSpec abundance setting or elemental abundance. get_xschatter
()Return the chatter level used by XSpec. get_xscosmo
()Return the XSpec 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 XSpec model library in use. get_xsxsect
()Return the cross sections used by XSpec models. get_xsxset
(name)Return the XSpec model setting. read_xstable_model
(modelname, filename)Create a XSPEC table model. set_xsabund
(abundance)Set the elemental abundances used by XSpec models. set_xschatter
(level)Set the chatter level used by XSpec. set_xscosmo
(h0, q0, l0)Set the cosmological parameters used by XSpec 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 XSpec models. set_xsxset
(name, value)Set a XSpec 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 “bestfit” parameter settings.
A simple example is the leastsquares 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 onedimensional gaussian with normallydistributed 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('statexample', x, y)
The statistic value, along with perbin 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: leastsquare for when there is no knowledge about the errors on each point, a variety of chisquare statistics for when the errors are assumed to be gaussian, and a variety of maximumlikelihood estimators for Poissondistributed 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 usersupplied 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
LevenbergMarquardt
optimiser are:
>>> from sherpa.optmethods.LevMar
>>> lm = LevMar()
>>> print(lm)
name = levmar
ftol = 1.19209289551e07
xtol = 1.19209289551e07
gtol = 1.19209289551e07
maxfev = None
epsfcn = 1.19209289551e07
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, nonisolated 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 physicallymeaningful 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 runtime 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 “singleshot” routines, which start from a guessed set of parameters, and then try to improve the parameters in a continuous fashion, and the “scattershot” 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.
Singleshot techniques¶
As the reader might expect, the singleshot 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 singleshot 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 singleshot 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 singleshot 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 singleshot 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 highdimension 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 welldefined topology. Since it works in a different way than LevenbergMarquardt 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.
LevenbergMarquardt
 This can be considered to be a censored maximumgradients 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 minimumrefiners 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.
Scattershot 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 incrediblytediously 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
method
attribute to one of the singleshot 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 gridsearcher will miss it completely. This is a good technique for finding an approximation to the minimum for a slowlyvarying 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 overkill for relatively simple problems.
Summary and bestbuy strategies¶
Overall, the singleshot methods are best regarded as ways of refining minima located in other ways: from good starting guesses, or from the scattershot techniques. Using intelligence to come up with a good firstguess solution is the best approach, when the singleshot refiners can be used to get accurate values for the parameters at the minimum. However, I would certainly recommend running at least a second singleshot 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  singleshot  fast  OK for refining minima 
LevenbergMarquardt  singleshot  fast  OK for refining minima 
GridSearch  scattershot  slow  OK for smooth functions 
Monte Carlo  scattershot  very slow  Good in many cases 
Reference/API¶
The sherpa.optmethods module¶
Classes
OptMethod (name, optfunc) 
Base class for the optimisers. 
LevMar ([name]) 
LevenbergMarquardt optimization method. 
NelderMead ([name]) 
NelderMead 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
Fit
object; call its
fit()
method one or more times, potentially varying themethod
attribute to change the optimiser; inspect the
FitResults
object returned byfit()
to extract information about the fit; review the fit quality, perhaps refitting with a different set of parameters or using a different optimiser;
 once the “bestfit” solution has been found, calculate error estimates by calling the
est_errors()
method, which returns aErrorEstResults
object detailing the results; visualize the parameter errors with classes such as
IntervalProjection
andRegionProjection
.
The following discussion uses a onedimensional 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 overriden 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 perbin chisquared value (for chisquare
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.46616165292e180
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 chisquare. 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
chisquare 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 bestfit 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 [Qvalue] = 4.67533e20
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.67533320771e20
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 refit.
>>> 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 previouslyapplied 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 (LevenbergMarquardt) 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 [Qvalue] = 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 bestfit 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 [Qvalue] = 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.782733e01
6.000000e+00 4.016824e+00 9.381649e+00 2.416915e+00 4.782733e01
7.000000e+00 4.016833e+00 9.384890e+00 2.416081e+00 4.782733e01
8.000000e+00 4.016879e+00 9.384890e+00 2.416915e+00 4.784385e01
9.000000e+00 4.016822e+00 9.384890e+00 2.416915e+00 4.782733e01
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.7763568394002505e15
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.52e11
1.47e13
8.87e14
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. Realworld examples often require more indepth 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 chisquare 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.7788876325820937e22
This provides evidence that the threeterm 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 1sigma (68.2689%) bounds:
Param BestFit 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 onesigma (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 chisquare distributed
so the change in statistic as a statistic varies can be related to
a confidence bound. One requirement is that  for chisquare
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 onesigma 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.6sigma (89.0401%) bounds:
Param BestFit 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.39950074e02 6.51777887e02]
[ 4.39950074e02 6.29454494e02 6.59925111e04]
[ 6.51777887e02 6.59925111e04 9.77666831e04]]
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 chisquared
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 1sigma (68.2689%) bounds:
Param BestFit 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 realworld
situations the error analysis can often take significantlylonger
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 bestfit
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 replotting, the
calc()
method needs to be rerun to calculate the new data.
The onesigma 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
intervalprojection 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 onesigma 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 onedimensional 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 twodimensional 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 errorestimation 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 (datamodel) 
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 chisquare value per point. 
ComponentModelPlot () 

DelchiPlot () 
Create plots of the deltachi 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 chisquare value per point. 
BkgChisqrPlot () 
Derived class for creating background plots of 1D chi**2 ((datamodel)/error)**2 
DataPlot () 
Create 1D plots of data values. 
BkgDataPlot () 
Derived class for creating plots of background counts 
DelchiPlot () 
Create plots of the deltachi value per point. 
BkgDelchiPlot () 
Derived class for creating background plots of 1D delchi chi ((datamodel)/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 (datamodel) 
Class Inheritance Diagram¶
Markov Chain Monte Carlo and Poisson data¶
Sherpa provides a Markov Chain Monte Carlo (MCMC) method designed for Poissondistributed data. It was originally developed as the Bayesian LowCount Xray 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 Poissondistributed 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 NelderMead 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 errorestimation method instance 'covariance'>
>>> eres = f.est_errors()
>>> print(eres.format())
Confidence Method = covariance
Iterative Fit Method = None
Fitting Method = neldermead
Statistic = cash
covariance 1sigma (68.2689%) bounds:
Param BestFit 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¶
MonteCarlo Markov Chain support for lowcount data (Poisson statistics).
The sherpa.sim
module provides support for exploring the posterior
probability density of parameters in a fit to lowcount data, for
which Poisson statistics hold, using a Bayesian algorithm and a
MonteCarlo Markov Chain (MCMC). It was originally known as the
pyBLoCXS (python Bayesian LowCount Xray 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 AstroStatistics collaboration http://heawww.harvard.edu/AstroStat/, and was called pyBLoCXS. It has since been developed by the Chandra Xray 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 pvalues 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 singlecomponent 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 MonteCarlo 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:
MH
uses a MetropolisHastings jumping rule that is a multivariate tdistribution with userspecified degrees of freedom centered on the bestfit parameters, and with multivariate scale determined by thecovar
function applied at the bestfit location.MetropolisMH
mixes this Metropolis Hastings jumping rule with a Metropolis jumping rule centered at the current draw, in both cases drawing from the same tdistribution as used withMH
. The probability of using the bestfit location as the start of the jump is given by thep_M
parameter of the rule (useget_sampler
orget_sampler_opt
to view andset_sampler_opt
to 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
priorparameter 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 2D 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://heawww.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 Xray 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 bestfit 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 burnin 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 burnin 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 () 
Highlevel 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 LowCount Xray 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 highenergy Xray 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 pvalues 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 singlecomponent spectral models. It should be used with great care in more complex settings. Readers interested in Bayesian lowcount 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 MetropolisHastings jumping rule that is a multivariate tdistribution 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 tdistribution 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://heawww.harvard.edu/AstroStat/pyBLoCXS/
Classes
LimitError 

MH (fcn, sigma, mu, dof, *args) 
The Metropolis Hastings Sampler 
MetropolisMH (fcn, sigma, mu, dof, *args) 
The Metropolis MetropolisHastings Sampler 
Sampler () 

Walk ([sampler, niter]) 
Functions
dmvnorm (x, mu, sigma[, log]) 
Probability Density of a multivariate Normal distribution 
dmvt (x, mu, sigma, dof[, log, norm]) 
Probability Density of a multivariate 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 tdistribution. 
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 ObjectOriented 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 floatingpoint 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 rootfinding 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]) 
Onedimensional 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 onedimensional interpolation. 
multinormal_pdf (x, mu, sigma) 
The PDF of a multivariatenormal. 
multit_pdf (x, mu, sigma, dof) 
The PDF of a multivariate studentt. 
nearest_interp (xout, xin, yin) 
Nearestneighbor onedimensional interpolation. 
neville (xout, xin, yin) 
Polynomial onedimensional interpolation using Neville’s method. 
neville2d (xinterp, yinterp, x, y, fval) 
Polynomial twodimensional 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 onedimensional 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 secondorder 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 secondorder 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 leastsquares
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 leastsquares 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 bestfit 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 secondorder 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 recreated 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 independentaxis 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 bestfit 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 leastsquared 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 onedimensional 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 suboptimal  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 nonAstronomy 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 hardcoded 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
usermodel 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 Chisquare 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 [Qvalue] = 5.64518e19
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 morenoticeable 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 chisquare 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 refit; 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
lesslikely 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 reuse 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 [Qvalue] = 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 divisionbyzero 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 bestfit 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 1sigma (68.2689%) bounds:
Param BestFit 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.6sigma (89.0401%) bounds:
Param BestFit 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 morerobust, but often significantlyslower, 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 confidencederived 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 1sigma (68.2689%) bounds:
Param BestFit 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 onedimensional 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 bestfit location for the parameter, and the horizontal
line the statistic value for the bestfit 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 onesigma 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 Chisquare statistic). The plot shows how, as the parameter moves away from its bestfit 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 bestfit value. As with the intervalprojection 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 threesigma 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 [Qvalue] = 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
bestfit 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 1sigma (68.2689%) bounds:
Param BestFit 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 erroranalysis plots, such as the regionprojection contour plot (where the limits have been chosen to cover the threesigma 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 objectbased 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 Astronomyspecific version adds domainspecific functionality; in this case support for Astronomical data analysis, with a strong focus on highenergy (Xray) 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 automaticallycreated 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 userdefined model class. add_user_pars
(modelname, parnames[, …])Add parameter information to a user model. calc_chisqr
([id])Calculate the perbin chisquared 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 confidenceinterval 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 intervalprojection object. get_int_unc
([par, id, otherids, recalc, …])Return the intervaluncertainty 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 iterativefitting 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 modelcomponent 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 confidenceinterval 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 regionprojection object. get_reg_unc
([par0, par1, id, otherids, …])Return the regionuncertainty 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 userdefined model. load_user_stat
(statname, calc_stat_func[, …])Create a userdefined 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 chisquared 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 (datamodel) 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 (datamodel) 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 iterativefitting scheme used in the fit. set_iter_method_opt
(optname, val)Set an option for the iterativefitting 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 logarithmicallyscaled X axis. set_ylinear
([plottype])New plots will display a linear Y axis. set_ylog
([plottype])New plots will display a logarithmicallyscaled 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 tdistribution. 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 reexports 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 userdefined model class. add_user_pars
(modelname, parnames[, …])Add parameter information to a user model. calc_chisqr
([id])Calculate the perbin chisquared 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