Last modified: November 2016

Jump to: Description · Examples · Bugs · See Also

AHELP for CIAO 4.12 Sherpa v1


Context: methods


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


get_draws(id=None, otherids=(), niter=1000)


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 - i.e. after a standard Sherpa fit.

The Sherpa get_draws() function runs a pyBLoCXS chain using fit information associated with the specified data set(s), the currently set sampler 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. It can only be used with a likelihood-based statistic (e.g. not any of the chi-square variants supported by Sherpa).

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 command. The available samplers are returned by the list_samplers command.

Jumping Rules

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

By default, pyBLoCXS uses a flat prior defined between the hardcoded parameter minima and maxima. The set_prior() command is used to associate a function or model ('prior') with a thawed fit parameter ('par'). The list of currently set prior-parameter pairs is returned by list_priors(); and the prior function associated with a particular Sherpa model parameter may be accessed with get_prior().

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.

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.

Refer to the pyBLoCXS documentation for additional information about the algorithm.


Example 1

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_source(xsphabs.abs1 * powlaw1d.p1)
sherpa> set_stat("cash")
sherpa> fit()
sherpa> covar()

sherpa> print(get_sampler_name())

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]

Example 2

By default, get_draws() runs on only one data set; for multiple data sets, "otherids" must be supplied to the function.

Here, the spectra from three observations are fit simultaneously, prior to running get_draws using the default flat prior for all parameters:

sherpa> load_pha(1, "obsid1.pi")
sherpa> load_pha(2, "obsid2,pi")
sherpa> load_pha(3, "obsid3.pi")
sherpa> load_pha(4, "obsid4.pi")
sherpa> set_full_model(1, rsp_1(xsphabs.abs1*xsapec.kt1))
sherpa> set_full_model(2, rsp_2(abs1*kt1))
sherpa> set_full_model(3, rsp_3(abs1*kt1))
sherpa> set_full_model(4, rsp_4(abs1*kt1))
sherpa> abs1.nh = 0.039
sherpa> freeze(abs1)
sherpa> thaw(kt1.Abundanc)
sherpa> fit()
sherpa> covar()

sherpa> stats, accept, params = get_draws(1, [2,3,4], niter=1e4)
Using Priors:
kt1.kT:  <function flat at 0x128b8758>
kt1.Abundanc:  <function flat at 0x128b8758>
kt1.norm:  <function flat at 0x128b8758>

sherpa> print(params)
[[  7.05838558e-01   7.47127420e-01   6.74755992e-01   6.48366180e-01
    6.48366180e-01   6.48366180e-01   6.48366180e-01   6.48366180e-01
    6.48366180e-01   6.48366180e-01   7.67033002e-01   7.67033002e-01
    7.33908084e-01   7.33908084e-01   7.66567682e-01   7.66567682e-01


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

See Also

get_conf, get_covar, get_int_proj, get_int_unc, get_proj, get_reg_proj, get_reg_unc
get_chart_spectrum, get_marx_spectrum
get_areascal, get_arf, get_arf_plot, get_axes, get_backscal, get_bkg, get_bkg_plot, get_bkg_scale, get_coord, get_counts, get_data, get_data_plot, get_dep, get_dims, get_error, get_exposure, get_grouping, get_indep, get_quality, get_rmf, get_specresp, get_staterror, get_syserror
calc_stat_info, get_fit, get_stat_info
get_default_id, list_stats
get_iter_method_name, get_iter_method_opt, get_method
get_model, get_model_component, get_model_component_image, get_model_component_plot, get_model_plot, get_num_par, get_order_plot, get_par, get_pileup_model, get_response, get_source, get_source_component_image, get_source_component_plot, image_source
get_kernel, get_psf
get_chisqr_plot, get_delchi_plot, get_prior, get_sampler, get_stat
get_analysis, get_rate
get_ratio, get_resid, image_getregion