In the last two notebooks I've created, I've shown how to create a user model (using
add_model rather than
load_user_model) and use the "low-level" interface$^\dagger$ to fit a data set (that is, without using the routines in
sherpa.astro.ui). One thing I haven't done is show how to create plots when using the "low-level" interface, so I thought I'd show how here.
The previous notebooks, which probably need to be skimmed for the following to make any sense, are:
$^\dagger$ we haven't come up with a better name for this API yet.
This was written by Douglas Burke on June 5 2015. This notebook, and others that may be of interest, can be found on GitHub at https://github.com/DougBurke/sherpa-standalone-notebooks.
The information in this document is placed into the Publc Domain. It is not an official product of the Chandra X-ray Center, and I make no guarantee that it is not without bugs or embarassing typos. Please contact me via the GitHub repository or on Twitter - at @doug_burke - if you have any questions.
I have updated the notebook for the Sherpa 4.8.2 release in September 2016; there are no changes to the notebook, but you can now use Python 3.5 as well as 2.7. Actually, there's one minor change to the notebook, but that is due to a change in the IPython/Jupyter code rather than Sherpa: it is no longer necessary to say
import logging logging.getLogger('sherpa').propagate = False
to stop seeing double messages from Sherpa.
When was this notebook last run?
import datetime datetime.datetime.now().strftime("%Y-%m-%d %H:%M")
Let's start with the data set I've been using recently:
from scipy import stats import numpy as np import matplotlib.pyplot as plt %matplotlib inline
np.random.seed(1) k_orig = 1.1 theta_orig = 2.6 ysim = stats.gamma.rvs(a=k_orig, scale=theta_orig, size=1000)
I am going to concentrate on the "CDF" version of the data, but will also use the "PDF" version to show how to plot
Data1DInt data sets. The "CDF" refers to the cumulative distribution function and "PDF" the probability density function.
xcdf = ysim.copy() xcdf.sort() ycdf = np.arange(1, xcdf.size+1) * 1.0 / xcdf.size
(ypdf, edges) = np.histogram(ysim, bins=31, density=False) xlo = edges[:-1] xhi = edges[1:]
The data is stored in objects; here the
from sherpa.data import Data1D, Data1DInt dcdf = Data1D('cdf', xcdf, ycdf) dpdf = Data1DInt('pdf', xlo, xhi, ypdf)
To plot these, I need the
from sherpa.plot import DataPlot
DataPlot object is given the data to plot via the
prepare method, and then actually plotted with the
plot method. Note that
prepare takes an optional statistics object, which is used to determine how to estimate error bars if none are given explicitly. For this example I do not want any error bars, so leave it at its default.
plotcdf = DataPlot() plotcdf.prepare(dcdf)
With all this I can create a plot of the data:
matplotlib commands can be used to adjust the display; for instance, to increase the Y range to better show the behavior of the data at large
plotcdf.plot() plt.ylim(-0.1, 1.1)
The approach for the "binned" data set is the same: call
prepare and then
plot. This time I also change the preferences for the plot -
plot_prefs is a dictionary - so that the Y axis is drawn in a log scale:
plotpdf = DataPlot() plotpdf.prepare(dpdf) plotpdf.plot_prefs['ylog'] = True plotpdf.plot() plt.ylim(0.8, 400)
I use the same models as in the previous notebooks, but this time I cut out the documentation, to save space. Note that as I'm using the low-level interface I do not need to call
add_model to use either the
from sherpa import utils from sherpa.models import model from sherpa.models.parameter import tinyval
def calc_gamma_cdf(x, k, theta): tval = theta * 1.0 kval = k * 1.0 x = np.asarray(x) out = np.zeros_like(x) for i,xi in enumerate(x): out[i] = utils.igam(kval, xi/tval) return out def calc_gamma_pdf(x, k, theta): tval = theta * 1.0 kval = k * 1.0 norm = utils.gamma(kval) * theta**kval return x**(kval-1.0) * np.exp(-x / tval) / norm def calc_gamma_pdf_int(xlo, xhi, k, theta): return calc_gamma_cdf(xhi, k, theta) - calc_gamma_cdf(xlo, k, theta)
class GammaCDF(model.ArithmeticModel): """A Gamma CDF.""" def __init__(self, name='gammacdf'): self.k = model.Parameter(name, 'k', 5, min=tinyval, hard_min=0) self.theta = model.Parameter(name, 'theta', 2, min=tinyval, hard_min=0) model.ArithmeticModel.__init__(self, name, (self.k, self.theta)) @model.modelCacher1d def calc(self, pars, x, *args, **kwargs): (k, theta) = pars if len(args) == 1: x = (x + args) / 2.0 return calc_gamma_cdf(x, k, theta) class GammaPDF(model.ArithmeticModel): """A Gamma PDF.""" def __init__(self, name='gammapdf'): self.k = model.Parameter(name, 'k', 5, min=tinyval, hard_min=0) self.theta = model.Parameter(name, 'theta', 2, min=tinyval, hard_min=0) self.norm = model.Parameter(name, 'norm', 1, min=0) # allow -ve if they realy want it model.ArithmeticModel.__init__(self, name, (self.k, self.theta, self.norm)) @model.modelCacher1d def calc(self, pars, x, *args, **kwargs): (k, theta, norm) = pars if len(args) == 0: return norm * calc_gamma_pdf(x, k, theta) elif len(args) == 1: return norm * calc_gamma_pdf_int(x, args, k, theta) else: raise ValueError("Expected x or xlo/xhi grid.") def guess(self, dep, *args, **kwargs): norm = dep.sum() self.norm.val = dep.sum() self.norm.frozen = True
To plot a model, I need an instance of the model to plot (in this case it's just a single component, but it could be some more-complicated expression).
mcdf = GammaCDF('cdf') mpdf = GammaPDF('pdf')
from sherpa.plot import ModelPlot
The approach is similar to the data case: create an instance of the
ModelPlot class, call
prepare and then
plot. The difference here is that
prepare needs both the data set (to know what the independent axis is) and the model expression.
mplotcdf = ModelPlot()
The help for the method is not particularly illuminating - it is on our list of things to improve! - but at least the argument names are somewhat meaningful!
Help on method prepare in module sherpa.plot: prepare(data, model, stat=None) method of sherpa.plot.ModelPlot instance
mplotcdf.prepare(dcdf, mcdf) mplotcdf.plot()
plot method also accepts optional arguments, and I am going to take advantage of
overplot to combine the data with the model:
Help on method plot in module sherpa.plot: plot(overplot=False, clearwindow=True) method of sherpa.plot.ModelPlot instance
Using this we can see that the current set of model parameters is not a good fit to the data!
plotcdf.plot() mplotcdf.plot(overplot=True) plt.ylim(-0.05, 1.05)
The above plot may also look a bit like the output of the
sherpa.ui.plot_fit command, which is not-too surprising as it essentially is! I'll show below how the "fit"-style plot is created if you do not want to do it manually, like this.
However, before worrying about fits, I'd like to just repeat the above for the "pdf" data set, since it is binned ("integrated" in Sherpa parlance). Note that it uses a log scale for Y axis due to the settings of
mplotpdf = ModelPlot() mplotpdf.prepare(dpdf, mpdf) plotpdf.plot() mplotpdf.plot(overplot=True)
To look at fit results, I want to first fit the data (so that the models are a better approximation to the data). For this we need the
Fit object. I am not going to explain this in depth since it is a repeat of the previous two notebooks, and I really just need this
done so that the models can be plotted!
from sherpa.optmethods import LevMar, NelderMead from sherpa.stats import Cash, LeastSq from sherpa.fit import Fit
For the CDF data, for which I am going to use a least-square statistic, I use the "default" Sherpa optimiser, which is the
LevMar class. For the PDF data set, which I am using a Maximum Likelihood estimator, I use
NelderMead (which is the Nelder-Mead Simplex algorithm), since, although it is often slower than
LevMar, is less-likely to stop in a local minima (this is a statement of my experience, depends on the data and model being fit, and comes with no warranty).
When using the "standard" Sherpa UI for fitting data, only one optimiser and statistic can be in use at a time. This restriction does not hold here, since you have to be explicit in your choices!
fcdf = Fit(dcdf, mcdf, LeastSq(), LevMar()) fpdf = Fit(dpdf, mpdf, Cash(), NelderMead())
Note that creating a
Fit object does not run the fit, so at this time the model parameters
are unchanged. The
Fit class comes with some useful methods, such as calculating the statistic for the current model and data:
More context is provided by the
name = ids = None bkg_ids = None statname = leastsq statval = 230.64654731128954 numpoints = 1000 dof = 998 qval = None rstat = None
name = ids = None bkg_ids = None statname = cash statval = 11016.758241819585 numpoints = 31 dof = 28 qval = None rstat = None
This is the same output as the
ui.calc_stat_info routine returns, except that some concepts - such as data set ids - are not relevant here. As neither of the chosen statistics are based on (i.e. derived from) the Chi Square class (
rstat fields are not filled in.
So, what happens when the data is fit?
<Fit results instance>
<Fit results instance>
As can be seen, the fits have reduced the statistic value for both data sets. The model parameter values, or the return value of the
fit call, can be viewed for more details, but let's look at the data.
To display both model and data in the same plot, I could do it manually, as I did earlier, but I want to try out the
from sherpa.plot import FitPlot
fplotcdf = FitPlot()
Help on method prepare in module sherpa.plot: prepare(dataplot, modelplot) method of sherpa.plot.FitPlot instance
So, prepare takes "plot" objects, which are the
ModelPlot instances I have already used. However, I need to "refresh" the model plot contents, since the model parameter values have changed (it may well be safer to always re-create these objects so you are sure what is being plotted, but for now I want to save some typing):
mplotcdf.prepare(dcdf, mcdf) fplotcdf.prepare(plotcdf, mplotcdf)
fplotcdf.plot() plt.ylim(-0.05, 1.05)
mplotpdf.prepare(dpdf, mpdf) fplotpdf = FitPlot() fplotpdf.prepare(plotpdf, mplotpdf) fplotpdf.plot()
So, the fits have worked!
The residuals could be calculated directly from the
ModelPlot instances, but Sherpa also provides a
ResidPlot class (there's also
DelchiPlot for when you want
from sherpa.plot import ResidPlot
rplotcdf = ResidPlot()
Help on method prepare in module sherpa.plot: prepare(data, model, stat) method of sherpa.plot.ResidPlot instance
Unlike previous plots, the
stat argument is not optional, so I use the same statistic as in the fit (but a different instance). This is only used for calculating the error bars, which I am going to turn off anyway, so the exact choice is unimportant.
rplotcdf.prepare(dcdf, mcdf, LeastSq())
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with leastsq
With the plot prepared, can create it:
As the defaults are not sensible here, I can change them by
prefs = rplotcdf.plot_prefs prefs['yerrorbars'] = False prefs['marker'] = '' prefs['linestyle'] = 'solid' rplotcdf.plot()
Earlier I mentioned that the residual plot could be created manually from the data and model
plots. This is shown below, and only holds if the
prepare method has been called on each
plt.plot(plotcdf.x, plotcdf.y - mplotcdf.y)
[<matplotlib.lines.Line2D at 0x7f7f5c3a4ba8>]
I can use
matpotlib commands to create a combined fit-and-residual plot (this time I
have to change the
clearwindow option to
False so that the subplots are not obliterated
plt.subplot(2,1,1) fplotcdf.plot(clearwindow=False) plt.subplot(2,1,2) rplotcdf.plot(clearwindow=False)
Perhaps unsurprisingly, given the presence of the
sherpa.ui.plot_fit_resid function, Sherpa
provides a class to handle displays with multiple plots:
from sherpa.plot import SplitPlot
splot = SplitPlot()
addplot with a plot object will add that plot in the next available element in the grid (which defaults to a 2-row, 1-column display).
Larger grids can also be created, as shown below, but for that I need to create a residual plot of the "PDF" case:
rplotpdf = ResidPlot() rplotpdf.prepare(dpdf, mpdf, Cash()) rplotpdf.plot_prefs['yerrorbars'] = False
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with cash
Multiple data sets can be displayed. Here I chose to switch to a two-by-two grid
and plot up the data using the
plot method which requires you to specify the column and
grid location of each plot (in this example I could have stuck with
addplot since the order of the plots matches that used by
addplot, but I wanted to show the direct approach here):
splot.reset(rows=2, cols=2) splot.plot(0, 0, fplotcdf) splot.plot(0, 1, fplotpdf) splot.plot(1, 0, rplotcdf) splot.plot(1, 1, rplotpdf)
Note that the two residual plots use different styles because of the different
plot_prefs settings for each object.