### A Simple Bayesian analysis of X-ray spectrum in Sherpa¶

First import Sherpa

In [1]:
from sherpa.astro.ui import *

WARNING: imaging routines will not be available,
failed to import sherpa.image.ds9_backend due to
'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH'


We are only working with spectra, so imaging with ds9 will not be shown. We will use "matplotlib" for visualization.

Set Poisson likelihood statistics called "cash"in Sherpa and select an optimization method for the initial fit - simplex method which uses nelder-mead algorithm.

In [2]:
set_stat('cash')


Next load the data, filter and setup a simple absorbed power law model. Note that the calibration files, ARF and RMF, which are needed for the X-ray analysis and a background fileare read and when the file "acis.pi" is loaded into the session.

In [3]:
load_data("acis.pi")
ignore(":0.5,7.0:")
set_model(xsphabs.abs1*xspowerlaw.p1)

read ARF file acis.arf


Run the initial fit using simplex and cash. This returns the best fit model parameters, absorption column, power law photon index and normalization.

In [4]:
fit()

WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
Dataset               = 1
Statistic             = cash
Initial fit statistic = 3.11699e+07
Final fit statistic   = -89796.9 at function evaluation 602
Data points           = 444
Degrees of freedom    = 441
Change in statistic   = 3.12597e+07
abs1.nH        0.0668095
p1.PhoIndex    1.17301
p1.norm        0.000531594


Next get the covariance matrix with covar(). The matrix is later used in the MCMC run.

In [5]:
covar()

WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
Dataset               = 1
Confidence Method     = covariance
Iterative Fit Method  = None
Statistic             = cash
covariance 1-sigma (68.2689%) bounds:
Param            Best-Fit  Lower Bound  Upper Bound
-----            --------  -----------  -----------
abs1.nH         0.0668133  -0.00756632   0.00756632
p1.PhoIndex       1.17303   -0.0239168    0.0239168
p1.norm       0.000531608 -1.34321e-05  1.34321e-05


Plot the model and the data.

In [5]:
plot_fit()

WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with cash


There are a few samplers in Sherpa that can be used to explore the parameter space. The standard "metropolismh" runs MCMC with Metropolis-Metropolis-Hastings criterium for accepting parameters.

In [7]:
list_samplers()

Out[7]:
['metropolismh', 'fullbayes', 'mh', 'pragbayes']

Select MetropolisMH and check the options:

In [8]:
set_sampler("MetropolisMH")
get_sampler()

Out[8]:
{'defaultprior': True,
'inv': False,
'log': False,
'originalscale': True,
'p_M': 0.5,
'priors': (),
'priorshape': False,
'scale': 1,
'sigma_m': False}

Here we list several options for running the MCMC sampler:

defaultprior – indicates that all parameters have the default flat prior

inv – indicates which parameters on the inverse scale.

log – indicates which parameters are on the logarithm scale (natural log).

originalscale – shows that the parameters are on the original scale.

p_M – the proportion of jumps generated by the MCMC jumping rule.

priorshape – indicates which parameters have user-defined prior functions.

scale – A scalar multiple of the output of covar() used in the scale of the t-distribution

We first run the MCMC with the default settings:

In [9]:
stats, accept, params = get_draws(niter=1e4)
# check the parameters at the minimum statistic value based on the sampler:

WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
Using Priors:
abs1.nH: <function flat at 0x10295cc80>
p1.PhoIndex: <function flat at 0x10295cc80>
p1.norm: <function flat at 0x10295cc80>

In [10]:
bestfit = params[::,stats.argmin()].T
# print bestfit values and standard deviations
print bestfit
[params[i].std() for i in [0,1,2]]

[  6.68133131e-02   1.17302824e+00   5.31608044e-04]

Out[10]:
[0.007682098051382813, 0.024372684959065014, 1.3781256833497648e-05]
In [11]:
# visualization of the MCMC run - Cash at each iteration
plot_trace(stats)

In [13]:
plot_trace(params[0])

In [14]:
plot_trace(params[1])

In [15]:
# dependence between parameters
plot_scatter(params[0],params[1])

In [16]:
plot_pdf(params[1])

In [17]:
plot_cdf(params[1])
print get_cdf_plot().median
print get_cdf_plot().upper
print get_cdf_plot().lower

1.17231462648
1.19718166318
1.14854738411


### Effects of Priors¶

Changing the prior - here an example of prior which has a gaussian shape. We use one of the predefined Sherpa models to set the prior.

In [18]:
normgauss1d.g1
g1.pos=2.
g1.fwhm=0.5
set_prior(p1.PhoIndex,g1)

In [19]:
set_sampler_opt('defaultprior',False)
set_sampler_opt('priorshape', [False, True, False])
set_sampler_opt('originalscale', [True, True, True])

In [20]:
#Run sampler
stats, accept, params = get_draws(niter=1e4)

WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
Using Priors:
abs1.nH: <function flat at 0x10295cc80>
p1.PhoIndex: normgauss1d.g1
Param        Type          Value          Min          Max      Units
-----        ----          -----          ---          ---      -----
g1.fwhm      thawed          0.5  1.17549e-38  3.40282e+38
g1.pos       thawed            2 -3.40282e+38  3.40282e+38
g1.ampl      thawed            1 -3.40282e+38  3.40282e+38
p1.norm: <function flat at 0x10295cc80>

In [31]:
plot_trace(params[1])
plot_scatter(params[0],params[1])

In [29]:
plot_scatter(params[0][1000:2000],params[1][1000:2000])

In [27]:
plot_cdf(params[1])

In [ ]: