psiphy.sbi

psiphy.sbi.rejection_abc.tqdm(iterable, **kwargs)

class psiphy.sbi.rejection_abc.rABC(simulator, distance, observation, prior, bounds, N=100, eps=0.1, verbose=True)

sample_prior(kk)

run(N=None)

class psiphy.sbi.sre.emceeSRE

ratio_estimation(classifier, X, y, test_size=0.1, normalizer=None)

set_obs(obs)

set_simulator(simulator)

learn_logL_with_classifier(classifier, mins, maxs, Nsamples=None, test_size=0.1)

LFIRE — Likelihood-Free Inference by Ratio Estimation (Thomas et al. 2022, Bayesian Analysis 17, 1).

For a parameter value theta, simulated data from p(x | theta) and from the marginal p(x) (prior predictive) are separated with a logistic-regression classifier. Its log-odds estimate the log density ratio

log r(x, theta) = log p(x | theta) - log p(x),

so p(theta | x_obs) ∝ prior(theta) * r(x_obs, theta). Each evaluation costs n_theta simulations plus one classifier fit, so LFIRE either evaluates a grid directly (low dimensions) or learns a GP surrogate of the log ratio with Bayesian optimisation and samples that, as in BOLFI.

Examples

>>> lfire = LFIRE(simulator, x_obs, {'mu': (0, 3), 'sigma': (0.1, 2.5)}, summary=my_summaries)
>>> lfire.learn_posterior(n_calls=60)
>>> samples = lfire.sample_posterior(5000)

class psiphy.sbi.lfire.LFIRE(simulator, observation, prior, n_m=500, n_theta=500, summary=None, penalty='l1', classifier=None, acquisition='UCB', clip_log_ratio=20.0, bounds=None, param_names=None, n_jobs=1, random_state=None, verbose=True)

Likelihood-free inference by ratio estimation.

Parameters:

• simulator (callable) – simulator(theta) -> x (array-like).

• observation (array-like) – Observed data, in the same format as the simulator output.

• prior (list of (low, high), dict name -> (low, high), or prior object)

• n_m (int) – Simulations from the marginal (prior predictive), drawn once and reused.

• n_theta (int) – Simulations at each parameter value evaluated.

• summary (callable, optional) – summary(x) -> 1-D features used by the classifier. Defaults to the flattened data. The log ratio is linear in these features, so include nonlinear terms (e.g. squares) when the data distribution changes shape with theta.

• penalty ({"l1", "l2"}) – Regularisation of the default classifier (cross-validated). "l1" follows the paper and selects among many candidate summaries; "l2" is several times faster when there are few summaries.

• classifier (sklearn classifier, optional) – Replaces the default; must provide decision_function (log-odds). It is cloned for every evaluation.

• acquisition ({"UCB", "EI", "PI"} or callable) – Acquisition for learn_posterior() (maximising the log ratio).

• clip_log_ratio (float or None) – The GP surrogate is fitted to max(log_r, max(log_r) - clip_log_ratio). Far from the posterior the classes separate perfectly and the log ratio becomes hugely negative (limited only by regularisation); without the floor those values dominate the GP and widen the posterior. Points below the floor carry less than exp(-clip_log_ratio) of the peak posterior density. None disables it.

• bounds (list of (low, high), optional) – Search box; defaults to the prior bounds.

• param_names (list of str, optional)

• n_jobs (int) – Workers for the classifier’s cross-validation.

• random_state (int, optional) – Seed for the classifier.

• verbose (bool)

simulate_marginal()

Simulate the marginal (prior-predictive) set once; reused by every evaluation.

log_ratio(theta, return_score=False)

Estimate log r(x_obs, theta) = log p(x_obs | theta) - log p(x_obs).

Parameters:

• theta (array-like, shape (n_dim,))

• return_score (bool) – Also return the classifier’s training accuracy. Values near 0.5 mean data at theta are indistinguishable from the marginal.

Return type:

float, or (float, float) with return_score

run_grid(n_grid=20, thetas=None)

Evaluate the log ratio directly at every grid point (no surrogate).

Costs n_theta simulations per point, so this suits one or two parameters. Results: self.thetas, self.log_ratios and the normalised grid posterior self.posterior (also used by psiphy.plotting.corner.plot_lfire()).

Return type:

thetas, posterior

log_likelihood(theta)

Surrogate log ratio (the log-likelihood up to a constant).

learn_posterior(n_calls=60, n_initial_points=None, batch_size=1, filename=None, reset_model=False, random_state=None, n_candidates=None, **acq_kwargs)

Learn a GP surrogate of the log ratio with Bayesian optimisation.

Evaluations concentrate where the log ratio is high, i.e. in the posterior bulk. Afterwards use sample_posterior() or posterior_grid().

Parameters:

• n_calls (int) – New log-ratio evaluations (each costs n_theta simulations).

• n_initial_points (int, optional) – Latin-hypercube initial design size; default 5 * n_dim.

• batch_size (int) – Points proposed per round (kriging believer).

• filename (str, optional) – Pickle file to resume from (unless reset_model) and save to. It also stores the marginal simulations.

• reset_model (bool)

• random_state (int, optional)

• n_candidates (int, optional)

• **acq_kwargs – e.g. kappa for UCB (default 1.96), xi for EI/PI.

Return type:

self

posterior_grid(n_grid=50)

Normalised posterior of the surrogate on a grid.

Sets self.thetas and self.posterior (for psiphy.plotting.corner.plot_lfire()).

Return type:

thetas, posterior

load_model(filename)

log_posterior(theta)

sample_posterior(n_samples=5000, method='MCMC', **kwargs)

Draw equal-weight samples from the surrogate posterior.

Parameters:

• n_samples (int)

• method ({"MCMC", "NestedSampling", "IS", "SMC"}) – See psiphy.sbi._surrogate.sample_posterior() for the options of each method.

Returns:

The full sampler output is kept in self.sampler_info.

Return type:

np.ndarray, shape (n_samples, n_dim)

save_model(filename)

Store the evaluated points (the GP is refitted on load).

BOLFI — Bayesian Optimisation for Likelihood-Free Inference.

A Gaussian-process surrogate of the discrepancy d(theta) between simulated and observed data is learned with Bayesian optimisation (Gutmann & Corander 2016). The synthetic log-likelihood log K(d_hat(theta)) is then sampled with any of the psiphy samplers (emcee, nested sampling, importance sampling, SMC).

Examples

>>> bolfi = BOLFI(distance_fn, {'mu': (-5, 5), 'sigma': (0.1, 5)})
>>> bolfi.learn_likelihood(n_calls=100)
>>> samples = bolfi.sample_posterior(5000, method='MCMC')

class psiphy.sbi.bolfi.BOLFI(distance, prior, distance_kernel='exp', acquisition='LCB', learn_log_dist=False, n_repeats=1, bounds=None, param_names=None, verbose=True)

Bayesian Optimisation for Likelihood-Free Inference.

Parameters:

• distance (callable) – distance(theta) -> float: simulate at theta and return the discrepancy to the observed data.

• prior (list of (low, high), dict name -> (low, high), or prior object) – Bounds are converted to a UniformPrior. A prior object must expose mins/maxs unless bounds is given.

• distance_kernel ({"exp", "gaussian"} or callable) – Maps the discrepancy to a log-likelihood. A callable must have signature log_kernel(d) -> log L.

• acquisition ({"LCB", "EI", "PI"} or callable) – Acquisition function choosing where to simulate next; see psiphy.utils.bayesian_opt. "LCB" follows Gutmann & Corander (2016).

• learn_log_dist (bool) – Model log d instead of d (helps when d spans orders of magnitude).

• n_repeats (int) – Distance evaluations (simulations) averaged per parameter value. With a stochastic simulator each d(theta) is noisy; averaging n of them cuts the noise variance by n and makes the GP surrogate, and the posterior widths, more stable from run to run. The simulations at one theta are independent, so they can run in parallel. The extra noise in a single-simulation d adds a roughly constant offset that cancels in the normalised posterior.

• bounds (list of (low, high), optional) – Search space for the optimisation; defaults to the prior bounds.

• param_names (list of str, optional)

• verbose (bool)

property n_simulations

Distance evaluations (simulations) used so far.

predict_distance(theta, return_std=False)

GP mean (and std) of the modelled discrepancy (log d if learn_log_dist).

log_likelihood(theta)

Synthetic log-likelihood log K(d_hat(theta)).

learn_likelihood(n_calls=100, n_initial_points=None, batch_size=1, filename=None, reset_model=False, random_state=None, n_candidates=None, **acq_kwargs)

Learn the discrepancy surrogate with Bayesian optimisation.

Parameters:

• n_calls (int) – Number of new parameter values evaluated (including initial design points); each costs n_repeats distance evaluations.

• n_initial_points (int, optional) – Size of the Latin-hypercube initial design; default 5 * n_dim. Points already loaded from filename count towards it.

• batch_size (int) – Points proposed per optimisation round (kriging believer). The points of one batch are independent, so their simulations can run in parallel; they are currently evaluated in sequence.

• filename (str, optional) – Pickle file to resume from (unless reset_model) and save to.

• reset_model (bool) – Ignore an existing filename.

• random_state (int, optional)

• n_candidates (int, optional) – Random candidates scored per pick (default 1000 * n_dim).

• **acq_kwargs – Acquisition options, e.g. kappa (LCB) or xi (EI/PI).

Return type:

self

save_likelihood_model(filename)

Store the evaluated points (the GP is refitted on load).

load_likelihood_model(filename)

load_model(filename)

log_posterior(theta)

sample_posterior(n_samples=5000, method='MCMC', **kwargs)

Draw equal-weight samples from the surrogate posterior.

Parameters:

• n_samples (int)

• method ({"MCMC", "NestedSampling", "IS", "SMC"}) – See psiphy.sbi._surrogate.sample_posterior() for the options of each method.

Returns:

The full sampler output is kept in self.sampler_info.

Return type:

np.ndarray, shape (n_samples, n_dim)

save_model(filename)

Store the evaluated points (the GP is refitted on load).

psiphy.sbi.abc_gp.tqdm(iterable, **kwargs)

class psiphy.sbi.abc_gp.ABC_gpL(simulator, distance, obs, theta_sampler=None, theta_range={}, n_train_init=100, mcmc_sampler=None, mcmc_sampler_info=None)

prepare_distance_model(model=None, kernel=None, alpha=1e-10, optimizer='fmin_l_bfgs_b', n_restarts_optimizer=3)

create_dataset(n_train=None)

learn_distance(n_train=None)

learn_distance_BO(tol=0.01)

run_mcmc(n_samples=None)