MCMC

class sherpa.sim.MCMC[source] [edit on github]

Bases: NoNewAttributesAfterInit

High-level UI to pyBLoCXS that joins the loop in ‘Walk’ with the jumping rule in ‘Sampler’. Implements a user interface for configuration. This class implements a calc_stat() function using the Sherpa interface to ‘Fit’.

Attributes Summary

`sampler`
`walker`

Methods Summary

`get_draws`(fit, sigma[, niter, cache, rng])	Run the pyBLoCXS MCMC algorithm.
`get_prior`(par)	Return the prior function for a parameter.
`get_sampler`()	Return the current MCMC sampler options.
`get_sampler_name`()	Return the name of the current MCMC sampler.
`get_sampler_opt`(opt)	Return an option of the current MCMC sampler.
`list_priors`()	Return the priors set for model parameters, if any.
`list_samplers`()	List the samplers available for MCMC analysis with `get_draws`.
`set_prior`(par, prior)	Set the prior function to use with a parameter.
`set_sampler`(sampler)	Set the MCMC sampler.
`set_sampler_opt`(opt, value)	Set an option for the current MCMC sampler.

Attributes Documentation

sampler

walker

Methods Documentation

get_draws(fit, sigma, niter=1000, cache=True, rng=None)[source] [edit on github]

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. Unlike many MCMC algorithms, it is designed to explore the parameter space at the suspected statistic minimum (i.e. after using fit). The return values include the statistic value, parameter values, and a flag indicating whether the row represents a jump from the current location or not.

Added in version 4.16.0: The rng parameter was added.

Parameters:

fit – The Sherpa fit object to use.
sigma – The covariance matrix, defined at the best-fit parameter values.
niter (int, optional) – The number of draws to use. The default is 1000.
rng (numpy.random.Generator, numpy.random.RandomState, or None, optional) – Determines how random numbers are created. If set to None then the routines from numpy.random are used, and so can be controlled by calling numpy.random.seed.

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.

Return type:

stats, accept, params

get_prior(par)[source] [edit on github]

Return the prior function for a parameter.

Parameters:: par (sherpa.models.parameter.Parameter) – A parameter of a model instance.
Returns:: The function or parameter instance set by a previous call to set_prior.
Return type:: prior
Raises:: ValueError – If a prior has not been set for the parameter.

See also

set_prior: Set the prior function to use with a parameter.

Examples

>>> pfunc = get_prior(bgnd.c0)

get_sampler()[source] [edit on github]

Return the current MCMC sampler options.

Returns:: options – A copy of the options for the chosen sampler. Use set_sampler_opt to change these values. The fields depend on the current sampler.
Return type:: dict

See also

get_sampler_name: Return the name of the current MCMC sampler.
get_sampler_opt: Return an option of the current MCMC sampler.
set_sampler: Set the MCMC sampler.
set_sampler_opt: Set an option for the current MCMC sampler.

get_sampler_name()[source] [edit on github]

Return the name of the current MCMC sampler.

Returns:: name
Return type:: str

See also

get_sampler: Return the current MCMC sampler options.
set_sampler: Set the MCMC sampler.

Examples

>>> get_sampler_name()
'MetropolisMH'

get_sampler_opt(opt)[source] [edit on github]

Return an option of the current MCMC sampler.

Returns:: opt – The name of the option. The fields depend on the current sampler.
Return type:: str

See also

get_sampler: Return the current MCMC sampler options.
set_sampler_opt: Set an option for the current MCMC sampler.

Examples

>>> get_sampler_opt('log')
False

list_priors()[source] [edit on github]

Return the priors set for model parameters, if any.

Returns:: priors – The dictionary of mappings between parameters (keys) and prior functions (values) created by set_prior.
Return type:: dict

See also

get_prior: Return the prior function for a parameter.
set_prior: Set the prior function to use with a parameter.

list_samplers()[source] [edit on github]

List the samplers available for MCMC analysis with get_draws.

Returns:: samplers – A list of the names (in lower case) that can be used with set_sampler.
Return type:: list of str

See also

get_sampler_name: Return the name of the current MCMC sampler.
set_sampler: Set the MCMC sampler.

Notes

The available samplers depends on what modules have been loaded, and may be more extensive than the example output below.

Examples

>>> list_samplers()
['metropolismh', 'mh']

set_prior(par, prior)[source] [edit on github]

Set the prior function to use with a parameter.

The default prior used by get_draws for each parameter is flat, varying between the hard minimum and maximum values of the parameter (as given by the hard_min and hard_max attributes of the parameter object). The set_prior function is used to change the form of the prior for a parameter.

Parameters:

par (sherpa.models.parameter.Parameter instance) – A parameter of a model instance.
prior (function or sherpa.models.model.Model instance) – The function to use for a prior. It must accept a single argument and return a value of the same size as the input.

See also

get_draws: Run the MCMC algorithm.
get_prior: Set the prior function to use with a parameter.
set_sampler: Set the MCMC sampler.

Examples

Set the prior for the kT parameter of the therm instance to be a gaussian, centered on 1.7 keV and with a FWHM of 0.35 keV:

>>> create_model_component('xsapec', 'therm')
>>> create_model_component('gauss1d', 'p_temp')
>>> p_temp.pos = 1.7
>>> p_temp.fwhm = 0.35
>>> set_prior(therm.kT, p_temp)

Create a function (lognorm) and use it as the prior of the nH parameter of the abs1 instance:

>>> def lognorm(x):
...     nh = 20
...     sigma = 0.5  # use a sigma of 0.5
...     # nH is in units of 10^-22 so convert
...     dx = np.log10(x) + 22 - nh
...     norm = sigma / np.sqrt(2 * np.pi)
...     return norm * np.exp(-0.5 * dx * dx / (sigma * sigma))
...
>>> create_model_component('xsphabs', 'abs1')
>>> set_prior(abs1.nH, lognorm)

set_sampler(sampler)[source] [edit on github]

Set the MCMC sampler.

The sampler determines the type of jumping rule to be used by get_draws.

Parameters:: sampler (str or sherpa.sim.Sampler instance) – When a string, the name of the sampler to use (case insensitive). The supported options are given by the list_samplers function.

See also

get_draws: Run the MCMC algorithm.
list_samplers: List the MCMC samplers.
set_sampler: Set the MCMC sampler.
set_sampler_opt: Set an option for the current MCMC sampler.

Notes

The jumping rules are:

MH: The Metropolis-Hastings rule, which always jumps from the best-fit location, even if the previous iteration had moved away from it.
MetropolisMH: This is the Metropolis with Metropolis-Hastings algorithm, that jumps from the best-fit with probability p_M, otherwise it jumps from the last accepted jump. The value of p_M can be changed using set_sampler_opt.
PragBayes: This is used when the 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 (the nsubiters option) MCMC sub-iteration steps are carried out, choosing between Metropolis and Metropolis-Hastings types of samplers with probability p_M. Only the last of these sub-iterations are kept in the chain. The nsubiters and p_M values can be changed using set_sampler_opt.
FullBayes: Another sampler for use when including uncertainties due to the effective area.

Examples

>>> set_sampler('metropolismh')

set_sampler_opt(opt, value)[source] [edit on github]

Set an option for the current MCMC sampler.

Parameters:

opt (str) – The option to change. Use get_sampler to view the available options for the current sampler.
value – The value for the option.

See also

get_sampler: Return the current MCMC sampler options.
set_prior: Set the prior function to use with a parameter.
set_sampler: Set the MCMC sampler.

Notes

The options depend on the sampler. The options include:

defaultprior: Set to False when the default prior (flat, between the parameter’s soft limits) should not be used. Use set_prior to set the form of the prior for each parameter.
inv: A bool, or array of bools, to indicate which parameter is on the inverse scale.
log: A bool, or array of bools, to indicate which parameter is on the logarithm (natural log) scale.
original: A bool, or array of bools, to indicate which parameter is on the original scale.
p_M: The proportion of jumps generatd by the Metropolis jumping rule.
priorshape: An array of bools indicating which parameters have a user-defined prior functions set with set_prior.
scale: Multiply the output of covar by this factor and use the result as the scale of the t-distribution.

Examples

>>> set_sampler_opt('scale', 3)