Samplers

Typically, samplers are instantiated this way:

sample:
  # Specify the sampler class
  class: cosmofit.samplers.EmceeSampler
  init:
    # Any sampler option
    chains: 4  # number of chains or paths to existing chains to resume from
    seed: None  # random seed
    ref_scale: 1.  # rescale all parameter reference distribution (from which they are initially sampled from) by this factor
    max_tries: 1000  # maximum number of calls to get finite posterior
  run:
    check_every: 200  # save samples every 200 iterations
    min_iterations: 100  # minimum number of iterations to run (useful if convergence criteria below are satisfied by chance at the beginning of the run)
    max_iterations: sys.maxsize  # maximum number of iterations to run
    check:
      # Fraction of samples to remove for convergence tests
      burnin: 0.5
      # For how many "checks" the criteria must be fulfilled for the inference to stop
      stable_over: 2
      # Gelman-Rubin criterion (on eigenvalues of the parameter covariance matrix) < 0.03
      max_eigen_gr: 0.03
      # Gelman-Rubin criterion (on parameter covariance matrix diagonal) < 0.03
      max_diag_gr: 0.03
      # Gelman-Rubin criterion on variance of nsigmas_cl_diag_gr = 1 interval limits < 0.03
      max_cl_diag_gr: 0.03
      nsigmas_cl_diag_gr: 1
      # Minimal number of iterations over integrated auto-correlation time (~ # of independent samples)
      min_iterations_over_iact: 1000
      reliable_iterations_over_iact: 50  # after how many samples is auto-correlation time estimation reliable
      # All max_* have a min_* counterpart, and vice-versa
  save: 'chain_*.npy'  # where to save chains

Note however that for nested samplers (dynesty, polychord, “check” options are not the same). Sampling can be interrupted anytime, and resumed by providing the path to the saved chains in “chains” argument of “init”.

One will typically run sampling nchains * nprocs_per_chain + 1 processes, with nchains >= 1 the number of chains and nprocs_per_chain >= 1 the number of processes per chain — plus 1 root process to distribute the work.

In the following we present samplers’ default options (in addition to those presented above).

dynesty

class: StaticDynestySampler

info:
  version: 0.0.1
  date: 05/04/2022
  maintainer: Arnaud de Mattia
  description: Wrapper for dynesty (static) sampler
  url: https://github.com/joshspeagle/dynesty
  doi: [10.1093/mnras/staa278, 10.5281/zenodo.7215695]
  long_description: 'Nested Sampling. Proper priors only are supported.
                     Using less “informative” priors will increase the expected number of nested sampling iterations.
                     Static nested sampling is designed to estimate the evidence. For posterior estimation, rather use dynamic nested sampling.'

init:
  # Number of "live" points. Larger numbers result in a more finely sampled posterior (more accurate evidence), but also a larger
  # number of iterations required to converge. Default is `500`.
  nlive: 500
  # Method used to approximately bound the prior using the current
  # set of live points. Conditions the sampling methods used to
  # propose new live points. Choices are no bound ('none'), a single
  # bounding ellipsoid ('single'), multiple bounding ellipsoids
  # ('multi'), balls centered on each live point ('balls'), and
  # cubes centered on each live point ('cubes'). Default is 'multi'.
  bound: 'multi'
  # Method used to sample uniformly within the likelihood constraint,
  # conditioned on the provided bounds. Unique methods available are:
  # uniform sampling within the bounds('unif'),
  # random walks with fixed proposals ('rwalk'),
  # multivariate slice sampling along preferred orientations ('slice'),
  # "random" slice sampling along all orientations ('rslice'),
  # "Hamiltonian" slices along random trajectories ('hslice'), and
  # any callable function which follows the pattern of the sample methods
  # defined in dynesty.sampling.
  # 'auto' selects the sampling method based on the dimensionality of the problem (ndim).
  # When ndim < 10, this defaults to 'unif'.
  # When 10 <= ndim <= 20, this defaults to 'rwalk'.
  # When 'ndim > 20', this defaults to 'hslice' if a 'gradient' is provided and 'rslice' otherwise.
  # 'slice' is provided as alternatives for 'rslice'.
  # Default is 'auto'.
  sample: 'auto'
  # If an integer is passed, only update the proposal distribution
  # every update_interval-th likelihood call. If a float is passed,
  # update the proposal after every round(update_interval * nlive)-th likelihood call.
  # Larger update intervals larger can be more efficient when the likelihood function is quick to evaluate.
  # Default behavior is to target a roughly constant change in prior volume, with
  # 1.5 for 'unif', 0.15 * walks for 'rwalk', 0.9 * ndim * slices for 'slice', 2.0 * slices for 'rslice',
  # and 25.0 * slices for 'hslice'.
  update_interval: None

run:
  check:
    # Iteration will stop when the estimated contribution of the
    # remaining prior volume to the total evidence falls below
    # this threshold. Explicitly, the stopping criterion is
    # ln(z + z_est) - ln(z) < dlogz, where z is the current
    # evidence from all saved samples and z_est is the estimated
    # contribution from the remaining volume.
    # The default is 1e-3 * (nlive - 1) + 0.01.
    #dlogz: 0.01

    # Minimum number of effective posterior samples. If the estimated
    # "effective sample size" (ESS) exceeds this number,
    # sampling will terminate. Default is inf.
    #n_effective:

---

class: DynamicDynestySampler

info:
  version: 0.0.1
  date: 05/04/2022
  maintainer: Arnaud de Mattia
  description: Wrapper for dynesty (static) sampler
  url: https://github.com/joshspeagle/dynesty
  doi: [10.1093/mnras/staa278, 10.5281/zenodo.7215695]
  long_description: 'Nested Sampling. Proper priors only are supported.
                     Using less “informative” priors will increase the expected number of nested sampling iterations.
                     Dynamic nested sampling allocate live points dynamically to sample preferentially in the posterior mass.'

init:
  # Number of "live" points. Larger numbers result in a more finely sampled posterior (more accurate evidence), but also a larger
  # number of iterations required to converge. Default is `500`.
  nlive: 500
  # Method used to approximately bound the prior using the current
  # set of live points. Conditions the sampling methods used to
  # propose new live points. Choices are no bound ('none'), a single
  # bounding ellipsoid ('single'), multiple bounding ellipsoids
  # ('multi'), balls centered on each live point ('balls'), and
  # cubes centered on each live point ('cubes'). Default is 'multi'.
  bound: 'multi'
  # Method used to sample uniformly within the likelihood constraint,
  # conditioned on the provided bounds. Unique methods available are:
  # uniform sampling within the bounds('unif'),
  # random walks with fixed proposals ('rwalk'),
  # multivariate slice sampling along preferred orientations ('slice'),
  # "random" slice sampling along all orientations ('rslice'),
  # "Hamiltonian" slices along random trajectories ('hslice'), and
  # any callable function which follows the pattern of the sample methods
  # defined in dynesty.sampling.
  # 'auto' selects the sampling method based on the dimensionality of the problem (ndim).
  # When ndim < 10, this defaults to 'unif'.
  # When 10 <= ndim <= 20, this defaults to 'rwalk'.
  # When 'ndim > 20', this defaults to 'hslice' if a 'gradient' is provided and 'rslice' otherwise.
  # 'slice' is provided as alternatives for 'rslice'.
  # Default is 'auto'.
  sample: 'auto'
  # If an integer is passed, only update the proposal distribution
  # every update_interval-th likelihood call. If a float is passed,
  # update the proposal after every round(update_interval * nlive)-th likelihood call.
  # Larger update intervals larger can be more efficient when the likelihood function is quick to evaluate.
  # Default behavior is to target a roughly constant change in prior volume, with
  # 1.5 for 'unif', 0.15 * walks for 'rwalk', 0.9 * ndim * slices for 'slice', 2.0 * slices for 'rslice',
  # and 25.0 * slices for 'hslice'.
  update_interval: None

run:
  check:
    # The baseline run will stop when the estimated contribution of the
    # remaining prior volume to the total evidence falls below
    # this threshold. Explicitly, the stopping criterion is
    # ln(z + z_est) - ln(z) < dlogz, where z is the current
    # evidence from all saved samples and z_est is the estimated
    # contribution from the remaining volume.
    # The default is 0.01
    #dlogz_init: 0.01

    # Minimum number of effective posterior samples. If the estimated
    # "effective sample size" (ESS) exceeds this number,
    # sampling will terminate. Default is max(10000, ndim^2).
    #n_effective:

polychord

class: PolychordSampler

info:
  version: 0.0.1
  date: 05/04/2022
  author: Arnaud de Mattia
  maintainer: Arnaud de Mattia
  description: Wrapper for polychord sampler
  url: https://github.com/PolyChord/PolyChordLite
  doi: [10.1093/mnrasl/slv047, 10.1093/mnras/stv1911]
  long_description: 'Nested sampling for cosmology'

init:
  # Parameter blocks are groups of parameters which are updated alltogether with a frequency proportional to oversample_factor
  # Typically, parameter blocks are chosen such that parameters in a given block require the same evaluation time of the likelihood when updated
  # By default these blocks are defined at runtime, based on (measured) speeds and oversample_power (below), but can be specified there in the format
  # - [[param1, param2], oversample_factor1]
  # - [[param3, param4], oversample_factor2]
  blocks:
  # Oversample factors are speed^oversample_power
  oversample_power: 0.4
  # Number of live points. Increasing nlive increases the accuracy of posteriors and evidences,
  # and proportionally increases runtime ~ O(nlive).
  nlive: '25 * ndim'
  # The number of prior samples to draw before starting compression.
  nprior: '10*nlive'
  # The number of failed spawns before stopping nested sampling.
  nfail: '1*nlive'
  # The number of slice slice-sampling steps to generate a new point.
  # Increasing nrepeats increases the reliability of the algorithm.
  # Typically:
  # * for reliable evidences need num_repeats ~ O(5*nDims).
  # * for reliable posteriors need num_repeats ~ O(nDims)
  nrepeats: '2*ndim'
  # Variable number of live points option. This dictionary is a mapping
  # between loglike contours and nlive.
  # You should still set nlive to be a sensible number, as this indicates
  # how often to update the clustering, and to define the default value.
  nlives: {}
  # Whether or not to explore multi-modality on the posterior
  do_clustering: True
  # Increase the number of posterior samples produced. This can be set
  # arbitrarily high, but you won't be able to boost by more than nrepeats
  # Warning: in high dimensions PolyChord produces _a lot_ of posterior
  # samples. You probably don't need to change this
  boost_posterior: 0.
  # Parallelise with synchronous workers, rather than asynchronous ones.
  # This can be set to False if the likelihood speed is known to be
  # approximately constant across the parameter space. Synchronous
  # parallelisation is less effective than asynchronous by a factor ~O(1)
  # for large parallelisation.
  synchronous: True

run:
  # How often to update the files and do clustering
  compression_factor: e'np.exp(-1)'
  check:
    # Nested sampling terminates when the evidence contained in the live points is precision_criterion fraction of the total evidence.
    precision_criterion: 0.001

zeus

class: ZeusSampler

info:
  version: 0.0.1
  date: 05/04/2022
  author: Arnaud de Mattia
  maintainer: Arnaud de Mattia
  description: Wrapper for zeus sampler
  url: https://github.com/minaskar/zeus
  doi: [10.1007/s11222-021-10038-2, 10.1093/mnras/stab2867]
  long_description: Ensemble Slice Sampling method

init:
  # Number of walkers, defaults to 2 * max((int(2.5 * ndim) + 1) // 2, 2)
  # Can be given in dimension units, e.g. 3 * ndim
  nwalkers: None
  # If True (default is False) then no expansions are performed after the tuning phase.
  # This can significantly reduce the number of log likelihood evaluations but works best in target distributions that are apprroximately Gaussian.
  light_mode: False

run:
  # Thin samples by this factor
  thin_by: 1

emcee

class: EmceeSampler

info:
  version: 0.0.1
  date: 05/04/2022
  maintainer: Arnaud de Mattia
  description: Wrapper for emcee sampler
  url: https://github.com/dfm/emcee
  doi: 10.1086/670067
  long_description: 'Python implementation of the affine-invariant ensemble sampler
                     for Markov chain Monte Carlo (MCMC) proposed by Goodman & Weare (2010)'

init:
  # Number of walkers, defaults to 2 * max((int(2.5 * ndim) + 1) // 2, 2)
  # Can be given in dimension units, e.g. 3 * ndim
  nwalkers: None

run:
  # Thin samples by this factor
  thin_by: 1

mcmc

class: MCMCSampler

info:
  version: 0.0.1
  date: 05/04/2022
  maintainer: Arnaud de Mattia
  description: Blocked fast-slow Metropolis sampler
  url: https://github.com/CobayaSampler/cobaya/tree/master/cobaya/samplers/mcmc
  doi: [10.1103/PhysRevD.66.103511, 10.1103/PhysRevD.87.103529, 10.48550/arXiv.math/0502099]
  long_description: 'Antony Lewis CosmoMC sampler, wrapped for cobaya by Jesus Torrado, reimplemented in cosmofit'

init:
  # Parameter blocks are groups of parameters which are updated alltogether with a frequency proportional to oversample_factor
  # Typically, parameter blocks are chosen such that parameters in a given block require the same evaluation time of the likelihood when updated
  # By default these blocks are defined at runtime, based on (measured) speeds and oversample_power (below), but can be specified there in the format
  # - [[param1, param2], oversample_factor1]
  # - [[param3, param4], oversample_factor2]
  blocks:
  # Oversample factors are speed^oversample_power
  oversample_power: 0.4
  # (Initial) proposal covariance (to draw parameter jumps)
  # Can be previous samples ({fn: chain.npy, burnin: 0.5}), or profiles (containing covariance matrix)
  # If variance for a given parameter is not provided, use parameter 'proposal' squared
  covariance:
  # Scale proposal by this value when drawing jumps
  proposal_scale: 2.4
  # Learn proposal covariance matrix?
  # Can be a dictionary, specifying when to update covariance matrix, with same options as check
  # e.g. to update proposal when Gelman-Rubin is between 0.03 and 0.1: {'max_eigen_gr': 0.1, 'min_eigen_gr': 0.03}
  learn: True
  # Use dragging ("integrating out" fast parameters)
  drag: False

run:
  # Thin samples by this factor
  thin_by: 1

pocomc

class: PocoMCSampler

info:
  version: 0.0.1
  date: 05/04/2022
  maintainer: Arnaud de Mattia
  description: Wrapper for PocoMC sampler
  url: https://github.com/minaskar/pocomc
  doi: [10.1093/mnras/stac2272, 10.48550/arXiv.2207.05660]
  long_description: 'Preconditioned Monte Carlo method for accelerated Bayesian inference'

init:
  # Number of walkers, defaults to 2 * max((int(2.5 * ndim) + 1) // 2, 2)
  nwalkers: None
  # The threshold value for the (normalised) proposal scale parameter below which
  # normalising flow preconditioning (NFP) is enabled.
  # Default is threshold = 1.0, meaning that NFP is used all the time
  threshold: 1.0
  # Whether to scale the distribution of particles to
  # have zero mean and unit variance.
  scale: True
  # Whether to rescale the distribution of particles to
  # have zero mean and unit variance in every iteration.
  rescale: False
  # Use a diagonal covariance matrix when rescaling instead of a full covariance.
  diagonal: True
  # Configuration of the normalizing flow
  flow_config: None
  # Configuration for training the normalizing flow
  train_config: None