gpmp.mcmc module

The mcmc module provides posterior covariance-parameter sampling helpers and lower-level samplers. Samples, particles, and log-probability traces are backend-native gpmp.num objects unless a function explicitly documents conversion. For GP parameter posterior exploration, prefer the sample_from_selection_criterion_* functions: they consume an optimization info object returned by gpmp.kernel selection functions, or a standalone negative log-posterior callable.

Target convention

The posterior helpers treat the selection criterion J(theta) as a negative log-density. The target log-density is -J(theta) for NUTS and -J(theta) / temperature for tempered MH/SMC. sampling_box imposes hard bounds on the target support. init_box is only an initialization box.

Initialization convention

For MH and NUTS, pass param_initial_states explicitly or set random_init=True with init_box. For SMC, init_box is mandatory. The parameter dimension is inferred from info.covparam, initial states, or the supplied box.

Experimental samplers are intentionally not documented here until their API is stable.

Posterior parameter sampling

gpmp.mcmc.sample_from_selection_criterion_mh(info: object = None, selection_criterion: Callable = None, param_initial_states: ndarray[tuple[Any, ...], dtype[floating]] = None, random_init: bool = False, init_box: list = None, sampling_box: list = None, temperature: float = 1.0, n_steps_total: int = 10000, burnin_period: int = 4000, n_chains: int = 2, n_pool: int = 2, silent: bool = False, show_progress: bool = True, plot_chains: bool = True, plot_empirical_distributions: bool = True)

Sample from a parameter posterior with adaptive Metropolis-Hastings.

The target log-density is defined as log_target(theta) = -J(theta) / temperature, where J is the selection criterion. When sampling_box is provided, points outside the box are assigned -inf log-density.

Parameters:

• info (object, optional) – Optimization info object carrying at least covparam and a criterion callable (typically selection_criterion_nograd or selection_criterion). Provide exactly one of info and selection_criterion.

• selection_criterion (callable, optional) – Negative log-posterior (or any negative objective) with signature f(theta) -> scalar. Provide exactly one of info and selection_criterion.

• param_initial_states (array_like, optional) – Initial parameter states. Accepted shapes are scalar, (dim,), (n_chains, dim), (1, dim), and (dim, n_chains).

• random_init (bool, default=False) – If True, initialize chains uniformly in init_box.

• init_box (list, optional) – Initialization box [lower, upper] used only when random_init=True.

• sampling_box (list, optional) – Hard sampling bounds [lower, upper] used to truncate the target density during MH transitions.

• temperature (float, default=1.0) – Tempering factor applied to the criterion in log-space. Must be strictly positive. temperature=1.0 gives the nominal target; larger values flatten the target density.

• n_steps_total (int, default=10000) – Total number of MH steps per chain, including burn-in.

• burnin_period (int, default=4000) – Number of initial steps discarded in the returned samples.

• n_chains (int, default=2) – Number of MH chains.

• n_pool (int, default=2) – Chain pooling factor used by covariance adaptation.

• silent (bool, default=False) – If True, suppress MH textual diagnostics.

• show_progress (bool, default=True) – If True and silent=False, print global progress messages.

• plot_chains (bool, default=True) – If True, display MH chain traces after sampling.

• plot_empirical_distributions (bool, default=True) – If True, display marginal empirical distributions after sampling.

Returns:

• samples_post_burnin (ndarray) – Posterior samples with shape (n_chains, n_steps_total - burnin_period, dim).

• mh (MetropolisHastings) – Sampler instance containing full trajectories, acceptance history, and diagnostics.

gpmp.mcmc.sample_from_selection_criterion_nuts(info: object = None, selection_criterion: Callable = None, param_initial_states: ndarray[tuple[Any, ...], dtype[floating]] = None, random_init: bool = False, init_box: list = None, sampling_box: list = None, num_samples: int = 2000, num_warmup: int = 1000, n_chains: int = 2, target_accept: float = 0.8, max_depth: int = 10, delta_max: float = 1000.0, jitter: float = 0.0001, init_step_size: float | None = None, init_mass_diag: ndarray[tuple[Any, ...], dtype[floating]] = None, seed: int = None, progress: bool = True, verbose: int = 1, log_every: int = 50, options: NUTSOptions = None, plot_diagnostics: bool = False, diagnostics_window: int = 50, diagnostics_show: bool = True, diagnostics_save_dir: str = None)

Sample from a parameter posterior with the No-U-Turn Sampler (NUTS).

The target log-density is defined as log_prob(theta) = -J(theta), where J is the selection criterion. When sampling_box is provided, points outside the box are assigned -inf log-density.

Parameters:

• info (object, optional) – Optimization info object carrying at least covparam and a differentiable criterion callable. Provide exactly one of info and selection_criterion.

• selection_criterion (callable, optional) – Differentiable negative log-posterior with signature f(theta) -> scalar. Provide exactly one of info and selection_criterion.

• param_initial_states (array_like, optional) – Initial parameter states. Accepted shapes are scalar, (dim,), (n_chains, dim), (1, dim), and (dim, n_chains).

• random_init (bool, default=False) – If True, initialize chains uniformly in init_box.

• init_box (list, optional) – Initialization box [lower, upper] used only when random_init=True.

• sampling_box (list, optional) – Hard sampling bounds [lower, upper] used to truncate the target density during NUTS transitions.

• num_samples (int, default=2000) – Number of post-warmup samples per chain.

• num_warmup (int, default=1000) – Number of warmup iterations per chain.

• n_chains (int, default=2) – Number of chains.

• target_accept (float, default=0.8) – Target acceptance probability for dual averaging.

• max_depth (int, default=10) – Maximum tree depth per NUTS iteration.

• delta_max (float, default=1000.0) – Divergence threshold on Hamiltonian error.

• jitter (float, default=1e-4) – Lower clipping value for diagonal mass entries.

• init_step_size (float, optional) – User-provided initial step size. If None, a heuristic is used.

• init_mass_diag (array_like, optional) – User-provided initial diagonal mass vector of shape (dim,).

• seed (int, optional) – Random seed used by NUTS.

• progress (bool, default=True) – If True, display progress bars.

• verbose (int, default=1) – Logging verbosity level.

• log_every (int, default=50) – Logging frequency in iterations.

• options (NUTSOptions, optional) – Optional fully specified NUTS options object.

• plot_diagnostics (bool, default=False) – If True, generate NUTS diagnostic figures.

• diagnostics_window (int, default=50) – Moving-average window for diagnostics.

• diagnostics_show (bool, default=True) – If True, show diagnostic plots interactively.

• diagnostics_save_dir (str, optional) – Directory where diagnostic plots are saved when provided.

Returns:

• samples (ndarray) – Posterior samples with shape (n_chains, num_samples, dim).

• info_nuts (dict) – NUTS diagnostics and adaptation outputs (acceptance, divergences, tree depth, leapfrog counts, final step size, and final mass diagonal).

gpmp.mcmc.sample_from_selection_criterion_smc(info: object = None, selection_criterion: Callable = None, init_box: list = None, sampling_box: list = None, n_particles: int = 1000, initial_temperature: float = 1000000.0, final_temperature: float = 1.0, min_ess_ratio: float = 0.5, mh_steps: int = 20, max_stages: int = 50, debug: bool = False, plot_marginals: bool = False, plot_particles: bool = False)

Sample from a parameter posterior with Sequential Monte Carlo tempering.

The sampler targets a sequence of tempered densities of the form exp(-J(theta) / temperature) from initial_temperature down to final_temperature. When sampling_box is provided, points outside the box are assigned -inf log-density.

Parameters:

• info (object, optional) – Optimization info object carrying at least covparam and a criterion callable. Provide exactly one of info and selection_criterion.

• selection_criterion (callable, optional) – Negative log-posterior with signature f(theta) -> scalar. Provide exactly one of info and selection_criterion.

• init_box (list) – Mandatory initialization box [lower, upper] used to initialize SMC particles.

• sampling_box (list, optional) – Hard sampling bounds [lower, upper] used to truncate the target density during SMC.

• n_particles (int, default=1000) – Number of SMC particles.

• initial_temperature (float, default=1e6) – Initial tempering parameter (high temperature, easier target).

• final_temperature (float, default=1.0) – Final tempering parameter (target posterior scale).

• min_ess_ratio (float, default=0.5) – Minimum effective sample size ratio used by the tempering scheduler.

• mh_steps (int, default=20) – Number of MH move steps per SMC stage.

• max_stages (int, default=50) – Reserved stage budget parameter for compatibility.

• debug (bool, default=False) – If True, print detailed SMC diagnostics.

• plot_marginals (bool, default=False) – If True, plot empirical marginal distributions from SMC.

• plot_particles (bool, default=False) – Reserved plotting flag for compatibility.

Returns:

• particles (ParticlesSet) – Final particle set after tempering.

• smc_instance (SMC) – SMC driver instance containing execution logs and diagnostics.

gpmp.mcmc.get_log_target_values(mh: MetropolisHastings, *, discard_burnin: bool = False) → ndarray[tuple[Any, ...], dtype[floating]]

Return stored MH log-target values.

Parameters:

• mh (MetropolisHastings) – MH sampler instance returned by sample_from_selection_criterion_mh.

• discard_burnin (bool, default=False) – If True, return only values after mh.burnin_period.

Returns:

Log-target values with shape (n_chains, n_steps) when discard_burnin=False, or (n_chains, n_steps - burnin_period) when discard_burnin=True.

Return type:

ndarray

Raises:

ValueError – If mh.log_target_values is unavailable or burn-in indices are invalid.

Metropolis-Hastings

class gpmp.mcmc.MHOptions(dim: int = 1, n_chains: int = 1, symmetric: bool = True, target_acceptance: float = 0.3, acceptance_tol: float = 0.15, adaptation_method: str = 'Haario', proposal_distribution_param_init: ndarray[tuple[Any, ...], dtype[floating]] | None = None, adaptation_interval: int = 50, freeze_adaptation: bool = True, discard_burnin: bool = False, n_pool: int = 1, RM_adapt_factor: float = 1.0, RM_diminishing: bool = True, haario_adapt_factor_burnin_phase: float = 1.0, haario_adapt_factor_sampling_phase: float = 0.5, haario_initial_scaling_factor: float = 1.0, sliding_rate_width: int = 200, show_global_progress: bool = False, progress_interval: int = 200, init_msg: str | None = 'Sampling from target distribution...')

Configuration for the Metropolis–Hastings sampler.

class gpmp.mcmc.MetropolisHastings(log_target: Callable[[ndarray[tuple[Any, ...], dtype[floating]]], float], prop_rnd: Callable[[ndarray[tuple[Any, ...], dtype[floating]], int], ndarray[tuple[Any, ...], dtype[floating]]] = None, options: MHOptions = None)

Metropolis–Hastings sampler with two adaptive strategies:

• RM (Robbins–Monro) for diagonal or scalar proposals proposal_params[i] *= exp(step * (rate - target))

• Haario for full covariance: self.haario_scaling_factors[i] gets updated similarly, and each chain uses proposal covariance scaling * EmpCov + eps * I. The default 2.38**2 / dim factor comes from empirical results on multivariate Gaussians suggesting optimal acceptance near 0.23.

gpmp.mcmc.sample_multivariate_normal_with_jitter(mean, cov, initial_jitter=1e-08, max_attempts=5)

NUTS

class gpmp.mcmc.NUTSOptions(num_warmup: int = 1000, target_accept: float64 = 0.8, max_depth: int = 10, delta_max: float64 = 1000.0, jitter: float64 = 0.0001, init_step_size: float64 | None = None, init_mass_diag: any | None = None, seed: int | None = None, progress: bool = True, verbose: int = 1, log_every: int = 50, dual_averaging_gamma: float64 = 0.05, dual_averaging_t0: float64 = 10.0, dual_averaging_kappa: float64 = 0.75, dual_averaging_mu_factor: float64 = 10.0, warmup_min_no_window: int = 20, warmup_large_threshold: int = 150, warmup_large_init_buffer: int = 75, warmup_large_term_buffer: int = 50, warmup_large_base_window: int = 25, warmup_init_buffer_ratio: float64 = 0.15, warmup_term_buffer_ratio: float64 = 0.1, warmup_base_window_divisor: float64 = 3.0, find_eps_init: float64 = 1.0, find_eps_target_accept: float64 = 0.5, find_eps_scale_base: float64 = 2.0, find_eps_min: float64 = 1e-06, find_eps_max: float64 = 100.0)

Configuration object for NUTS sampling and warmup adaptation.

gpmp.mcmc.nuts_sample(log_prob: Callable[[any], any], q_init: any, num_samples: int, num_warmup: int = 1000, target_accept: float64 = 0.8, max_depth: int = 10, delta_max: float64 = 1000.0, jitter: float64 = 0.0001, init_step_size: float64 | None = None, init_mass_diag: any | None = None, seed: int | None = None, progress: bool = True, verbose: int = 1, log_every: int = 50, options: NUTSOptions | None = None) → Tuple[any, Dict[str, any]]

q_init: (chains, dim) log_prob: takes (dim,) and returns scalar

options:

Optional NUTSOptions object. Default keyword arguments override options when set to non-default values.

verbose:

0: silent 1: phase + window events + periodic summaries 2: more frequent summaries

gpmp.mcmc.nuts_transition(log_prob: Callable[[any], any], q0: any, step_size: float64, inv_mass_diag: any, max_depth: int, delta_max: float64) → Tuple[any, float64, int, int, bool]

gpmp.mcmc.plot_nuts_diagnostics(info: Dict[str, any], window: int = 50, show: bool = True, save_dir: str | None = None)

Sequential Monte Carlo

class gpmp.mcmc.ParticlesSetConfig(initial_distribution_type: str = 'randunif', resample_scheme: str = 'multinomial', param_s_initial_value: float = 0.5, param_s_upper_bound: float = 100000.0, param_s_lower_bound: float = 0.001, jitter_initial_value: float = 1e-16, jitter_max_iterations: int = 10, covariance_method: str = 'normal', covariance_knn_n_random: int = 20, covariance_knn_n_neighbors: int = 200)

class gpmp.mcmc.SMCConfig(compute_next_logpdf_param_method: str = 'p0', mh_steps: int = 20, mh_acceptation_rate_min: float = 0.15, mh_acceptation_rate_max: float = 0.3, mh_adjustment_factor: float = 1.4, mh_adjustment_max_iterations: int = 50)

gpmp.mcmc.run_smc_sampling(logpdf_parameterized_function, initial_logpdf_param: float, target_logpdf_param: float, compute_next_logpdf_param_method, min_ess_ratio: float, p0: float = None, init_box: list = None, n_particles: int = 1000, mh_steps: int = 20, smc_config: SMCConfig = None, particles_config: ParticlesSetConfig = None, debug: bool = False, plot_particles: bool = False, plot_empirical_distributions: bool = False)

Run a full Sequential Monte Carlo (SMC) simulation.

Parameters:

• logpdf_parameterized_function (callable) – Function of the form f(x, param) returning the log-density at x for a given param.

• initial_logpdf_param (float) – Initial value of the logpdf parameter.

• target_logpdf_param (float) – Final value of the logpdf parameter.

• compute_next_logpdf_param_method (str) – Method to compute the next logpdf parameter (‘p0’ or ‘ess’).

• min_ess_ratio (float) – Minimum acceptable ratio ESS / n. Below this threshold, a restart is triggered.

• p0 (float, optional) – Prescribed probability used if compute_next_logpdf_param_method is ‘p0’.

• init_box (array_like) – Domain box for particle initialization.

• n_particles (int, optional) – Number of particles. Default is 1000.

• mh_steps (int, optional) – Number of MH steps for moving the particles. Default is 20.

• smc_config (SMCConfig, optional) – An instance of SMCConfig to set SMC options.

• particles_config (ParticlesSetConfig, optional) – An instance of ParticlesSetConfig to set particle options. (Takes precedence on mh_steps setting.)

• debug (bool, optional) – If True, print debug information during execution.

• plot_particles (bool, optional) – If True, display a matrix plot of the final particle distribution.

• plot_empirical_distributions (bool, optional) – If True, display the final particle distribution.

Returns:

• particles (ndarray) – Final particle positions.

• smc (SMC) – The SMC instance containing diagnostics and logs.

gpmp.mcmc.log_indicator_density(f, threshold, log_px, tail='lower')

Return logpdf(x) = log(1_{f(x) ? threshold} * p_X(x)) where ? depends on tail.

gpmp.mcmc.run_subset_simulation(f, thresholds, init_box, log_px, tail='upper', n_particles=1000, mh_steps=20, min_acceptation=0.15, max_acceptation=0.3, resample_scheme='residual', debug=False)

Estimate P(f(X) ? u_T), with ? = < or >, via Subset Simulation.

Parameters:

• f (callable) – Function from R^d to R (performance or score function).

• thresholds (list of float) – Monotonic threshold sequence: decreasing for ‘<’, increasing for ‘>’.

• init_box (list of [lower_bounds, upper_bounds]) – Sampling domain for initial distribution.

• log_px (callable) – Log-density of the base distribution p_X.

• tail (str) – Either ‘lower’ (f < u_i) or ‘upper’ (f > u_i).

Returns:

• p_estimate (float) – Final estimate of P(f(X) ? u_T).

• stage_probs (list of float) – Estimated conditional probabilities p_{u_i | u_{i-1}}.

• smc (SMC) – The SMC object with diagnostics.