VBMC

Note

The VBMC class implements the Variational Bayesian Monte Carlo (VBMC) algorithm.

VBMC computes a variational approximation of the full posterior and a lower bound on the log normalization constant (log marginal likelihood or log model evidence) for a provided unnormalized log posterior.

To perform inference, first initialize a VBMC object and then call vbmc.optimize() on the instance.

By default VBMC assumes noiseless evaluations of the log posterior, but noisy likelihoods can also be handled. See PyVBMC Example 6: Noisy log-likelihood evaluations for more details.

See below for more details on the VBMC class methods and interface. The primary entry-points for users are the VBMC class, which initializes the algorithm, and the VariationalPosterior class, which represents the returned variational solution. The Basic options may also be useful.

The iteration history (vbmc.iteration_history) keeps each iteration’s Gaussian process surrogate as its training data and hyperparameters, without the posterior factors; vbmc.get_gp(iteration) returns the surrogate of an iteration with its posterior restored, ready for prediction.

Reusing evaluations and charging initialization work

Use precomputed_evaluations to give VBMC target evaluations that are already available when the run starts. These observations are independent of the starting point and plausible box: they neither become additional rows of x0 nor change the plausible bounds. Pass (X, y) for an exact target or one with unknown noise, where X has shape (N, D) in the target’s original input coordinates and y has shape (N,):

vbmc = VBMC(
    log_density,
    x0,
    lower_bounds,
    upper_bounds,
    plausible_lower_bounds,
    plausible_upper_bounds,
    precomputed_evaluations=(X, y),
)

The values in y must be outputs of log_density at the corresponding rows of X. When prior or log_prior is supplied separately, log_density and y contain log-likelihood values; VBMC adds the deterministic log-prior once. Precomputed points must lie strictly inside the hard bounds, but they need not lie inside the plausible box.

For a target with user-provided noise estimates, set options["specify_target_noise"] and pass (X, y, y_sd). y_sd has shape (N,) and contains positive standard deviations of the log-density observations, not uncertainties in X. Repeated rows of X remain independent noisy observations and are pooled by inverse-variance weighting:

vbmc = VBMC(
    noisy_log_density,
    x0,
    lower_bounds,
    upper_bounds,
    plausible_lower_bounds,
    plausible_upper_bounds,
    precomputed_evaluations=(X, y, y_sd),
    options={"specify_target_noise": True},
)

Precomputed observations do not count as fresh target calls. By default they also carry no cost for the current run, which is appropriate for evaluations obtained during earlier work. If preparing this run incurred target or derivative work, pass its nonnegative integer function-equivalent charge as initialization_cost. The charge is independent of the number of retained observations and counts once against the total max_fun_evals allowance:

vbmc = VBMC(
    log_density,
    x0,
    lower_bounds,
    upper_bounds,
    plausible_lower_bounds,
    plausible_upper_bounds,
    precomputed_evaluations=(X, y),
    initialization_cost=13,
    options={"max_fun_evals": 100},
)

Here VBMC can make at most 87 fresh target calls. results["func_count"] remains the literal number of fresh calls made by VBMC. When initialization_cost is positive, results["evaluation_budget"] additionally reports the function-equivalent limit, initialization, new_calls and used total. Zero-cost runs retain the ordinary result and display without this additional budget report.

Performance calibration

New runs reuse compatible machine-local chunk settings produced by an explicit calibrate campaign. If no compatible record exists, they use historical defaults. Optimization never starts a campaign. Use the performance_calibration option to select "cached" (default), "off", or an explicit CalibrationProfile.

Settings are fixed before the first relevant kernel call and survive final boost and save/resume. A resolved save rejects load-time "cached" or different explicit settings; start a new run to use a different profile.

Startup tips

With show_tips=True (default), a new run may print one tip before the iteration display. Resumed and continued runs do not show tips. PyVBMC may use the same space for a performance-calibration suggestion or an old-release reminder.

Set show_tips=False to disable tips and the old-release reminder. Set display="off" to suppress tips and both reminders. Pass either setting in the options dictionary when constructing VBMC. The FAQ also describes the old-release reminder.

class pyvbmc.VBMC(log_density: callable, x0: ndarray = None, lower_bounds: ndarray = None, upper_bounds: ndarray = None, plausible_lower_bounds: ndarray = None, plausible_upper_bounds: ndarray = None, options: dict = None, options_path: PathLike | None = None, prior: Prior = None, log_prior: callable = None, sample_prior: callable = None, seed=None, *, precomputed_evaluations=None, initialization_cost=0)

Posterior and model inference via Variational Bayesian Monte Carlo (VBMC).

VBMC computes a variational approximation of the full posterior and a lower bound on the normalization constant (marginal likelihood or model evidence) for a provided unnormalized log posterior.

Initialize a VBMC object to set up the inference problem, then run optimize(). See the examples for more details.

Parameters:

log_densitycallable or pyvbmc.pymc.PyMCTarget

A given target log-posterior or log-likelihood, or a constructed pyvbmc.pymc.PyMCTarget. A PyMCTarget supplies its target function, starting point, bounds, setup evaluations, and setup cost; explicitly supplied constructor arguments override its starting point, plausible bounds, evaluations, and cost. Supplied hard bounds must equal target.lb and target.ub, which are in the target’s flat coordinates: the support of a variable the adapter leaves untransformed, and infinite for a variable that keeps its PyMC transform. An additional prior or log_prior cannot be combined with a PyMCTarget. For an ordinary callable, if log_prior is None, log_density accepts input x and returns the value of the target log-joint, that is, the unnormalized log-posterior density at x. If log_prior is not None, log_density should return the unnormalized log-likelihood. In either case, if options["specify_target_noise"] is true, log_density should return a tuple where the first element is the noisy log-density, and the second is an estimate of the standard deviation of the noise. If options["vectorized_target"] is true, it instead receives a NumPy array with shape (N, D) and returns values with shape (N,) or (N, 1). A noisy vectorized target returns a pair of arrays with one value and positive standard deviation per row. The initial design is evaluated in one batch; subsequent calls contain one point, with shape (1, D). A separately supplied prior is still evaluated one point at a time.

x0array_like, optional

Starting point for the inference. Ideally x0 is a point in the proximity of the mode of the posterior. One starting point is a row of D coordinates (a 1-D array or a list; for D = 1, also a number), and several are the rows of an array of shape (n0, D) or a list of lists. Default is None, in which case the number of variables is read from the plausible bounds, one of which then needs one entry per variable.

lower_bounds, upper_boundsnp.ndarray, optional

lower_bounds (LB) and upper_bounds (UB) define a set of strict lower and upper bounds for the coordinate vector, x, so that the posterior has support on LB < x < UB. If scalars, the bound is replicated in each dimension. Use None for LB and UB if no bounds exist. Set LB [i] = -inf and UB [i] = inf if the i-th coordinate is unbounded (while other coordinates may be bounded). Note that if LB and UB contain unbounded variables, the respective values of PLB and PUB need to be specified (see below), by default None.

plausible_lower_bounds, plausible_upper_boundsnp.ndarray, optional

Specifies a set of plausible_lower_bounds (PLB) and plausible_upper_bounds (PUB) such that LB < PLB < PUB < UB. Both PLB and PUB need to be finite. PLB and PUB represent a “plausible” range, which should denote a region of high posterior probability mass. If scalars, the bound is replicated in each dimension. Among other things, the plausible box is used to draw initial samples and to set priors over hyperparameters of the algorithm. When in doubt, we found that setting PLB and PUB using the topmost ~68% percentile range of the prior (e.g, mean +/- 1 SD for a Gaussian prior) works well in many cases (but note that additional information might afford a better guess). Both are by default None.

optionsdict, optional

Additional options can be passed as a dict. Please refer to the VBMC options page for the default options. If no options are passed, the default options are used.

options_pathpathlike, optional

Additional options can also be specified by a .ini file. The options dict keyword takes precedence over options set in the .ini file. options_path should be either absolute or relative to the pyvbmc/vbmc/ directory.

prior, optional

An optional separate prior. It can be a pyvbmc.priors.Prior subclass, an appropriate scipy.stats distribution, or a list of one-dimensional PyVBMC Prior and/or one-dimensional continuous scipy.stats distributions. (see the documentation on priors for more details). If prior is not None, the argument log_density is assumed to represent the log-likelihood (otherwise it is assumed to represent the log-joint). The hard bounds lower_bounds and upper_bounds must lie inside the support of the prior. Within them the prior is used as it is given: the model evidence reported is that of the prior restricted to the hard bounds, not that of a prior truncated to them and normalized again.

log_priorcallable, optional

An optional separate log-prior function, which should accept a single argument x and return the log-density of the prior at x. If log_prior is not None, the argument log_density is assumed to represent the log-likelihood (otherwise it is assumed to represent the log-joint).

sample_priorcallable, optional

An optional function which accepts a single argument n and returns an array of samples from the prior, of shape (n, D), where D is the problem dimension. Currently unused.

seedNone, int, SeedSequence or np.random.Generator, optional

Seed (or generator) for the random generator vbmc.rng used by this instance, by its variational posterior and by the GP hyperparameter fit, for reproducible runs. Anything accepted by numpy.random.default_rng. If None (default), a generator is derived from NumPy’s global random state, so that calling np.random.seed beforehand still makes a run reproducible; deriving it draws four integers from that state and so advances it, which a construction with a seed or a generator does not do. Every random draw of a run comes from vbmc.rng, and after the construction NumPy’s global random state is neither read nor written. Draws from vbmc.rng before optimize (e.g. vbmc.vp.sample) change the run like any other change of the generator’s state. This seed does not control random draws made by log_density; seed a stochastic target separately.

precomputed_evaluationstuple, optional

Evaluations available before the run, independently of x0. Pass (X, y) for exact or unknown-noise targets, or (X, y, y_sd) when options["specify_target_noise"] is true. X has shape (N, D) in the target’s original coordinates and the other arrays have shape (N,). If a separate prior is supplied, y contains log-likelihood values; VBMC adds the prior once. (The values of options["f_vals"], given at x0, are log-joint values, with the prior already added.) Inputs are copied to float64 storage. Exact duplicate points must agree within four float64 ULPs at their value scale and are retained once; noisy repeats remain independent observations and are pooled by the logger.

initialization_costint, optional

Nonnegative function-equivalent cost incurred while preparing this run. It is independent of the number of precomputed rows and counts once against options["max_fun_evals"], leaving the difference for fresh target calls. Default is zero.

Attributes:

rngnp.random.Generator

The random generator used by this instance. It is shared with vbmc.vp and with the variational posteriors returned by optimize.

targetpyvbmc.pymc.PyMCTarget or None

PyMC target supplied directly at construction, if any.

x0_orignp.ndarray, shape (n0, D)

The starting points in the coordinates the caller gave them in. x0 holds the same points in the inference space the instance was constructed with, which a later warp of that space leaves behind.

Raises:

ValueError

When neither x0 or (plausible_lower_bounds and plausible_upper_bounds) are specified.

ValueError

When various checks for the bounds (LB, UB, PLB, PUB) of VBMC fail.

ValueError

When an option value, precomputed_evaluations or initialization_cost is refused, when initialization_cost leaves too few evaluations of max_fun_evals for the initial design, when the hard bounds reach outside the support of prior, or when at least fun_eval_start starting points are given and the first fun_eval_start of them take one value in a coordinate.

NotImplementedError

When the options select a feature of MATLAB VBMC that is not ported (noise_shaping, acq_hedge, a gp_hyp_sampler other than "slicesample", an acquisition function that asks for the MCMC step of the importance sampler).

References

[1]

Huggins, B., Li, C., Tobaben, M., Aarnos, M., & Acerbi, L. (2023). “PyVBMC: Efficient Bayesian inference in Python”. Journal of Open Source Software 8(86), 5428, https://doi.org/10.21105/joss.05428.

[2]

Acerbi, L. (2018). “Variational Bayesian Monte Carlo”. In Advances in Neural Information Processing Systems 31 (NeurIPS 2018), pp. 8213-8223.

[3]

Acerbi, L. (2020). “Variational Bayesian Monte Carlo with Noisy Likelihoods”. In Advances in Neural Information Processing Systems 33 (NeurIPS 2020).

Examples

For VBMC usage examples, please look up the Jupyter notebook tutorials in the PyVBMC documentation: https://acerbilab.org/pyvbmc/_examples/pyvbmc_example_1_basic_usage.html

Previously computed target values can seed the logger without becoming starting points or changing the plausible box. Preparation performed for the current run can be charged against the total function-equivalent budget independently of the number of retained observations:

cached = (X_evaluated, log_density_values)
vbmc = VBMC(
    log_density,
    x0,
    lower_bounds,
    upper_bounds,
    plausible_lower_bounds,
    plausible_upper_bounds,
    precomputed_evaluations=cached,
    initialization_cost=setup_cost,
    options={"max_fun_evals": 100},
)

determine_best_vp(max_idx: int = None, safe_sd: float = None, frac_back: float = None, rank_criterion_flag: bool = None)

Return the best VariationalPosterior found during the optimization of VBMC as well as its ELBO, ELBO_SD and the index of the iteration.

Called without the selection arguments, the method uses the run’s options, so that on a finished run it returns the iteration the run selected.

Parameters:

max_idxint, optional

Check up to this iteration, by default None which means last iter.

safe_sdfloat, optional

Penalization for uncertainty: the multiple of the ELBO_SD subtracted from the ELBO to give the ELCBO. By default None, which means the option best_safe_sd (5 unless set).

frac_backfloat, optional

If no past stable iteration, go back up to this fraction of iterations. Used only without the ranking criterion. By default None, which means the option best_frac_back (0.25 unless set).

rank_criterion_flagbool, optional

If True use new ranking criterion method to pick best solution. It finds a solution that combines ELCBO, stability, and recency. By default None, which means the option rank_criterion (True unless set).

Returns:

vpVariationalPosterior

A copy of the VariationalPosterior of the selected iteration, carrying that iteration’s stability flag in its statistics. It shares the random generator with the recorded posterior, as every copy of a VariationalPosterior does, and the recorded posterior itself is left untouched.

elbofloat

The ELBO of the iteration with the best VariationalPosterior.

elbo_sdfloat

The ELBO_SD of the iteration with the best VariationalPosterior.

idx_bestint

The index of the iteration with the best VariationalPosterior.

Notes

A NaN ELCBO (elbo - safe_sd * elbo_sd) or reliability index ranks last in the ranking criterion, and a NaN ELCBO is skipped without it. If the ELCBO of every candidate iteration is NaN, the last iteration considered is selected, which with the default max_idx is the posterior the run ended on.

final_boost(vp: VariationalPosterior, gp: GP)

Perform a final boost of variational components.

Parameters:

vpVariationalPosterior

The VariationalPosterior that should be boosted.

gpgpyreg.GP

The Gaussian process the VariationalPosterior was fitted with.

Returns:

vpVariationalPosterior

The VariationalPosterior resulting from the final boost.

elbofloat

The ELBO of the VariationalPosterior resulting from the final boost.

elbo_sdfloat

The ELBO_SD of the VariationalPosterior resulting from the final boost.

changed_flagbool

Whether the returned VariationalPosterior is the boosted one: False when no boost was needed and when the guard rejected the boosted posterior.

Raises:

ValueError

If the option tol_elcbo_boost is neither None nor a finite nonnegative number; or, with variable_means off, when gp has fewer training inputs than vp has components.

RuntimeError

If, with the guard on, neither vp nor the boosted posterior has a finite ELBO and a finite nonnegative ELBO SD.

Notes

A boost is needed when vp has fewer components than the boost gives it (min_final_components with variable_means on, the number of training inputs of gp with it off), or when ns_ent_boost or ns_ent_fine_boost changes the corresponding entropy sample size.

With the option tol_elcbo_boost set (0.1 by default), a guard keeps the boosted posterior only if neither its ELBO nor its ELBO minus five times its ELBO SD falls more than tol_elcbo_boost below that of vp; otherwise it warns and returns vp. A boosted posterior whose ELBO or ELBO SD is not finite, or whose SD is negative, is rejected in the same way, and one whose statistics are valid is kept when those of vp are not. With tol_elcbo_boost=None the boosted posterior is always returned. The guard compares the optimizer’s stored pre- and post-boost ELBO and GP-based ELBO SD. It performs no diagnostic rescoring, and the SD does not include Monte Carlo uncertainty from the entropy estimate.

The boost optimizes with warm-up over, without pruning of components, with the entropy sample sizes of the boost options (ns_ent_boost and its kin), with the stochastic optimization allowed its ceiling of 10000 iterations whatever max_iter_stochastic says and, under the guard, without the weight penalty (weight_penalty). It runs on copies of the options and of the optimization state; neither the instance’s options nor its optimization state is changed.

With variable_means off the components of the boosted posterior sit at the training inputs of gp, one each, and min_final_components does not apply.

get_gp(iteration: int)

Return the Gaussian process surrogate of a recorded iteration.

The iteration history stores each iteration’s GP without its posterior factors: the training data, the hyperparameter samples and the model (covariance, mean and noise functions, hyperparameter priors and bounds), but not the Cholesky factors and the other per-sample quantities that a GP prediction needs, which would make the history grow with the square of the training-set size. This method returns a copy of that record with the factors recomputed from the stored data and hyperparameters. They are a deterministic function of those, so the returned GP is the one the algorithm held at that iteration. The stored record is left without them, and a record that does carry its factors (an instance saved by a version of PyVBMC that stored complete GPs) is copied as it is.

Parameters:

iterationint

The iteration whose GP to return, from 0 to the last recorded iteration.

Returns:

gpgpyreg.GP

The Gaussian process of that iteration, ready for prediction.

Raises:

ValueError

If no GP has been recorded (the instance has not been optimized), or if iteration is not the index of a recorded iteration.

classmethod load(file, new_options=None, iteration=None, set_random_state=False)

Load a VBMC instance from a file.

Note

A saved instance holds Python bytecode: always that of the options whose value is a function (ns_ent, k_fun_max and the like), and that of the target function, and of the function that PyVBMC builds from the target and a separate prior, whenever they cannot be pickled by name (a lambda, a function defined in a script or inside another function). Bytecode belongs to the minor version of Python that wrote the file. Under another minor version the file can be loaded and inspected (the variational posterior, the iteration history and the evaluated points are data), but the run should not be continued or saved again there: both reach that bytecode, which can end the interpreter. To move a result between Python versions, save the variational posterior on its own (vp.save); its file holds no bytecode.

Loading a run made with VBMC(target) from a PyMCTarget checks the arrays that the target could not copy and evaluates its log density once, at the starting point, against the value recorded during setup. If the target no longer matches the run, loading warns; the run’s results remain valid, and optimize refuses to continue the run.

Parameters:

filepath-like

The file name or path to read from. Default file extension .pkl will be added if no extension is specified.

new_optionsdict or None

Options to apply to the loaded run. This is useful for increasing max_fun_evals or max_iter before continuing. Names and values are checked as they are during construction. Options used only to construct a run, such as integer_vars, gp_mean_fun and warmup, cannot be changed on load; start a new run to change them. Supplying their stored value is allowed. Changing ns_gp_max updates GP hyperparameter sampling unless sampling has already stopped in the stable regime. Options with no effect, and noise_size when the target supplies its own noise estimate, produce the same warnings as at construction. See the changelog for compatibility details about files saved by earlier releases.

iterationint or None

The iteration at which to initialize the stored VBMC instance. Default is None, meaning initialize to the last recorded iteration.

set_random_statebool

Whether to restore the state of the instance’s generator (vbmc.rng) stored at the specified iteration, for reproducibility. Default False, in which case vbmc.rng continues from the state it had when the file was saved. Files saved before vbmc.rng existed only stored NumPy’s global random state; loading them with set_random_state=True restores it, derives a fresh generator from it (leaving the global state where it was saved), and warns. Files saved while the GP hyperparameter fit still drew from the global state also hold that state; it is ignored, so a run resumed from such a file does not replay the original exactly.

Returns:

vbmcVBMC

The loaded VBMC instance, initialized at the specified iteration.

Raises:

ValueError

If the specified iteration is less than zero or larger than the last stored iteration; if an option name or value is invalid; or if new_options tries to change an option that is fixed at construction.

NotImplementedError

If the options select a feature of MATLAB VBMC that is not ported (noise_shaping, acq_hedge, a gp_hyp_sampler other than "slicesample", an acquisition function that asks for the MCMC step of the importance sampler).

OSError

If the file cannot be found, or cannot be opened for other reasons.

optimize()

Run inference on an initialized VBMC object.

VBMC computes a variational approximation of the full posterior and the ELBO (evidence lower bound), a lower bound on the log normalization constant (log marginal likelihood or log model evidence) for the provided unnormalized log posterior.

Returns:

vpVariationalPosterior

The VariationalPosterior computed by VBMC.

resultsdict

A dictionary with the following information about the run:

• "success_flag" (bool): whether the returned posterior was marked stable.

• "message" (str): reason the optimization stopped.

• "elbo" (float): estimated evidence lower bound for the returned posterior.

• "elbo_sd" (float): estimated standard deviation of the ELBO.

• "iterations" (int): number of iterations performed.

• "func_count" (int): number of fresh target evaluations, counting each row of a vectorized call.

• "best_iter" (int): zero-based index of the iteration whose posterior was selected for return, before any final boost.

• "train_set_size" (int): number of observations represented in that iteration’s training set; pooled repeated observations are counted separately.

• "components" (int): number of mixture components in the returned posterior.

• "r_index" (float): reliability index at best_iter.

• "convergence_status" (str): "probable" if the selected iteration was stable, otherwise "no".

• "problem_type" (str): "bounded" if any hard bound is finite, otherwise "unconstrained".

• "function" (str): string representation of the function logged by VBMC.

• "algorithm" (str): algorithm name.

• "performance_calibration" (dict): serialized calibration profile used by the run.

• "version" (str or None): installed PyVBMC version, or None if package metadata is unavailable.

• "rng_state" (dict): state of vbmc.rng after the run.

• "overhead" (float): reserved field, currently np.nan.

When at least one precomputed observation was supplied, the dictionary also contains "precomputed_observations" (int), the number supplied, and "precomputed_locations" (int), the number of distinct input points among them.

When initialization_cost is positive, it also contains "evaluation_budget" (dict), with these entries:

• "unit" (str): "function_equivalent_evaluations".

• "limit" (int or float): the total max_fun_evals allowance, possibly np.inf.

• "initialization" (int): the initialization charge.

• "new_calls" (int): fresh target evaluations, equal to func_count.

• "used" (int): initialization charge plus fresh target evaluations.

Raises:

ValueError

If no fresh function-equivalent budget remains for the run.

RuntimeError

If the PyMCTarget of a run made with VBMC(target) no longer matches the run: an array that the model reads inside the body of a pytensor.scan loop or of an OpFromGraph was changed in place, or, for a loaded run, the target no longer computed the log density recorded during setup. Restoring the array lets a run that was not loaded continue; otherwise, build a new target and run to fit the changed data.

Notes

Every inference draw of the run comes from vbmc.rng (see the seed parameter); NumPy’s global random state is neither read nor written. Optional startup tips use a separate private random stream.

save(file, overwrite=False)

Save the VBMC instance to a file.

Note

A saved instance holds Python bytecode: always that of the options whose value is a function (ns_ent, k_fun_max and the like), and that of the target function, and of the function that PyVBMC builds from the target and a separate prior, whenever they cannot be pickled by name (a lambda, a function defined in a script or inside another function). Bytecode belongs to the minor version of Python that wrote the file. Under another minor version the file can be loaded and inspected (the variational posterior, the iteration history and the evaluated points are data), but the run should not be continued or saved again there: both reach that bytecode, which can end the interpreter. To move a result between Python versions, save the variational posterior on its own (vp.save); its file holds no bytecode.

Saving a run made with VBMC(target) whose PyMCTarget no longer matches the run (an array that the model reads inside the body of a pytensor.scan loop or of an OpFromGraph was changed in place, or loading marked the run) warns: the file records the target as it is, and optimize refuses to continue the run once it is loaded.

Parameters:

filepath-like

The file name or path to write to. Default file extension .pkl will be added if no extension is specified.

overwritebool

Whether to allow overwriting existing files. Default False.

Raises:

FileExistsError

If the file already exists and overwrite is False.

OSError

If the file cannot be opened for other reasons (e.g., the directory is not found, the disk is full, etc.).

property target

PyMC target supplied directly at construction, if any.

The lower-level functions used by the algorithm are listed under VBMC algorithm helpers.