Welcome to the Bayesian X-ray Analysis (BXA) tutorial page.

BXA connects PyXspec/Sherpa to nested sampling algorithms for Bayesian parameter estimation and model comparison with the Bayesian evidence.

If you have any questions/comments/suggestions concerning the tutorial, please get in touch!

Why use BXA?

Key strengths:

  • Systematic fitting with minimal user input of initial parameters
  • Efficient to systematically fit large datasets with many sources or compare many models
  • Provides accurate parameter constraints, even for (very) low counts or for complex parameter spaces with strong degeneracies
For the full list of features and strengths available, see BXA's documentation. Here, we focus on two specific scenarios.

Parameter space constraints in the (very) low count regime

BXA is capable of accurately estimating parameter confidence contours, even when no source counts are present. In this figure, a series of obscured accreting supermassive black hole (ak AGN) spectra are simulated (left panel). The right panel then shows the BXA_derived confidence contours for intrinsic line of sight obscuring column density (log NH) vs. intrinsic (i.e. unabsorbed) X-ray luminosity in the 2-10 keV band.
By simultaneously fitting the background and source spectra, BXA is able to derive a confidence contour for low count spectra even down to zero source counts (though the constraint is very poor, as expected). The contour found makes sense, since an unobscured AGN (i.e. log NH ~ 20) should have been detected if it were very bright. This is why the upper left portion of the contour plane is always ruled out.

No need for fine-tuning initial parameters

BXA does not require a single pre-defined position in the parameter space to start exploring. In this figure, a series of X-ray spectra are simulated at incrementally brighter fluxes. The right panel shows the confidence contour between two parameters in the simulations - black hole spin vs. the height of the X-ray emitting source above the black hole (assuming the lamp post geometry). Both parameters were assigned uninformative priors before fitting over a wide range of possible values. The size and shape of the confidence contours found are thus independent of any single initial parameter values.


⚠️ WARNING! ⚠️ BXA requires a working form of PyXspec and/or Sherpa.

Make sure either of these are installed before doing the exercises! Also note the exercises were written for BXA v4.0.5 - if some commands are not working as expected, check which version you are running.

For PyXspec users

First, install PyXspec: We advise using a conda environment (e.g. with Miniconda) to ensure no conflicts occur when installing the Python packages required for BXA. If you have not created a dedicated BXA conda environment before:
  • Create a new conda environment:
    $ conda create -n xspecbxa python=3.8 
  • Activate the new environment:
    $ conda activate xspecbxa 
  • Install BXA:
    • With Anaconda:
      $ conda install -c conda-forge bxa 
    • With pip (see warning below):
      $ pip install requests corner astropy h5py cython scipy tqdm pandas
      $ pip install bxa
Once the install process has finished, check everthing is installed and works as required: 
  • Check PyXspec works:
    $ python -c 'import xspec'  
  • Check UltraNest works:
    $ python -c 'import ultranest'  
  • Check BXA works:
    $ python -c 'import bxa.xspec as bxa'  
If no error messages come up from these commands, BXA and PyXspec should work.

For Sherpa users

  • First, install CIAO and Sherpa in a conda environment:
  • Once the install process has finished, check Sherpa works:
    $ conda activate ciao-environment-name
    $ sherpa

    If there are no errors, Sherpa should be working as expected.
Then install BXA inside the CIAO conda environment:
  • Activate the environment:
    $ conda activate ciao-environment-name  
  • Install BXA:
    • With Anaconda:
      $ conda install -c conda-forge bxa  
    • With pip (see warning below):
      $ pip install requests corner astropy h5py cython scipy tqdm pandas
      $ pip install bxa
Once the install process has finished, check everthing works as required: 
  • Check UltraNest works:
    $ python -c 'import ultranest'  
  • Check BXA works:
    $ python -c 'import bxa.sherpa as bxa'  
If no error messages come up from these commands, BXA should work with Sherpa.

*⚠️ Warning!
Installing packages with conda and pip can lead to package conflicts, such as:
ValueError: numpy.ndarray size changed, may indicate binary incompatibility. Expected 96 from C header, got 88 from PyObject

If this occurs, from inside your conda environment use:
pip uninstall ultranest
conda install -c cond-forge ultranest
This session aims to provide a sufficient background with Xspec/PyXspec and/or Sherpa to complete the remaining BXA tutorial exercises. For a complete and in-depth understanding of the features available in these software packages, see the useful links for the excellent online material provided by each team.

Key objectives:
  • Load a spectrum and fit a model to it in Xspec, PyXspec and/or Sherpa
  • Manually explore a parameter space and avoid local minima
  • Get to grips with plotting in Xspec, PyXspec and Sherpa

Getting started

If you don't have your own, download the following Chandra/ACIS spectral files from here. The zip file contains:
  • Source spectrum (binned with ftgrouppha grouptype=bmin groupscale=1): sim_chan_bmin.pi
  • Response: chan.rmf
  • Effective area: chan.arf
  • Background spectrum: sim_chan_bkg.pi
Here are some initial steps to set up a fit for BXA in PyXspec or Sherpa:
  • For interactive mode:
    • Xspec:
      $ xspec 
    • PyXspec:
      $ ipython
      from xspec import *
    • Sherpa:
      $ sherpa 

  • For scripting:
    • Xspec:
      $ xspec - xspec_script.tcl
      Remember to include exit or quit at the end of the script!
    • PyXspec:
      $ python3 pyxspec_script.py
    • Sherpa:
      $ sherpa sherpa_script.py
Before starting, the following commands are required:
  • Set the x-axis units (e.g., instrument channels or energy):

  • Plot device (Xspec & PyXspec only):
    • Xspec:
      • Set the plot device: cpd /xw 
    • PyXspec:
      • No plot device (useful to avoid plot window pop-ups): Plot.device = "/null" 
      • Set the plot device: Plot.device = "/xw" 
    (other plot devices are available, see here)

Next, we load in the data and ignore ineffective energies:
  • Load spectra:
    • Xspec:
      data 1:1 sim_chan_bmin1.pi 
    • PyXspec:
      s = Spectrum("sim_chan_bmin1.pi") 
      AllData("1:1 sim_chan_bmin1.pi") 
    • Sherpa:
      load_pha(1, "sim_chan_bmin1.pi") 

  • Ignore ineffective energies:
    • Xspec:
      ignore 1:0.-0.5 8.-** 
    • PyXspec:
      s.ignore("0.-0.5 8.-**") 
      AllData.ignore("1:0.-0.5 8.-**") 
    • Sherpa:
      ignore_id(1, "0.:0.5,8.:") 

There are a number of useful settings for plotting the data in Xspec & PyXspec. For a full list, see here. Some useful examples for these exercises include:
  • Xspec:
    • setplot area/noarea:  
      For plots, divide the observed data by the response effective area
    • setplot add/noadd:  
      Show all additive model components in plots
    • setplot background/nobackground:  
      Automatically plot the background in the same plot as the source
    • setplot rebin X Y N:  
      Adjacent bins are combined in dataset N until either: (1) the significance of the combined bins is ≥ X σ or (2) Y bins have been combined. To undo on all datasets, use setplot rebin 0 1 -1.
    • setplot command range x 0.4 10.  
      To append any IPLOT commands to a list that is executed when the user types plot ... . For a full list of available commands, see:
  • PyXspec (for command descriptions see above, for command documentation see here):
    • Plot.area=True/False  
    • Plot.add=True/False  
    • Plot.background=True/False  
    • Plot.setRebin(X, Y, N)  
    • Plot.addCommand("range x 0.4 10.")  
Plot the observed data inside Xspec, PyXspec and Sherpa:
  • Xspec (note you can try counts/data/lcounts in place of ldata):
    setplot rebin 2 1000 -1
    plot ldata
  • PyXspec (note you can try "counts"/"data"/"lcounts" in place of "ldata"):
    Plot.setRebin(2, 1000)
  • Sherpa (note you can try "counts" in place of "rate"):
    group_snr(1, 2.)
    set_analysis(1, "ener", "rate")
    plot_data(xlog=True, ylog=True)

Interactive PLT (Xspec only)

In Xspec, you can either use plot or iplot to produce plots. plot will make the plot and return the user to the Xspec prompt, whereas iplot starts the interactive plotting interface which can be used to edit a plot.
You can always view the full list of PLT commands being run by Xspec if you increase the verbosity to chatter 25 before typing iplot. This is useful for determining the plotting element numbers for different colours, etc. Here are some basic commands for using the interactive PLT interface:
  • LAB title text: Change the title of the plot to "text"  
  • time off: Remove the time label in the bottom right of the plot  
  • CS N: Change the character size for text on the plot (values of 0-5 are allowed, 1 is default)  
  • viewport xlo ylo xhi yhi: Change the location of the plot in the plotting window  
  • LAB y ylabel: Change the y axis label  
  • LAB x xlabel: Change the x axis label  
  • COL 0 on 1..2: Change the colour of elements 1..2 to 0 (i.e. black - for all available colours, see Fig. 5.1 here)  
  • LW 5 on 1..2: Change the line width of elements 1-2 to 5 (values can be ≥ 1)  
  • LAB 5 POS x y "This text is... \fiITALIC\fn NORMAL\uSUPERSCRIPT\d\dSUBSCRIPT": create a text label.
    • Greek alphabet in order (use \G for capitals):
      \ga \gb \gg \gd \ge \gz \gy \gh \gi \gk \gl \gm \gn \gc, \go, \gp, \gr \gs \gt \gu \gf \gx \gq \gw
  • LAB 5 COL 4 CS 1.4: Change the colour of label 5 to 4 (i.e. blue) and change the character size to 1.4  
IPLOT commands can be used from the Xspec prompt by appending all the commands to a list. Xspec then runs these commands only when you type plot. To interact with this list of commands:
  • setplot command ...: Type this followed by any interactive PLT command to append the command to the list  
  • setplot list: Show the currently appended plot commands  
  • setplot delete all: Remove all plot commands from the list  
A useful tip is to use the command wdata filename.dat inside interactive PLT (or use setplot command wdata filename.dat with Xspec plot) to save everything being plotted by Xspec to filename.dat (make sure to keep track of the order of the columns!). This file can then be read in with e.g., pandas to plot in Python.

Saving the data to be plotted by yourself can be done with the following:
  • Xspec:
    setplot rebin 2 1000 -1
    setplot command wdata ldata_plot.dat
    plot ldata
    Then the data can be read in to e.g., pandas with:
    import pandas as pd
    df = pd.read_csv("ldata_plot.dat", skiprows=3, delim_whitespace=True, names=["EkeV", "EkeV_err", "data", "data_err"])
  • PyXspec:
    import pandas as pd
    Plot.setRebin(2, 1000)
    EkeV = Plot.x()
    EkeV_err = Plot.xErr()
    data = Plot.y()
    data_err = Plot.yErr()
    df = pd.DataFrame(data={"EkeV": EkeV, "EkeV_err": EkeV_err, "data": data, "data_err": data_err})
  • Sherpa:
    import pandas as pd
    group_snr(1, 2.)
    set_analysis(0, "ener", "rate")
    srcdata = get_data_plot(1)
    EkeV = srcdata.x
    EkeV_err = srcdata.xerr
    data = srcdata.y
    data_err = srcdata.yerr
    df = pd.DataFrame(data={"EkeV": EkeV, "EkeV_err": EkeV_err, "data": data, "data_err": data_err})
Save the data being plotted and produce your own spectral plots in the plotting package of your choice. If using Python, here is an example to get started:
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
err_kwargs = {}
err_kwargs["fmt"] = "o"
err_kwargs["label"] = "Folded data"
err_kwargs["mec"] = "royalblue"
err_kwargs["ecolor"] = "royalblue"
err_kwargs["mfc"] = "white"
err_kwargs["markersize"] = 6.
err_kwargs["markeredgewidth"] = 1.
err_kwargs["linewidth"] = 1.
err_kwargs["capthick"] = 0.
err_kwargs["capsize"] = 0.
err_kwargs["zorder"] = -1.
ax.errorbar(EkeV, data, xerr=EkeV_err, yerr=data_err, **err_kwargs)
The count rate that we detect and plot as a spectrum with e.g., plot ldata in Xspec is not the true spectrum, but a form of the true data that has been folded with all the effects arising from detecting the source. This is why fitting is done with forward-folding.
Forward-folding is the process creating a model from some parameters, and then convolving it with the detector response and effective area. The convolved version of the model is then compared with the convolved data by some statistical measure, and is repeated until the parameters are optimised according to the statistical measure.

  • Set the fit statistic (here we use modified Poisson statistics, which are appropriate in relatively low counts - see figure below):
    • Xspec:
      set stat cstat 
    • PyXspec:
      Fit.statMethod = "cstat"  
    • Sherpa:

  • Here the Poisson and Gaussian statistics only become comparable in the high-count regime.

  • (Optional) turn off the interruption whilst fitting (Xspec & PyXspec only):
    • Xspec:
      query yes  
    • PyXspec:
      Fit.query = "yes"  

Then to start with, create a simple powerlaw model and set the model parameter values:
  • Create the model:
    • Xspec:
      model powerlaw 
    • PyXspec:
      mymod = Model("powerlaw") 
    • Sherpa:
      mymod = xspowerlaw.mypow
      set_source(1, mymod)
      Note xs is used to indicate Xspec models, and get_default_id() can be used to show the id that is used if an id is not specified

  • Set the parameter values:
    Note the min, max parameter bounds are used in the BXA priors and should be set to reasonable values.
    • Xspec:
      newpar 1 1.9, 0.1, -1., -1., 4., 4.
      newpar 2 1.e-3, 0.1, 1.e-8, 1.e-8, 1., 1.
    • PyXspec:
      mymod.powerlaw.PhoIndex.values = (1.9, 0.1, -1., -1., 4., 4.)
      mymod.powerlaw.norm.values = (1.e-3, 0.1, 1.e-8, 1.e-8, 1., 1.)
    • Or:
      AllModels(1)(1).values = 1.9, 0.1, -1., -1., 4., 4.
      AllModels(1)(2).values = 1.e-3, 0.1, 1.e-8, 1.e-8, 1., 1.
    • Sherpa:
      set_par(mypow.phoindex, val=1.9, min=-1., max=4.)
      set_par(mypow.norm, val=1.e-3, min=1.e-8, max=1.)
Then we plot the data + model and fit (try to always plot before you fit! - this helps to see if everything is in order):
  • Xspec:
    plot ldata del
  • PyXspec:
    Plot("ldata del")
  • Sherpa (note if you subtracted the data previously with subtract, you need to run unsubtract if fitting with wstat or cstat):
    plot_fit(xlog=True, ylog=True)

Beware of local minima!

Depending on where you start your fit, the Levenberg-Marquadt algorithm may struggle to find the global minimum (i.e. the true best-fit).

In this example, both plots show an example progression for a fit with the Levenberg-Marquadt algorithm in Xspec. The final fit in the left panel is the global minimum, whereas the right hand plot is a local minimum. Local minima can be quite hard to identify visually and even masquerade as a very good fit. A quick way to initially check is to look for physically-unrealistic parameter values from a fit. The more thorough way to combat this issue whilst using Levenberg-Marquadt is to try to explore the parameter space more fully - e.g., with steppar in Xspec/PyXspec or by manually starting the fit from different initial parameter values.

Have you found the best-fit parameter values for your powerlaw fit? Try starting the parameters at different initial values before fitting to check.

Plotting a model with data

When plotting a model with data, there are a few options inside Xspec, PyXspec and Sherpa:
  • Folded units:
    See the commands from Exercise 0.1
  • Unfolded units:
    Note that plotting in unfolded units is useful to see the possible unconvolved spectrum, but this is model dependent. Changing the model will change the plot.
    • Xspec:
      • plot ufspec  
      • plot eufspec  
      • plot eeufspec  
    • PyXspec:
      • Plot("ufspec")  
      • Plot("eufspec")  
      • Plot("eeufspec")  
      • Note the model y-axis values can be accessed after running the Plot(...) command with e.g., mod_data = Plot.model().
    • Sherpa:
      See Sherpa thread here.
  • Residuals:
    • Xspec:
      • plot residuals  
        Plot the data minus the folded model.
      • plot ratio  
        Plot the data divided by the folded model.
      • plot delchi  
        For cstat, this plots (data-model)/error where error is the square root of the model-predicted number of counts.
    • PyXspec (see equivalent descriptions above):
      • Plot("residuals")  
      • Plot("ratio")  
      • Plot("delchi")  
    • Sherpa:
  • Plot your model fit and data with a residual method of your choice. What can you conclude?
Lastly, we grid over the parameter space to check for local minima:
  • Xspec:
    steppar log/nolog parno1 lowerbound1 upperbound1 nsteps log/nolog parno2 lowerbound2 upperbound2 nsteps
    plot cont
  • PyXspec:
    Fit.steppar("log/nolog parno1 lowerbound1 upperbound1 nsteps log/nolog parno2 lowerbound2 upperbound2 nsteps")
  • Note: you will have to choose the lowerbound and upperbound in steppar to encompass the parameter values.
  • Sherpa:
    See Region Projection in this excellent Sherpa guide.

Key objectives:

  • Fit a spectrum with BXA in Sherpa and/or PyXspec
  • Use Bayes factors to perform model comparison
  • Visualise the confidence ranges on each parameter

To start with, make sure you import the relevant packages:
  • PyXspec:
    from xspec import *
    import bxa.xspec as bxa
  • Sherpa:
    import bxa.sherpa as bxa

Exercise 1.1 - prior predictive checks

Prior predictive checks are useful to see if the priors chosen can lead to realistic model realisations.
Note for these exercises, setting priors in BXA is done with the following commands:
Note: log uniform priors cannot be zero or negative. Make sure the parameter ranges you supply are positive.

To start these exercises, load the data, notice the correct energies and create a powerlaw model (see the commands listed above in Ex 0.1 for help). Then:
  • Set the parameter priors for the powerlaw model (note norm varies over many orders of magnitude, and PhoIndex can be assigned a uniform prior for now):
  • A note on log parameters (Sherpa only): Using a loguniform prior may not give the posterior values of that parameter in log units in Sherpa.
    Instead, the Parameter module can convert a parameter to log units, and then one can use a uniform prior for that parameter. E.g.:
    from sherpa.models.parameter import Parameter
    logpar = Parameter(modelname, parname, value, min, max, hard_min, hard_max)
    mymod.norm = 10**logpar
  • Note the parameter will only be valid for working in BXA. If you want to define a log parameter and use this with a standard Sherpa fit, the log parameter must appear in the model expression.
  • Create a solver:
    • PyXspec:
      solver = bxa.BXASolver(transformations=[prior1, prior2], outputfiles_basename="powerlaw_pyxspec")  
    • Sherpa:
      solver = bxa.BXASolver(prior=bxa.create_prior_function([prior1, prior2]),
                             parameters=[param1, param2],
    Make sure:
    1. outputfiles_basename is unique if running PyXspec and Sherpa in the same directory.
    2. (For Sherpa) the priors are in the same order for the prior and parameters attributes.

  • Perform prior predictive checks:
    1. Generate random samples of prior parameter values (remember to import numpy!):
      • PyXspec:
        values = solver.prior_function(numpy.random.uniform(size=len(solver.paramnames)))  
      • Sherpa:
        values = solver.prior_transform(numpy.random.uniform(size=len(solver.paramnames)))  
    2. Set the parameters to these values:
      • PyXspec:
        from bxa.xspec.solver import set_parameters
      • Sherpa:
        for i, p in enumerate(solver.parameters): p.val = values[i]  
    3. Plot the resulting prior model samples (with the data to help guide the eye):
      • PyXspec:
      • Sherpa:
        plot_fit(xlog=True, ylog=True)  
  • Produce 10 realisations of the model priors. Are the model ranges covered by the priors physically-acceptable?
Now that we have checked the priors look reasonable for the data we want to fit, we can fit the model to the data with BXA.
  • First fit the model using Levenberg-Marquadt included in Xspec & Sherpa:
    • PyXspec:
    • Sherpa:

  • Is the fit good? How can you tell? What does changing the parameters by a small amount do?

  • Next perform the fit with BXA (note this does not require fitting with Levenberg-Marquadt first):
    • PyXspec:
      results = solver.run(resume=True)  
    • Sherpa:
      results = solver.run(resume=True)  
    (Note the hyperlinks are different).
    results contains the posterior samples (results["samples"]) for each parameter and the Bayesian evidence (results["logZ"], results["logZerr"]). These quantities should also be printed to the screen after the fit is completed.
    To understand more about the information UltraNest prints to the screen during and after a fit, see here.

  • Examine the best-fit values and uncertainties from the posterior samples. This can be done with pandas (may require installing with e.g., conda install pandas but see warning at the end of Setup):
    import pandas as pd
    df = pd.DataFrame(data=results["samples"], columns=solver.paramnames)

  • How do the best-fits and confidence ranges compare with using nested sampling vs. Levenberg-Marquadt?

Exercise 1.3 - fit and compare multiple models

Here we use the Bayesian evidence to perform model comparison with Bayes factors.
  • Note to use TBabs, one should set the abundances to those from Wilms, Allen & McCray (2000):
  • Using your previous code, fit the following additional 2 models to the data with BXA (don't forget to give unique outputfiles_basename values!):
    1. Absorbed powerlaw:
      For nh, you can assume delta = 0.1, min = 1.e-2, bottom = 1.e-2, top = 5.e2, max = 5.e2.
    2. Absorbed powerlaw + Gaussian emission line (assume the emission line for this source is known to have line centroid energy LineE = 6.3 ± 0.15 keV, and that the line is narrow with fixed width Sigma = 1.e-3 keV):

  • Copy the model_compare.py script to your working directory to perform Bayesian evidence model comparison

  • Run the script:
    $ python3 model_compare.py outputfiles_basename1/ outputfiles_basename2/ outputfiles_basename3/  

  • Compute the Akaike Information Criterion (AIC) (given by AIC = 2×nfree+Cstat ≡ 2×nfree-2×log(Likelihood), where nfree is the number of free fit parameters) for each model in python:
    import json
    basenames = [outputfiles_basename1, outputfiles_basename2, outputfiles_basename3]
    loglikes = dict([(f, json.load(open(f + "/info/results.json"))["maximum_likelihood"]["logl"]) for f in basenames])
    nfree = dict([(f, len(json.load(open(f + "/info/results.json"))["paramnames"])) for f in basenames])
    aic = dict([(basename, (2. * nfree[basename] - 2. * loglike)) for basename, loglike in loglikes.items()])

  • Which model best explains the data? Why?

  • Try adding a new additive apec component to the model:
    • PyXspec:
      mod4 = Model("TBabs * powerlaw + gaussian + apec")
      mod4.apec.kT.values = 0.8, 0.01, 1.e-2, 1.e-2, 5., 5.
      mod4.apec.Abundanc.values = 1., -1
      mod4.apec.Redshift.values = 0., -1
      mod4.apec.norm.values = 1.e-4, 0.1, 1.e-10, 1.e-10, 1.e-1, 1.e-1
    • Sherpa:
      xstbabs.tbabs * xspowerlaw.powerlaw + xsgaussian.gaussian + xsapec.apec
      set_par(apec.kt, val=0.8, min=1.e-2, max=5.)
      set_par(apec.abundanc, val=1., frozen=True)
      set_par(apec.redshift, val=0., frozen=True)
      set_par(apec.norm, val=1.e-4, min=1.e-10, max=1.e-1)

  • How does the Bayesian Evidence from fitting the new apec model to the data compare to the previous models?

Exercise 1.4 - error propagation

By generating quantities directly from the posterior samples, we propagate the uncertainties and conserve any structure (e.g., degeneracies, multiple modes).
As an example, we will estimate the equivalent width of the emission line used in the Absorbed powerlaw + Gaussian model above.
  • The equivalent width (EW below) of an emission line is defined as the flux in the line divided by the flux of the continuum at the line energy (see the Xspec documentation for more information).

  • For the model we have used, the equivalent width of the emission line relative to the intrinsic (i.e. unabsorbed) powerlaw continuum is (in terms of the parameters defined):
    • PyXspec:
      solver3.set_paramnames(["log(nH)", "mypow.phoindex", "mypow.lognorm", "line.linee", "line.lognorm"])
      df = pd.DataFrame(data = results3["samples"], columns = solver3.paramnames)
      df.loc[:, "EW"] = (10**df["line.lognorm"] / 10**df["mypow.lognorm"]) * df["line.linee"] ** df["mypow.phoindex"]
    • Sherpa:
      df = pd.DataFrame(data = results3["samples"], columns = solver3.paramnames)
      df.loc[:, "EW"] = (df["line.norm"] / df["mypow.norm"]) * df["line.linee"] ** df["mypow.phoindex"]
    Generate posterior samples for equivalent width. If using pandas, this will involve creating a new column in the pandas dataframe of posterior samples for the equivalent width.

  • Use the results["samples"] to derive the posterior distribution on equivalent width for the emission line component. If using pandas:
    import corner
    figure = corner.corner(df, labels = df.columns, quantiles = [0.16, 0.5, 0.84], show_titles = True)

  • Is the equivalent width degenerate with any other parameters? Why or why not?
We want a simple way to check that our model makes sense, that the fit worked, and to infer information on the source being fit.

  • First, examine the corner and trace plots that BXA produces (located in outputfiles_basename/plots/). What do these plots show?

  • Are there any parameter degeneracies? Why or why not?

  • Create a Quantile-Quantile (also here for a more general explanation) plot for your fits:
    • PyXspec (produces plot):
      Warning: if you have used Plot.setRebin(), you should undo this before using bxa.qq with e.g., Plot.setRebin(0., 1.) (or restart without running Plot.setRebin())
      import matplotlib.pyplot as plt
      print("creating quantile-quantile plot ...")
      with bxa.XSilence():
          bxa.qq.qq(prefix=outputfiles_basename, markers=5, annotate=True)
      print("saving plot...")
      plt.savefig(outputfiles_basename + "qq_model_deviations.pdf", bbox_inches="tight")
    • Sherpa (produces file that needs to be plotted separately):
      from bxa.sherpa.qq import qq_export
      bxa.qq.qq_export(1, bkg=False, outfile=outputfiles_basename + "_qq", elow=0.2, ehigh=10)

  • What does the powerlaw model Q-Q plot show? Is this model lacking any components?

  • (PyXspec only) Produce a plot of posterior model spectra (convolved and unconvolved) using the code in the example_simplest.py script included in BXA. Which energies is the model best-constrained over?

Exercise 1.6 - fitting multiple datasets simultaneously with BXA

Sometimes it is advantageous to fit multiple spectra for the same source simultaneously with e.g., a cross-calibration constant to account for instrumental differences in the separate spectra.

For this example, we need another dataset. As an example, we will use a NuSTAR/FPMA spectrum which has a bandpass that extends to E > 10 keV. So including this data with our Chandra spectral analysis will extend the energy range we can study our source over to ~78 keV. To download the NuSTAR spectral files:
  • Source spectrum (sim_nuA_bmin1.pi) & background spectrum (sim_nuA_bkg.pi) here.
  • Response (nustar.rmf) and effective area (point_30arcsecRad_1arcminOA.arf) files from the Point source simulation files (distributed by the NuSTAR team, more info here).

We first assume that the source was not variable between the Chandra and NuSTAR observations, and thus only account for differences arising from the calibration of either instrument.
  • First, include a cross-calibration constant in your model between the Chandra and NuSTAR datasets:

  • Fit the new model with BXA and derive a posterior confidence range on the cross-calibration constant:
    • PyXspec:
      prior1 = bxa.create_gaussian_prior_for(mymod, mymod.powerlaw.PhoIndex, 1.9, 0.15)
      prior2 = bxa.create_loguniform_prior_for(mymod, mymod.powerlaw.norm)
      transformations = [prior1, prior2]
      transformations.append(bxa.create_loguniform_prior_for(AllModels(2), AllModels(2).constant.factor))
      solver = bxa.BXASolver(transformations=transformations, outputfiles_basename="xcal_pl_pyxspec")
      results = solver.run(resume=True)
    • Sherpa:
      prior1 = bxa.create_gaussian_prior_for(mymod, mymod.powerlaw.PhoIndex, 1.9, 0.15)
      prior2 = bxa.create_loguniform_prior_for(mymod, mymod.powerlaw.norm)
      priors = [prior1, prior2]
      solver = bxa.BXASolver(prior=bxa.create_prior_function([prior1, prior2, ...]), parameters=[param1, param2, ...], outputfiles_basename="xcal_pl_sherpa")
      results = solver.run(resume=True)
    Note: if using mod4 from the end of Exercise 1.3, the fit may take a long time to complete (up to an hour).

  • How have the confidence ranges on all the parameters changed after including NuSTAR?

  • How have the Bayes factors between different model fits changed? Why?

So far, the exercises have used modified C-statistics (aka W-statistics) in which the background is modelled as a stepwise function with the same number of parameters as bins. This process typically requires some form of minimal binning (see discussion here).

Using C-statistics can be more flexible, since no binning is required by simultaneously fitting a background model instead. Here we use the auto_background function in BXA, which provides a powerful way to implement background models in Sherpa, and simultaneously fit with a source model using BXA.

  • Load the data, source model and priors as in Exercise 1.1 (make sure to also change the statistic from wstat to cstat!)

  • Load the automatic PCA background model (see Simmonds et al., 2018 for more information):
    from bxa.sherpa.background.pca import auto_background
    set_model(1, model)
    convmodel = get_model(1)
    bkg_model = auto_background(1)
    set_full_model(1, convmodel + bkg_model*get_bkg_scale(1))
  • Remember to include the zeroth order component from the PCA background model in the priors and parameters list so that BXA knows to vary it as a free parameter:
    parameters += [bkg_model.pars[0]]
    priors += [bxa.create_uniform_prior_for(bkg_model.pars[0])]

  • Fit the source + background model with BXA

  • Check the background fit (e.g., with get_bkg_fit, get_bkg_fit_ratio, get_bkg_fit_resid).

  • Have the uncertainty ranges changed at all for the parameters? Why or why not?
Key objectives:
  • Combine individual posteriors to constrain a population
  • Calibrate model selection thresholds with simulations
  • Perform model verification to test the Goodness-of-Fit
These exercises involve simulating spectra. To do this you will need a model, response file (rmf), effective area file (arf) and (optionally) a background spectrum file:
  • PyXspec (using fakeit):
    mymod = Model("powerlaw")
    mymod.powerlaw.PhoIndex.values = (1.9, 0.01, -3., -3., 3., 3.)
    mymod.powerlaw.norm.values = (1.e-3, 0.01, 1.e-8, 1.e-8, 1., 1.)
    fakeit_kwargs = {}
    fakeit_kwargs["response"] = "chan.rmf"
    fakeit_kwargs["arf"] = "chan.arf"
    fakeit_kwargs["background"] = "sim_chan_bkg.pi"
    fakeit_kwargs["exposure"] = 30.e3
    fakeit_kwargs["correction"] = "1."
    fakeit_kwargs["backExposure"] = 30.e3
    fakeit_kwargs["fileName"] = "sim_chan.pi"
    AllData.fakeit(1, FakeitSettings(**fakeit_kwargs))
  • Sherpa (using fake_pha):
    mymod = xspowerlaw.mypow
    set_par(mymod.phoindex, val = 1.9, min = -3., max = 3.)
    set_par(mymod.norm, val = 1.e-3, min = 1.e-8, max = 1.)
    load_pha(1, "chan_src.pi") #loading pre-exising source spectrum with rmf, arf & bkg present
    set_source(1, mymod)
    fakepha_kwargs = {}
    fakepha_kwargs["rmf"] = get_rmf()
    fakepha_kwargs["arf"] = get_arf()
    fakepha_kwargs["bkg"] = get_bkg()
    fakepha_kwargs["exposure"] = 30.e3
    fake_pha(1, **fakepha_kwargs)
    save_pha(1, "sim_chan.pi", clobber = True)

Exercise 2.1 - Bayesian Hierarchical Modelling introduction

Here we will combine several individual source posterior distributions to derive a sample distribution for the photon index .

  • Install PosteriorStacker:
    $ conda install -c conda-forge posteriorstacker  
    Alternatively, you can just download the script and use it in your directory.

  • To generate the posterior samples used in this exercise, download and run the following PyXspec ($ python3 ex21_pyxspec.py ) or Sherpa ($ sherpa ex21_sherpa.py ) script. The script will generate 30 UltraNest output folders by:
    1. Simulating NuSTAR spectra from an absorbed powerlaw model (see Exercise 1.3), with PhoIndex sampled from a Gaussian distribution 1.8 +/- 0.3. Note you will need the Point source simulation files distributed by the NuSTAR team to run the scripts (more info here).
    2. Each spectrum will then be fit with BXA using non-informative priors for line-of-sight column density, powerlaw photon index and powerlaw normalisation.
  • The script may take some time (~ a few hours) to complete. Alternatively, you can download some pre-generated UltraNest folders using a similar model here.

  • From inside the ultranest_posteriors directory, run:
    $ load_ultranest_outputs.py fitsim0/ fitsim1/ fitsim2/ fitsim3/ fitsim4/ fitsim5/ fitsim6/ ... \
                                --samples 1000 \
                                --parameter PhoIndex \
                                --out posterior_samples.txt
  • The posterior_samples.txt file contains 1000 sampled rows of each posterior.

  • Important: if the posterior samples were obtained with a non-uniform prior, the posterior samples should first be resampled according to the inverse of the prior that was used.

  • Visualise the data by calculating the 16th, 50th & 84th quantiles of each posterior sample and plotting as an errorbar. E.g., in Python:
    import numpy as np
    import pandas as pd
    quantiles = [16, 50, 84]
    x = np.loadtxt("posterior_samples.txt")
    q = np.percentile(x, quantiles, axis = 1)
    df = pd.DataFrame(data = {"q%d" %(qvalue): q[i] for i, qvalue in enumerate(quantiles)})
    import matplotlib.pyplot as plt
    plt.errorbar(x=df["q50"], xerr=[df["q50"]-df["q16"], df["q84"]-df["q50"]], y=range(len(df)), marker="o", ls=" ", color="orange")
    plt.xlim(df["q16"].min() - 0.5, df["q84"].max() + 0.5)
  • Run PosteriorStacker:
    $ posteriorstacker.py posterior_samples.txt 0.5 3. 10 --name="PhoIndex"  
  • PosteriorStacker then fits two models for the parent sample distribution: a histogram model with each bin height as the free parameters (using a Dirichlet prior) and a Gaussian model with mean and sigma as the free parameters.

  • Examine the posteriorsamples.txt_out.pdf file.

  • What can you say about the distribution of PhoIndex? What would happen if you acquired more posterior samples (i.e. by increasing population_size in the simulation scripts) and/or used higher signal-to-noise spectra to fit?

Exercise 2.2 - calibrating model comparison thresholds

Here we will estimate the false-positive and false-negative rates associated with selecting models via the Bayes Factors (logZmodel 1/Zmodel 2) & Δ fit-statistic for the models and data we are using. This helps us determine a good selection threshold to use for these model selection methods.

The underlying theme of this exercise is null hypothesis significance testing. We try to find a distribution of model comparison value improvements that occur by chance (the red distribution in this figure). If it is much larger (e.g., than 99.9% of cases), the other model is likely to be correct.

Here, you will estimate the false-negative and false-positive rates for model selection between the powerlaw and absorbed powerlaw models from Exercises 1.2 and 1.3.

First, estimate the false-negative rate:
  • Simulate 100 Chandra spectra from the absorbed powerlaw (i.e. the more complex) best-fit model that you found in Exercise 1.3, with the original Chandra spectrum files.
  • Fit both the powerlaw and absorbed powerlaw models to each of the 100 simulated spectra with BXA.
  • For a series of different values of Δ logZ and Δ W-stat, calculate how many of the 100 model fits were correctly identified, and how many times the Bayes Factors did not make a decision (i.e. logZ1/Z2 < threshold).
  • Use the same process to estimate an appropriate Δ W-stat threshold that will have a small false-negative rate.
  • Repeat these step with a powerlaw normalisation 5× fainter than the best-fit.
  • How does the false-negative rate vary with signal-to-noise in the source spectrum?

Then, estimate the false-positive rate:
  • Simulate 100 Chandra spectra from the powerlaw (i.e. the simplest) best-fit model that you found in Exercise 1.2, with the original Chandra spectrum files.
  • Fit both the powerlaw and absorbed powerlaw models to each of the 100 simulated spectra with BXA.
  • For a series of different values of Δ logZ and Δ W-stat, calculate how many of the 100 model fits selected a best-fit model that was not powerlaw (i.e. a model was selected that was not the model being simulated).
  • Use the same process to estimate an appropriate Δ W-stat threshold that will have a small false-positive rate.
  • Repeat these step with a powerlaw normalisation 5× fainter than the best-fit.
  • How does the false-positive rate vary with signal-to-noise in the source spectrum?
Finally do the same thing for the more complex model:
  • Repeat the process to estimate the false-negative and false-positive rates between powerlaw and absorbed powerlaw + gaussian models.
  • Are the rates comparable? Why or why not?

Exercise 2.3 - posterior predictive checks

Model selection can only tell us how likely a model is relative to another. Here we will estimate the goodness-of-fit and test the ability of a given model to reproduce the nuances in the source spectrum data with Monte Carlo simulations.

  • Load the BXA fit (with e.g., results = solver.run(resume=True) ) for the powerlaw model from Exercise 1.2.

  • Note the best-fit C-statistic & degrees of freedom (dofs):
    stat = Fit.statistic
    dof = Fit.dof
    statinfo = get_stat_info()
    stat = statinfo[0].statval
    dof = statinfo[0].dof

  • Perform many times:
    1. Simulate the best-fit model (make sure to use the same spectral files and exposure time as the original).
    2. Without re-fitting*, store the simulated C-statistic & dofs of the simulated data and model.
      * It would be ideal to re-fit, but this can be inefficient and/or lead to incorrect results with a standard (i.e. Levenberg-Marquadt) fit in PyXspec and/or Sherpa due to local minima. Note also that to perform a Levenberg-Marquadt fit in Sherpa will require all free parameters, including specially-defined log parameters, to appear in the model expression.

  • Plot the distribution of simulated C-statistic / dof values and compare to the real value. Is the fit acceptable?

  • Now repeat the process with the absorbed powerlaw + gaussian model. How does the goodness-of-fit for this model compare to the powerlaw model?

Note this form of posterior predictive checks is similar to the goodness command in Xspec. An alternative would be to run the same simulations, but plot the range of fluxes simulated for every spectral energy bin instead.

Exercise 2.4 (extension) - Bayesian Hierarchical Modelling continued

Now you have created a parent population for photon index using PosteriorStacker, this exercise will focus on the NH distribution.

Context: The column density (NH) distribution is a useful quantity for studies of accreting supermassive black holes (aka Active Galactic Nuclei; AGN). Of particular interest in this regard is the "Compton-thick" fraction - i.e. the fraction of sources with log NH / cm-2 > 24.

The goal of this exercise is to estimate the uncertainty on each bin in a log NH distribution from a survey of 30 AGN.

For the AGN being simulated, we can assume:
  • Every AGN follows the same baseline model:
    model = TBabs[Gal] * TBabs[AGN] * powerlaw
    • TBabs[Gal] is the absorbing column density arising from the Milky Way, assumed to be NH = 1020 cm-2 for every AGN.
    • TBabs[AGN] is the absorbing column density at the source. This is the parameter we want to generate a population NH distribution from.
    • powerlaw is the intrinsic emission from the AGN, with photon index within 1.9 ± 0.1 and normalisation within 10-4 to 10-1 photons/keV/cm2/s at 1 keV.
    • Note this assumes that every AGN is at redshift zero.
  • Every AGN has a column density NH in the range 1020 - 1025 cm-2. (Remember the units of NH in TBabs is 1022 cm-2!)
  • Every AGN observation will be carried out by NuSTAR. Note the energy range for spectral analysis will be 3-78 keV, and you can assume observation exposures of 100 ks with no background present. Response (nustar.rmf) and effective area (point_30arcsecRad_1arcminOA.arf) files can be found in the point source simulation files (distributed by the NuSTAR team, more info here).
  • The parent NH distribution follows this shape:

    To generate N random log NH values from this distribution, you can use the following code:
    import scipy.interpolate
    import numpy as np
    fractions = [0.25, 0.05, 0.15, 0.2, 0.3]
    bins = [20., 21., 22., 23., 24., 25.]
    cdf = np.append(0., np.cumsum(fractions) / np.sum(fractions))
    interp_cdf = scipy.interpolate.interp1d(cdf, bins, bounds_error = False)
    rands = np.random.uniform(size = N)
    random_lognh_values = interp_cdf(rands)
Your challenge:
  • Simulate 30 AGN spectra from the model described in the assumptions above, each with a randomly generated log NH value.
  • Re-fit each spectrum with BXA to generate log NH posteriors for each source.
  • Combine the individual log NH posteriors with PosteriorStacker to estimate the NH distribution with the same bins as in the above plot.
  • Estimate the fractions for each log NH bin from the histogram model results["samples"], e.g.:
    lo_hist = np.quantile(result["samples"], 0.16, axis = 0)
    mid_hist = np.quantile(result["samples"], 0.5, axis = 0)
    hi_hist = np.quantile(result["samples"], 0.84, axis = 0)
  • Can your NH distribution constraints be improved? How could spectra below 3 keV help constrain the NH distribution?
  • To find the version of your bxa install:
    import bxa, pkg_resources

  • Speeding up BXA: For complex fits (e.g., multiple datasets simultaneously, many parameters, high signal-to-noise spectra, ...), the run time can be decreased with e.g.,:
    • Parallelisation with MPI (Message Passing Interface):
      mpiexec -np 4 python3 myscript.py
    • Arguments for solver.run()
      • frac_remain=0.5 sets the cut-off condition for the integration. A low value is optimal for finding peaks but higher values (e.g., 0.5) can be used if the posterior distribution is a relatively simple shape.
      • max_num_improvement_loops=0 sets the number of times run() will try to iteratively assess how many more samples are needed.
      • speed
        • speed="safe"
          Reactive Nested Sampler, recommended.
        • speed="auto"
          Step sampler with adaptive steps.
        • speed=int
          Step sampler with an integer number of steps. Faster, but the number of steps should be tested to ensure the log Z estimation is accurate.
      • To see the other arguments of run(), see here.

    • If you're sure you want to, you can use BXA with chi-squared statistics:

    • You can create custom priors in BXA with create_custom_prior_for() (example taken from example_advanced_prior.py):
      def my_custom_prior(u):
          # prior distributions transform from 0:1 to the parameter range
          # here: a gaussian prior distribution, cut below 1/above 3
          x = scipy.stats.norm(1.9, 0.15).ppf(u)
          if x < 1.:
              x = 1
          if x > 3:
              x = 3
          return x
      transformations = [bxa.create_custom_prior_for(model, parameter, my_custom_prior)]

    • You can use BXA with Sherpa to fit 2D models to images:
      import bxa.sherpa as bxa
      model = gauss2d.g1 + gauss2d.g2 + const2d.bg
      set_par(g1.ampl, val = 20., min = 1., max = 1.e2)
      parameters = [g1.ampl, ...]
      priors = [bxa.create_loguniform_prior_for(g1.ampl), ...]
      priorfunction = bxa.create_prior_function(priors)
      solver = bxa.BXASolver(prior=priorfunction, parameters=parameters, outputfiles_basename = "fit2d_sherpa")
      results = solver.run(resume=True)

      Here is an example fit to the image used in this Sherpa tutorial for 2D fitting (left, middle and right panels show the image, best-fit model and residuals respectively):

  • Using other sampling algorithms: Once the solver object has been created, you can acquire the log-Likelihood, parameter names and priors as follows:
    • PyXspec:
      loglike = solver.log_likelihood
      paramnames = solver.paramnames
      prior = solver.prior_function
    • Sherpa:
      loglike = solver.log_likelihood
      paramnames = solver.paramnames
      prior = solver.prior_transform

    These can then be passed to other sampling algorithms (including nested sampling) easily:
    • MCMC algorithms (e.g., emcee or zeus):
      from autoemcee import ReactiveAffineInvariantSampler
      sampler = ReactiveAffineInvariantSampler(paramnames, loglike, prior)
      result = sampler.run()
    • Laplace approximation, Importance Sampling:
      from snowline import ReactiveImportanceSampler
      sampler = ReactiveImportanceSampler(paramnames, loglike, transform)
      result = sampler.run()
    • Nested sampling:
      from ultranest import ReactiveNestedSampler
      sampler = ReactiveNestedSampler(paramnames, loglike, prior)
      result = sampler.run()
  • Plotting contours:
    There are a number of packages available to produce contour plots in Python. Here are a few examples, with some short code snippets to plot BXA posterior corner plots. Note the code snippets assume the posteriors are a pandas dataframe. To do this, you can type:
    import pandas as pd
    outputfiles_basename = "example_basename"
    posterior = pd.read_csv("./%s/chains/equal_weighted_post.txt" %(example_basename), delim_whitespace=True)

  • Note the confidence regions in 2D histograms are not the same as in 1D histograms. See this webpage from corner for more information.