psiphy.mcmc

psiphy.mcmc.importance_sampling.tqdm(iterable, **kwargs)[source]

psiphy.mcmc.importance_sampling.eff_sample_size(weights)[source]

Effective sample size (ESS) of a weighted particle set.

ESS = 1 / sum(w_i^2) for normalised weights w_i.

Parameters:: weights (array-like, shape (n,)) – Raw or normalised importance weights.
Returns:: ESS in [1, n].
Return type:: float

psiphy.mcmc.importance_sampling.importance_sampling(log_target_fn, proposal_samples, log_proposal_fn=None, vectorized=False)[source]

Basic importance sampling.

Parameters:

log_target_fn (callable) – Log of the (unnormalised) target density. Signature: log_target_fn(theta) -> float
proposal_samples (np.ndarray, shape (n, d)) – Samples already drawn from the proposal distribution.
log_proposal_fn (callable or None) – Log of the proposal density evaluated at each sample. None assumes a uniform (constant) proposal.
vectorized (bool) – If True, log_target_fn (and log_proposal_fn) accept a batch of shape (n, d) and return shape (n,), so all particles are evaluated in one call.

Returns:

samples : (n, d) final weighted particles. weights : (n,) normalised weights, aligned with samples. equal_weight_samples : (n, d) particles resampled by weight; use these where unweighted draws are needed. ess : effective sample size of weights. log_evidence : log marginal likelihood estimate, or None. samples are the proposal samples themselves. log_evidence is given only when log_proposal_fn is supplied (it needs the normalised proposal density).

Return type:

dict

Examples

>>> prior = UniformPrior([(-3, 3), (0.1, 3)])
>>> res = importance_sampling(log_posterior, prior.sample(2000), log_proposal_fn=prior.log_prob)
>>> np.average(res["samples"], weights=res["weights"], axis=0)

psiphy.mcmc.importance_sampling.sequential_importance_sampling(log_target_fn, prior, n_particles=500, n_steps=5, kernel='EmpiricalCovariance', ess_threshold=0.5, verbose=True, vectorized=False)[source]

Sequential importance sampling with MCMC move kernel.

Each step re-weights the current particles under the target and applies a Gaussian random-walk MCMC move to diversify them.

Parameters:

log_target_fn (callable) – Log of the (unnormalised) target.
prior (UniformPrior or object with .sample(n) and .log_prob(theta)) – Used to draw the initial particle set.
n_particles (int)
n_steps (int) – Number of SIS iterations.
kernel (str or sklearn covariance estimator) – Kernel used to estimate the MCMC move covariance. Accepts 'EmpiricalCovariance' or 'LedoitWolf'.
ess_threshold (float) – Resample when ESS / n_particles falls below this fraction.
verbose (bool)
vectorized (bool) – If True, log_target_fn accepts a batch (n, d) and returns (n,).

Returns:

samples : (n, d) final weighted particles. weights : (n,) normalised weights, aligned with samples. equal_weight_samples : (n, d) particles resampled by weight; use these where unweighted draws are needed. ess : effective sample size of weights. log_evidence : log marginal likelihood estimate, or None. log_evidence is always None for SIS.

Return type:

dict

psiphy.mcmc.importance_sampling.sequential_monte_carlo(log_likelihood_fn, prior, n_particles=1000, n_steps=10, kernel='EmpiricalCovariance', ess_threshold=0.5, verbose=True, vectorized=False)[source]

Sequential Monte Carlo via likelihood tempering.

Anneals from the prior to the posterior through a sequence of tempered targets:

p_t(theta) ∝ prior(theta) * likelihood(theta)^beta_t

where beta_t = (t / n_steps)^2 for t = 0, …, n_steps.

Parameters:

log_likelihood_fn (callable) – Log-likelihood. log_likelihood_fn(theta) -> float
prior (UniformPrior or object with .sample(n) and .log_prob(theta))
n_particles (int)
n_steps (int) – Number of temperature steps.
kernel (str or sklearn covariance estimator)
ess_threshold (float) – Fraction of n_particles below which resampling + MCMC move fires.
verbose (bool)
vectorized (bool) – If True, log_likelihood_fn accepts a batch (n, d) and returns (n,).

Returns:

samples : (n, d) final weighted particles. weights : (n,) normalised weights, aligned with samples. equal_weight_samples : (n, d) particles resampled by weight; use these where unweighted draws are needed. ess : effective sample size of weights. log_evidence : log marginal likelihood estimate, or None. log_evidence is the standard SMC estimate (sum over temperature steps of the log mean incremental weight); it requires a normalised log_likelihood_fn and prior. Also returns log_likelihood (n,).

Return type:

dict

Notes

Effective sample size at each step is reported when verbose=True.

psiphy.mcmc.importance_sampling.SMC(*args, **kwargs)[source]: Deprecated alias of sequential_monte_carlo(); will be removed in psiphy 0.3.

Ensemble Metropolis-Hastings style sampling with a pluggable backend.

Currently supported backends: "emcee".

Examples

>>> sampler = MHSampler(log_prob_fn, n_walkers=32, backend="emcee")
>>> sampler.run(theta_init, n_steps=1000)
>>> chain = sampler.get_chain(flat=True, discard=200, thin=10)

class psiphy.mcmc.mh_sampler.MHSampler(log_prob_fn, n_walkers=32, backend='emcee', filename=None, param_names=None, moves=None, pool=None, **kwargs)[source]

Wrapper around ensemble MCMC samplers.

Parameters:

log_prob_fn (callable) – log_prob_fn(theta) -> float; return -np.inf outside the support.
n_walkers (int) – Number of walkers in the ensemble.
backend (str) – Sampling library. Only "emcee" is available.
filename (str or None) – If given, the chain is stored in this HDF5 file (requires h5py) and can be reloaded with load() or extended with further run() calls.
param_names (list of str, optional) – Names used by summary().
moves (emcee move or list, optional) – Proposal moves forwarded to emcee.EnsembleSampler.
pool (object with a map method, optional) – For parallel log-prob evaluation.
**kwargs – Extra keyword arguments forwarded to emcee.EnsembleSampler.

run(theta_init=None, n_steps=1000, init_scale=0.001, reset=False, progress=False)[source]

Advance the chain by n_steps iterations.

Parameters:

theta_init (array-like, shape (n_dim,) or (n_walkers, n_dim), or None) – Starting point. A single point is jittered into a small Gaussian ball of relative width init_scale. None continues from the last state (of this run, or of the chain stored in filename).
n_steps (int) – Number of iterations to add.
reset (bool) – Discard any chain already stored in filename.
progress (bool) – Show emcee’s progress bar.

Return type:

self

classmethod load(filename, log_prob_fn=None, **kwargs)[source]

Reload a chain saved to an HDF5 file.

Pass log_prob_fn to be able to continue sampling with run(None, n_steps).

property n_steps

get_chain(flat=False, discard=0, thin=1)[source]

Chain of shape (n_steps, n_walkers, n_dim), or (n_samples, n_dim) if flat.

discard="auto" drops three times the longest autocorrelation time.

get_log_prob(flat=False, discard=0, thin=1)[source]

autocorr_time(quiet=True)[source]: Integrated autocorrelation time per parameter.

property acceptance_fraction

summary(discard='auto', thin=1, print_out=True)[source]

Per-parameter posterior summary.

Returns:: name -> {mean, std, median, lo68, hi68, tau} plus "acceptance_fraction" and "n_samples".
Return type:: dict

Nested sampling with interchangeable backends.

Supported backends: "dynesty" (static nested sampling) and "nautilus" (neural-network-boosted importance nested sampling).

Examples

>>> sampler = NestedSampler(log_likelihood_fn, prior_transform, n_dim=3, backend="dynesty")
>>> sampler.run()
>>> result = sampler.get_result()   # log_evidence, samples, weights

class psiphy.mcmc.nested_sampler.NestedSampler(log_likelihood_fn, prior_transform, n_dim=None, backend='dynesty', n_live=500, param_names=None, pool=None, seed=None, **kwargs)[source]

Backend-agnostic nested sampling wrapper.

Parameters:

log_likelihood_fn (callable) – log_likelihood_fn(theta) -> float.
prior_transform (callable or psiphy.utils.prior.BasePrior) – Maps a point in the unit hypercube to parameter space. A prior object is accepted directly; its prior_transform, ndim and names are used.
n_dim (int, optional) – Number of parameters. Inferred when prior_transform is a prior object.
backend ({"dynesty", "nautilus"})
n_live (int) – Number of live points.
param_names (list of str, optional)
pool (object with a map method, optional) – For parallel likelihood evaluation. With dynesty also pass queue_size.
seed (int, optional) – Random seed for the sampler and for equal-weight resampling.
**kwargs – Extra keyword arguments for the backend’s sampler constructor (dynesty.NestedSampler or nautilus.Sampler).

run(dlogz=0.01, progress=False, **run_kwargs)[source]

Run the sampler to convergence.

Parameters:

dlogz (float) – Stopping tolerance. For dynesty this is the remaining evidence dlogz; for nautilus it is used as f_live (fraction of evidence in the live set).
progress (bool) – Print backend progress.
**run_kwargs – Extra keyword arguments for run_nested (dynesty) or run (nautilus).

Return type:

self

get_result()[source]

Unified results dictionary (same keys for every backend).

Returns:: log_evidence, log_evidence_err (None for nautilus), samples (n, n_dim), normalised weights (n,), log_likelihood (n,), equal_weight_samples, ess, backend.
Return type:: dict

to_getdist(names=None, labels=None, label=None)[source]: Weighted getdist.MCSamples of the posterior.