• 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

• 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

• 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

• Set the x-axis units (e.g., instrument channels or energy):
  • Xspec:
    setplot energy 
  • PyXspec:
    Plot.xAxis = "keV"# or "channel" 
  • Sherpa:
    set_analysis("ener")# or "chan" 


• 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)

• Load spectra:
  • Xspec:
    data 1:1 sim_chan_bmin1.pi 
    
  • PyXspec:
    s = Spectrum("sim_chan_bmin1.pi") 
    Or:
    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.-**") 
    Or:
    AllData.ignore("1:0.-0.5 8.-**") 
  • Sherpa:
    ignore_id(1, "0.:0.5,8.:") 

• 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.")  

• 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)
  Plot("ldata")
  

   
• Sherpa (note you can try "counts" in place of "rate"):
  

  group_snr(1, 2.)
  subtract(1)
  set_analysis(1, "ener", "rate")
  plot_data(xlog=True, ylog=True)
  

   

Interactive PLT (Xspec only)

• 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  

• 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  

• Xspec:
  

  setplot rebin 2 1000 -1
  setplot command wdata ldata_plot.dat
  plot ldata
  exit
  

   
  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.)
  subtract(1)
  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})
  

   

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)
ax.set_xscale("log")
ax.set_yscale("log")
plt.show()
fig.savefig("folded_data.png")

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

• PyXspec:
  

  from xspec import *
  import bxa.xspec as bxa
  

   
• Sherpa:
  

  import bxa.sherpa as bxa
  

   

• PyXspec:
  

  bxa.create_uniform_prior_for(model, param)
  bxa.create_loguniform_prior_for(model, param)
  bxa.create_gaussian_prior_for(model, param, mean, std)
  

   
• Sherpa:
  

  bxa.create_uniform_prior_for(param)
  bxa.create_loguniform_prior_for(param)
  bxa.create_gaussian_prior_for(param, mean, std)
  

   

• 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):
  • PyXspec:
    

    prior1 = bxa.create_uniform_prior_for(mymod, mymod.powerlaw.PhoIndex)
    prior2 = bxa.create_loguniform_prior_for(mymod, mymod.powerlaw.norm)
    

     
  • Sherpa:
    

    prior1 = bxa.create_uniform_prior_for(mymod.PhoIndex)
    prior2 = bxa.create_loguniform_prior_for(mymod.norm)
    

     
  
  
• 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],
                           outputfiles_basename="powerlaw_sherpa")
    

     
  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
       set_parameters(transformations=solver.transformations,values=values)
       

        
     • 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:
       Plot("ldata")  
     • Sherpa:
       plot_fit(xlog=True, ylog=True)  

• Produce 10 realisations of the model priors. Are the model ranges covered by the priors physically-acceptable?

• First fit the model using Levenberg-Marquadt included in Xspec & Sherpa:
  • PyXspec:
    Fit.perform()  
  • Sherpa:
    fit()  


• 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)
  df.head()

   


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

• Note to use TBabs, one should set the abundances to those from Wilms, Allen & McCray (2000):
  • PyXspec:
    Xset.abund = "wilm"  
  • Sherpa:
    set_xsabund("wilm")  
  
  
• 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.
     
     • PyXspec:
       TBabs * powerlaw  
     • Sherpa:
       xstbabs * xspowerlaw  
  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):
     
     • PyXspec (note to fix a parameter, you can set delta = -1):
       TBabs * powerlaw + gaussian  
     • Sherpa:
       xstbabs * xspowerlaw + xsgaussian  
     
     
• 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?

• 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?

• 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 ...")
    plt.figure(figsize=(7,7))
    with bxa.XSilence():
        solver.set_best_fit()
        bxa.qq.qq(prefix=outputfiles_basename, markers=5, annotate=True)
    print("saving plot...")
    plt.savefig(outputfiles_basename + "qq_model_deviations.pdf", bbox_inches="tight")
    plt.close()

     
  • 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?

• 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).

• First, include a cross-calibration constant in your model between the Chandra and NuSTAR datasets:
  
  • PyXspec:
    

    s = AllData("1:1 sim_nuA_bmin1.pi 2:2 sim_chan_bmin1.pi")
    AllData.ignore("1:0.-3. 78.-** 2:0.-0.5 8.-**")
    mymod = Model("constant * powerlaw")
    AllModels(1)(1).values = (1., -0.01)
    AllModels(2)(1).values = (1., 0.01, 0.01, 0.01, 10., 10.)

     
  • Sherpa:
    

    load_pha(1, "sim_nuA_bmin1.pi")
    ignore_id(1, "0.:3.,78.:")
    load_pha(2, "sim_chan_bmin1.pi")
    ignore_id(2, "0.:0.5,8.:")
    mymod = xspowerlaw.mypow
    mymod1 = xsconstant("xcal1") * mymod
    set_par(xcal1.factor, val=1, frozen=True)
    mymod2 = xsconstant("xcal2") * mymod
    set_par(xcal2.factor, val=1, min = 1.e-2, max = 1.e2)
    set_source(1, mymod1)
    set_source(2, mymod2)

     


• 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]
    priors.append(bxa.create_loguniform_prior_for(xcal2))
    parameters.append(xcal2)
    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?

• 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?

• Combine individual posteriors to constrain a population
• Calibrate model selection thresholds with simulations
• Perform model verification to test the Goodness-of-Fit

• 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)
  

   

• 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.xlabel("PhoIndex")
  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)
  plt.show()
  

   
  
• 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?

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.

• 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. logZ₁/Z₂ < 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?

• 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?

• 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?

• 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):
  PyXspec:
  

  stat = Fit.statistic
  dof = Fit.dof
  

   
  Sherpa:
  

  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?

• Every AGN follows the same baseline model:
  model = TBabs[Gal] * TBabs[AGN] * powerlaw
  Where:
  • TBabs[Gal] is the absorbing column density arising from the Milky Way, assumed to be N_H = 10²⁰ 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 N_H 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/cm²/s at 1 keV.
  Note this assumes that every AGN is at redshift zero.
• Every AGN has a column density N_H in the range 10²⁰ - 10²⁵ cm^-2. (Remember the units of N_H in TBabs is 10²² 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 N_H distribution follows this shape:
  
  
  
  
  To generate N random log N_H values from this distribution, you can use the following code:

  import scipy.interpolate
  import numpy as np
  np.random.seed(42)
  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)
  

   

• Simulate 30 AGN spectra from the model described in the assumptions above, each with a randomly generated log N_H value.
• Re-fit each spectrum with BXA to generate log N_H posteriors for each source.
• Combine the individual log N_H posteriors with PosteriorStacker to estimate the N_H distribution with the same bins as in the above plot.
• Estimate the fractions for each log N_H 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 N_H distribution constraints be improved? How could spectra below 3 keV help constrain the N_H distribution?

• To find the version of your bxa install:
  

  import bxa, pkg_resources
  print(pkg_resources.get_distribution("bxa").version)
  

   


• 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:
    bxa.BXASolver.allowed_stats.append("chi2") 
  
  
  • You can create custom priors in BXA with create_custom_prior_for() (example taken from example_advanced_prior.py):
    PyXspec:
    

    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
    load_image("image.fits")
    image_data()
    set_coord("physical")
    set_stat("cash")
    model = gauss2d.g1 + gauss2d.g2 + const2d.bg
    set_source(model)
    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)
    image_data()
    image_model(newframe=True)
    image_resid(newframe=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)
  

   
  • corner:
    corner.readthedocs.io/en/latest/pages/quickstart.html
    

    import corner
    fig = corner.corner(posterior.values, labels=posterior.columns, quantiles=[0.16, 0.5, 0.84], show_titles=True)
    fig.savefig("/path/to/your/corner.png")
    

     
  • ChainConsumer:
    samreay.github.io/ChainConsumer/examples/plot_introduction.html
    

    from chainconsumer import ChainConsumer
    c=ChainConsumer()
    c.add_chain(posterior.values, parameters = [r"%s" %(p) for p in posterior.columns])
    fig=c.plotter.plot(figsize=(10, 10))
    fig.savefig("/path/to/your/corner.png")
    

     
  • GetDist:
    getdist.readthedocs.io/en/latest/plot_gallery.html
    

    from getdist import plots, MCSamples
    samples1 = MCSamples(samples=posterior.values, names=posterior.columns)
    g = plots.get_subplot_plotter()
    g.triangle_plot(samples1, filled=True)
    g.export("/path/to/your/corner.png")
    

     

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