psiphy.utils

psiphy.utils.distances.kl_divergence(p_probs, q_probs)

“KL (p || q)

psiphy.utils.distances.js_distance(p_probs, q_probs)

Kernel Density Estimation

Written around the KDE code in scikit-learn package.

psiphy.utils.kernel_density.bandwidth_kde_silverman(X)

h_ii = (4/(d+2))**(1/(d+4)) n**(-1/(d+4)) sigma_i

psiphy.utils.kernel_density.bandwidth_kde_scott(X)

h_ii = n**(-1/(d+4)) sigma_i

psiphy.utils.kernel_density.bandwidth_kdeCV(X, kde=None, bw=array([0.01, 0.01151395, 0.01325711, 0.01526418, 0.01757511, 0.0202359, 0.02329952, 0.02682696, 0.03088844, 0.0355648, 0.04094915, 0.04714866, 0.05428675, 0.06250552, 0.07196857, 0.08286428, 0.09540955, 0.10985411, 0.12648552, 0.14563485, 0.16768329, 0.19306977, 0.22229965, 0.25595479, 0.29470517, 0.33932218, 0.39069399, 0.44984327, 0.51794747, 0.59636233, 0.68664885, 0.79060432, 0.91029818, 1.04811313, 1.20679264, 1.38949549, 1.59985872, 1.84206997, 2.12095089, 2.44205309, 2.8117687, 3.23745754, 3.72759372, 4.29193426, 4.94171336, 5.68986603, 6.55128557, 7.54312006, 8.68511374, 10.]), cv=1, verbose=True, kernel='gaussian', metric='euclidean', atol=0, leaf_size=40)

Estimate the bandwidth using cross validation.

psiphy.utils.kernel_density.cross_val_loss_kdeCV_loo(h)=\int \hat{f}^2_n(x)dx - 2\int \hat{f}(x)f(x)dx \hat{J}(h) = \int \hat{f}^2_n(x)dx - \frac{2}{n}\sum \hat{f}_{-i}(x_i)

psiphy.utils.kernel_density.cross_val_loss_kdeCV_kFold(h)=\int \hat{f}^2_n(x)dx - 2\int \hat{f}(x)f(x)dx \hat{J}(h) = \int \hat{f}^2_n(x)dx - \frac{2}{n}\sum \hat{f}_{-i}(x_i)

psiphy.utils.helpers.loading_verbose(string)

psiphy.utils.helpers.spherical_to_cartesian(theta)

psiphy.utils.helpers.save_model(filename, data)

psiphy.utils.helpers.load_model(filename)

psiphy.utils.sampling_space.LH_sampling(n_params=2, samples=10, mins=0, maxs=1, outfile=None)

Parameters:

• n_params (int) – Number of parameters.

• samples (int) – Number of points in the parameter space.

• mins (float or list) – Lower and upper bounds. A single float is used for all parameters.

• maxs (float or list) – Lower and upper bounds. A single float is used for all parameters.

• outfile (str or None) – Text file to save the points to. None does not save them.

Returns:

The parameter values.

Return type:

np.ndarray, shape (samples, n_params)

psiphy.utils.sampling_space.MC_sampling(n_params=2, samples=10, mins=0, maxs=1, outfile=None)

Parameters:

• n_params (int) – Number of parameters.

• samples (int) – Number of points in the parameter space.

• mins (float or list) – Lower and upper bounds. A single float is used for all parameters.

• maxs (float or list) – Lower and upper bounds. A single float is used for all parameters.

• outfile (str or None) – Text file to save the points to. None does not save them.

Returns:

The parameter values.

Return type:

np.ndarray, shape (samples, n_params)

psiphy.utils.sampling_space.MCS_nsphere(n_params=2, samples=10, mins=0, maxs=1, outfile=None)

Parameters:

• n_params (int) – Number of parameters.

• samples (int) – Number of points in the parameter space.

• mins (float or list) – Lower and upper bounds. A single float is used for all parameters.

• maxs (float or list) – Lower and upper bounds. A single float is used for all parameters.

• outfile (str or None) – Text file to save the points to. None does not save them.

Returns:

The parameter values.

Return type:

np.ndarray, shape (samples, n_params)

psiphy.utils.sampling_space.LHS_nsphere(n_params=2, samples=10, mins=0, maxs=1, outfile=None)

Parameters:

• n_params (int) – Number of parameters.

• samples (int) – Number of points in the parameter space.

• mins (float or list) – Lower and upper bounds. A single float is used for all parameters.

• maxs (float or list) – Lower and upper bounds. A single float is used for all parameters.

• outfile (str or None) – Text file to save the points to. None does not save them.

Returns:

The parameter values.

Return type:

np.ndarray, shape (samples, n_params)

Framework-agnostic prior specifications.

Every prior exposes the same interface so it can be passed to the importance sampling / SMC routines, MCMC samplers, nested samplers (via BasePrior.prior_transform()) and, through the conversion helpers, to scipy, torch and sbi.

Examples

>>> prior = UniformPrior({'mu': (-5, 5), 'sigma': (0.1, 5)})
>>> theta = prior.sample(1000)                  # (1000, 2)
>>> lp = prior.log_prob(theta)                  # (1000,)
>>> theta = prior.prior_transform(np.random.rand(2))   # unit cube -> theta

class psiphy.utils.prior.BasePrior

Common interface for all priors.

Subclasses implement sample(), log_prob() and, where it exists in closed form, prior_transform().

log_prob accepts a single point of shape (ndim,) (returns a float) or a batch of shape (n, ndim) (returns an array of shape (n,)).

names = None

property ndim

sample(n=1)

Draw n samples; returns shape (n, ndim).

log_prob(theta)

prior_transform(u)

Map points from the unit hypercube to parameter space (for nested sampling).

in_support(theta)

to_scipy()

to_torch()

to_sbi()

Return a distribution usable as an sbi prior (defaults to to_torch()).

class psiphy.utils.prior.UniformPrior(bounds)

Box-uniform prior.

Parameters:

bounds (list of (low, high) tuples, or dict mapping name -> (low, high))

Examples

>>> prior = UniformPrior([(0, 1), (-2, 2)])
>>> prior = UniformPrior({'mu': (-5, 5), 'sigma': (0.1, 5)})

property ndim

sample(n=1)

Draw n samples; returns shape (n, ndim).

log_prob(theta)

prior_transform(u)

Map points from the unit hypercube to parameter space (for nested sampling).

to_scipy()

List of independent frozen scipy.stats.uniform, one per parameter.

to_torch()

to_sbi()

Return a distribution usable as an sbi prior (defaults to to_torch()).

class psiphy.utils.prior.LogUniformPrior(bounds)

Independent log-uniform (Jeffreys-like) prior: log(theta) is uniform.

Parameters:

bounds (list of (low, high) tuples, or dict mapping name -> (low, high).) – All bounds must be strictly positive.

property ndim

sample(n=1)

Draw n samples; returns shape (n, ndim).

log_prob(theta)

prior_transform(u)

Map points from the unit hypercube to parameter space (for nested sampling).

to_scipy()

List of independent frozen scipy.stats.loguniform, one per parameter.

to_torch()

class psiphy.utils.prior.GaussianPrior(mean, cov, names=None)

Multivariate Gaussian prior.

Parameters:

• mean (array-like, shape (ndim,))

• cov (float, array-like of shape (ndim,) (diagonal variances), or (ndim, ndim))

• names (list of str, optional)

property ndim

sample(n=1)

Draw n samples; returns shape (n, ndim).

log_prob(theta)

prior_transform(u)

Map points from the unit hypercube to parameter space (for nested sampling).

to_scipy()

to_torch()

class psiphy.utils.prior.MixturePrior(components, weights=None)

Weighted mixture of priors with the same dimensionality.

Parameters:

• components (list of BasePrior)

• weights (array-like, optional) – Mixture weights (normalised internally). Equal weights by default.

property ndim

sample(n=1)

Draw n samples; returns shape (n, ndim).

log_prob(theta)

Bayesian-optimisation helpers built on scikit-learn Gaussian processes.

Acquisition functions return values to maximise. Each takes the candidate points, a fitted GaussianProcessRegressor, the incumbent y_best and a minimize flag stating whether the underlying objective is minimised (e.g. a BOLFI discrepancy) or maximised (e.g. an LFIRE posterior).

Examples

>>> gpr = GaussianProcessRegressor(Matern(nu=2.5), normalize_y=True).fit(X, y)
>>> x_next = propose_location(gpr, bounds, acquisition="EI", X_sample=X, y_sample=y)
>>> batch = propose_location(gpr, bounds, acquisition="LCB", X_sample=X, y_sample=y, batch_size=4)

psiphy.utils.bayesian_opt.expected_improvement(X, gpr, y_best, minimize=True, xi=0.01)

Expected improvement over y_best.

psiphy.utils.bayesian_opt.probability_of_improvement(X, gpr, y_best, minimize=True, xi=0.01)

Probability of improving on y_best by at least xi.

psiphy.utils.bayesian_opt.confidence_bound(X, gpr, y_best=None, minimize=True, kappa=1.96)

GP confidence bound: -(mu - kappa*sd) (LCB) when minimising, mu + kappa*sd (UCB) when maximising. Larger kappa explores more.

psiphy.utils.bayesian_opt.get_acquisition(acquisition)

Look up an acquisition by name, or pass a callable through.

psiphy.utils.bayesian_opt.propose_location(gpr, bounds, acquisition='EI', X_sample=None, y_sample=None, batch_size=1, minimize=True, n_candidates=None, n_restarts=5, random_state=None, **acq_kwargs)

Propose the next point(s) to evaluate.

Parameters:

• gpr (fitted sklearn GaussianProcessRegressor)

• bounds (array-like, shape (d, 2)) – Search box, in the same coordinates the GP was fitted in.

• acquisition ({"EI", "PI", "LCB", "UCB"} or callable) – A callable must have signature acq(X, gpr, y_best, minimize, **acq_kwargs) and return values to maximise.

• X_sample (array-like) – Evaluated points. Needed for EI/PI (incumbent) and for batch_size > 1.

• y_sample (array-like) – Evaluated points. Needed for EI/PI (incumbent) and for batch_size > 1.

• batch_size (int) – Number of points to return. Points after the first are chosen by the kriging-believer heuristic: each chosen point is added to the GP with its predicted mean as a fantasy observation, which shrinks the uncertainty there and pushes the next pick elsewhere. Use this to run several simulations in parallel per round.

• minimize (bool) – Whether the objective modelled by the GP is minimised.

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

• n_restarts (int) – Best candidates polished with L-BFGS-B.

• random_state (int or numpy Generator, optional)

• **acq_kwargs – Passed to the acquisition (xi for EI/PI, kappa for LCB/UCB).

Return type:

np.ndarray, shape (batch_size, d)

psiphy.utils.bayesian_opt.check_inside_nsphere(X, bound_min, bound_max)

psiphy.utils.bayesian_opt.draw_nsphere_surface(bound_min, bound_max, N=100)