get_draws
- sherpa.astro.ui.get_draws(id: IdType | None = None, otherids: Sequence[IdType] = (), niter=1000, covar_matrix=None)
Run the pyBLoCXS MCMC algorithm.
The function runs a Markov Chain Monte Carlo (MCMC) algorithm designed to carry out Bayesian Low-Count X-ray Spectral (BLoCXS) analysis. It explores the model parameter space at the suspected statistic minimum (i.e. after using
fit
). The return values include the statistic value, parameter values, and an acceptance flag indicating whether the row represents a jump from the current location or not. For more information see thesherpa.sim
module and the reference given below.- Parameters:
id (int, str, or None, optional) – The data set that provides the data. If not given then all data sets with an associated model are used simultaneously.
otherids (sequence of int or str, optional) – Other data sets to use in the calculation.
niter (int, optional) – The number of draws to use. The default is
1000
.covar_matrix (2D array, optional) – The covariance matrix to use. If
None
then the result fromget_covar_results().extra_output
is used.
- Returns:
The results of the MCMC chain. The stats and accept arrays contain
niter+1
elements, with the first row being the starting values. The params array has(nparams, niter+1)
elements, where nparams is the number of free parameters in the model expression, and the first column contains the values that the chain starts at. The accept array contains boolean values, indicating whether the jump, or step, was accepted (True
), so the parameter values and statistic change, or it wasn’t, in which case there is no change to the previous row. Thesherpa.utils.get_error_estimates
routine can be used to calculate the credible one-sigma interval from the params array.- Return type:
stats, accept, params
See also
covar
,fit
,get_sampler
,plot_cdf
,plot_pdf
,plot_scatter
,plot_trace
,set_prior
,set_rng
,set_sampler
Notes
The chain is run using fit information associated with the specified data set, or sets, the currently set sampler (
set_sampler
) and parameter priors (set_prior
), for a specified number of iterations. The model should have been fit to find the best-fit parameters, andcovar
called, before runningget_draws
. The results fromget_draws
is used to estimate the parameter distributions.The
set_rng
routine is used to control how the random numbers are generated.References
Examples
Fit a source and then run a chain to investigate the parameter distributions. The distribution of the stats values created by the chain is then displayed, using
plot_trace
, and the parameter distributions for the first two thawed parameters are displayed. The first one as a cumulative distribution usingplot_cdf
and the second one as a probability distribution usingplot_pdf
. Finally the acceptance fraction (number of draws where the chain moved) is displayed. Note that in a full analysis session a burn-in period would normally be removed from the chain before using the results.>>> fit() >>> covar() >>> stats, accept, params = get_draws(1, niter=1e4) >>> plot_trace(stats, name='stat') >>> names = [p.fullname for p in get_source().pars if not p.frozen] >>> plot_cdf(params[0,:], name=names[0], xlabel=names[0]) >>> plot_pdf(params[1,:], name=names[1], xlabel=names[1]) >>> accept[:-1].sum() * 1.0 / len(accept - 1) 0.4287
The following runs the chain on multiple data sets, with identifiers ‘core’, ‘jet1’, and ‘jet2’:
>>> stats, accept, params = get_draws('core', ['jet1', 'jet2'], niter=1e4)