Quickstart

Short examples of the main tools. Each runs in a few seconds; the Tutorials explain the methods and their options in detail.

Define a prior

All samplers take the same prior objects. Named bounds give parameter names to the outputs.

from psiphy.utils.prior import UniformPrior, GaussianPrior

prior = UniformPrior({'slope': (-3, 3), 'intercept': (0, 10)})
prior.sample(5)                  # (5, 2) array of draws
prior.log_prob([0.5, 4.0])       # log density

MCMC

Sample a posterior with an ensemble MCMC sampler (requires the mcmc extra).

from psiphy.mcmc import MHSampler
from psiphy.toy_models import NoisyLine

line = NoisyLine()
y_obs = line.simulate([-0.8, 5.0])

def log_posterior(theta):
    return prior.log_prob(theta) + line.log_prob(theta, y_obs)

sampler = MHSampler(log_posterior, n_walkers=32, param_names=prior.names)
sampler.run(prior.sample(32), n_steps=2000)
samples = sampler.get_chain(flat=True, discard='auto')
sampler.summary()

Nested sampling and evidence

Nested sampling returns the Bayesian evidence for model comparison together with posterior samples (requires the mcmc extra).

from psiphy.mcmc import NestedSampler

ns = NestedSampler(lambda t: line.log_prob(t, y_obs), prior, backend='nautilus')
result = ns.run().get_result()
result['log_evidence'], result['samples'], result['weights']

Likelihood-free inference

When only a simulator is available, BOLFI learns the posterior from a small number of simulations. distance simulates at theta and compares summary statistics of the simulated and observed data — here the least-squares slope and intercept, scaled by their expected scatter, so that exp(-distance / 2) behaves like a likelihood.

import numpy as np
from psiphy.sbi.bolfi import BOLFI

A = np.column_stack([line.x, np.ones_like(line.x)])
fit = lambda y: np.linalg.lstsq(A, y, rcond=None)[0]       # (slope, intercept)
inv_cov = A.T @ A / line.noise_sigma**2                      # inverse covariance of the fit
p_obs = fit(y_obs)

def distance(theta):
    dp = fit(line.simulate(theta)) - p_obs
    return dp @ inv_cov @ dp

bolfi = BOLFI(distance, prior, verbose=False)
bolfi.learn_likelihood(n_calls=60)
bolfi_samples = bolfi.sample_posterior(2000)

Fisher forecast

Forecast parameter constraints from a simulator with signature simulator(theta, seed, noise_seed).

from psiphy.forecasting.fisher import FisherMatrix

def simulator(theta, seed, noise_seed):
    rng = np.random.default_rng(1000 * seed + noise_seed)
    noise = rng.normal(0, line.noise_sigma, line.x.size) if noise_seed else 0.0
    return theta[0] * line.x + theta[1] + noise

fm = FisherMatrix(simulator, theta_fid=[-0.8, 5.0], param_names=['slope', 'intercept'],
                  verbose=False)
fm.estimate_covariance(seeds=range(1, 21), noise_seeds=range(1, 11))
fm.estimate_derivatives(delta_frac=0.01, seeds=range(1, 3), noise_seed=0)
fm.compute()
fm.sigma_1d()                    # marginalised 1-sigma forecasts

Plot and summarise

from psiphy.plotting.diagnostics import plot_triangle, print_chain_stats

plot_triangle(samples, labels=['m', 'b'], true_values=[-0.8, 5.0], engine='getdist')
print_chain_stats(samples, labels=['m', 'b'])

engine='corner' (the default) needs the corner package.