This tutorial shows some basic Sherpa features.
First of all, let's activate the inline matplotlib mode. Sherpa seamlessly uses matplotlib to provide immediate visual feedback to the user. Support for matplotlib in Sherpa requires the matplotlib package to be installed.
The following commands just avoid some unnecessary logging duplication when using Sherpa:
import logging logging.getLogger('sherpa').propagate = 0
Now, let's create a simple synthetic dataset, using numpy: a parabola between x=-5 and x=5, with some randomly generated noise
import numpy as np x = np.arange(-5, 5.1) y = x*x + 23.2 + np.random.normal(size=x.size) e = np.ones(x.size)
Let's import Sherpa:
from sherpa.astro import ui as sherpa
WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available
Depending on how you installed Sherpa, certain special features may be enabled or disabled. Sherpa prints some warning messages when it cannot find some of its modules, as shown above. These warnings are benign. You can refer to the Sherpa documentation to find out what additional features you can enable in Sherpa and how to enable them.
Let's load and plot the data we just created. Notice we are assigning the ID
mydata to the dataset we are loading. We will use this ID to refer to the same dataset in the rest of the tutorial. Sherpa can deal with multiple datasets, fit them simultaneously with the same model, and even link parameters between models. Sherpa can read ASCII table and FITS files (provided the
pyfits package is installed).
sherpa.load_arrays("mydata", x, y, e) sherpa.plot_data("mydata")
We can set the model we want to fit to the data using the
set_model call. There are different ways to instantiate the models: in this case, we just use the string
polynom1d to refer to a 1D polynomial. The name of the model will be
poly, and will be accessible as a Python variable. One can use more object oriented patterns to access and instantiate built-in models. Also, new models can be added by the user as Python functions or from tabular data.
Several Sherpa commands can be used to inspect the model. In this case we just use a simple
polynom1d.poly Param Type Value Min Max Units ----- ---- ----- --- --- ----- poly.c0 thawed 1 -3.40282e+38 3.40282e+38 poly.c1 frozen 0 -3.40282e+38 3.40282e+38 poly.c2 frozen 0 -3.40282e+38 3.40282e+38 poly.c3 frozen 0 -3.40282e+38 3.40282e+38 poly.c4 frozen 0 -3.40282e+38 3.40282e+38 poly.c5 frozen 0 -3.40282e+38 3.40282e+38 poly.c6 frozen 0 -3.40282e+38 3.40282e+38 poly.c7 frozen 0 -3.40282e+38 3.40282e+38 poly.c8 frozen 0 -3.40282e+38 3.40282e+38 poly.offset frozen 0 -3.40282e+38 3.40282e+38
By default, only the first component (the intercept) is thawed, i.e. is free to vary in the fit. This corresponds to a constant function. In order to fit a parabola, we need to thaw the coefficients of the first two orders in the polynomial, as shown below.
We are going to fit the dataset using the default settings. However, Sherpa has a number of optimization algorithms, each configurable by the user, and a number of statistics that can be used to take into account the error and other characteristics of data being fitted.
Dataset = mydata Method = levmar Statistic = chi2 Initial fit statistic = 11908.3 Final fit statistic = 7.5272 at function evaluation 8 Data points = 11 Degrees of freedom = 8 Probability [Q-value] = 0.48096 Reduced statistic = 0.9409 Change in statistic = 11900.8 poly.c0 22.6647 poly.c1 0.0391891 poly.c2 1.00155
Notice that Sherpa used a Levemberg Marquadt minimization strategy, and the $\chi^2$ error function. The best fit values of the parameters $c_0=22.9, c_1=-0.1, c_2=1.1$ are close to the ones defined when we generated the dataset, i.e. $23.2, 0, 1$
In order to get immediate feedback of the fit results, we can plot the fit and the residuals. Again, Sherpa facilitates the creation of the plots, but users can harvest the power of matplotlib directly if they want to.
We can now compute the confidence intervals for the free parameters in the fit:
Dataset = mydata Confidence Method = confidence Iterative Fit Method = None Fitting Method = levmar Statistic = chi2gehrels confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- poly.c0 22.6647 -0.455477 0.455477 poly.c1 0.0391891 -0.0953463 0.0953463 poly.c2 1.00155 -0.0341394 0.0341394 poly.c0 lower bound: -0.455477 poly.c1 lower bound: -0.0953463 poly.c2 lower bound: -0.0341394 poly.c0 upper bound: 0.455477 poly.c1 upper bound: 0.0953463 poly.c2 upper bound: 0.0341394
Sherpa allows to inspect the parameter space. In the cell below we ask sherpa to show us the projection of the confidence regions for the
c1 parameters. The contours are configurable by the user: by default they show the confidence curves at 1, 2, and 3 $\sigma$
We can also directlty inspect the parameter space. For instance, in the plot below Sherpa displays the Interval Projection of the
c0 parameter, i.e. a plot of the error for each value of the parameter, around the minimum found by the optimization method during the fit.