get_draws

Context: methods

Synopsis

Run the pyBLoCXS MCMC-based algorithm using the current sampler.

Syntax

get_draws([id, otherids, niter])

pyBLoCXS is a sophisticated Markov chain Monte Carlo (MCMC) based algorithm designed to carry out Bayesian Low-Count X-ray Spectral (BLoCXS) analysis in the Sherpa environment. The algorithm explores parameter space at a suspected minimum using a predefined Sherpa model to high-energy X-ray spectral data.

The Sherpa get_draws() function runs a pyBLoCXS chain using fit information associated with the specified data set(s), the currently set sampler (either Metropolis-Hastings , a mix of Metropolis with Metropolis-Hastings, or PragBayes) and parameter priors, for a specified number of iterations. It returns an array of statistic values, an array of acceptance Booleans, and a 2-D array of associated parameter values.

The sampler, or jumping rule, to be used by pyBLoCXS may be set with set_sampler(), and its configuration options modifed using the set_sampler_opt() fuction. The available samplers - "MH", "MetropolisMH", and "PragBayes" - are returned by the list_samplers() command ("ahelp list_samplers"). "MH" is Metropolis-Hastings, which always jumps from the best-fit, and "MetropolisMH" is Metropolis with Metropolis-Hastings that jumps from the best-fit with probability 'p_M', else it jumps from the last accepted jump. "PragBayes" is used when effective area calibration uncertainty is to be included in the calculation. (At each nominal MCMC iteration, a new calibration product is generated, and a series of N (option in set_sampler_opt) MCMC sub-iteration steps are carried out, choosing between Metropolis and Metropolis-Hastings types of samplers with probability p_M (option in set_sampler_opt). Only the last of these sub-iterations are kept in the chain.)

Note that before running get_draws(), a Sherpa fit must be run and the covariance matrix should be calculated at the resultant fit minimum; see the Examples section below for a demonstration.

id - ID of the data set which has been fit in preparation for the run of the pyBLoCXS algorithm; default=1
otherids - list of additional data set IDs with fit information to use in the run of the pyBLoCXS algorithm; default=None
niter - number of iterations to run; default = 1e3

For a detailed description of the algorithm, see:

van Dyk, D.A., Connors, A., Kashyap, V.L., & Siemiginowska, A. 2001, Ap.J., 548, 224, "Analysis of Energy Spectra with Low Photon Counts via Bayesian Posterior Simulation"

and

van Dyk, D.A., Connors, A., Kashyap, V.L., & Siemiginowska, A. 2001, Ap.J., 548, 224, "Analysis of Energy Spectra with Low Photon Counts via Bayesian Posterior Simulation"

pyBLoCXS documentation is also available here.

Example

Fit an absorbed power-law model to source data set "1", and then run the covariance function to calculate confidence intervals for the thawed parameters. Run the pyBLoCXS algorithm using the currently set Metropolis with Metropolis-Hastings jumping rule and fit information from data set 1, 1e4 times. Store the resulting arrays of statistic values, acceptance flags, and parameter values to variables.

sherpa> load_pha("pha.fits")
sherpa> set_model(xsphabs.abs1*powlaw1d.p1)
sherpa> set_stat("cash")
sherpa> fit()
sherpa> covar()

sherpa> print get_sampler_name()
        MetropolisMH

sherpa> stats, accept, params = get_draws(1, niter=1e4)

Using Priors:
abs1.nH: <function flat at 0x1b948e60>
p1.PhoIndex: <function flat at 0x1b948e60>
p1.norm: <function flat at 0x1b948e60>

sherpa> print stats
[-2449.71534006 -2449.71534006 -2449.71534006 -2447.74139152 -2447.74139152 -2447.74139152 -2441.5488066  -2441.5488066  -2441.5488066 -2446.84464469 ...  -2447.14055844]

Bugs

See the bugs pages on the Sherpa website for an up-to-date listing of known bugs.