`SVBMC`#

Note

The SVBMC class implements Stacking Variational Bayesian Monte Carlo (S-VBMC): it combines the variational posteriors of several completed VBMC runs on the same problem into one stacked posterior by optimizing the stacked ELBO with respect to the mixture weights only. The runs’ components and their stored expected log-joint statistics stay fixed, so stacking needs no further evaluations of the model. Use it when repeated VBMC runs on the same model and data disagree, or each captures part of a multimodal posterior.

S-VBMC needs the optional torch extra (see Installation); constructing an SVBMC object without it raises an ImportError. Importing PyVBMC does not import torch.

Basic usage#

Run VBMC several times on the same target, collect the returned posteriors, and stack them:

from pyvbmc import VBMC, SVBMC

vps = []
for seed in range(5):
    vbmc = VBMC(log_joint, x0[seed], lb, ub, plb, pub, seed=seed)
    vp, results = vbmc.optimize()
    vps.append(vp)

stacked = SVBMC(vps, seed=0)
stacked.optimize()

samples = stacked.sample(10000)          # (10000, D), original space
print(stacked.elbo)                       # headline stacked ELBO
print(stacked.elbo_sd)                    # uncertainty of the raw estimate
print(stacked.elbo_details["raw"])
print(stacked.elbo_details["shrunk_two_level"])

PyVBMC Example 7: Stacking the posteriors of several runs (S-VBMC) walks through this on a bimodal target.

ELBO reporting#

stacked.elbo is the headline estimate of the stacked ELBO, chosen by the noise status of the stack:

For a stack treated as noiseless, it is raw, the ELBO from the fresh final evaluation at the selected weights.
If any retained run is treated as noisy, it is shrunk_two_level, the two-level shrinkage estimate at the same weights and with the same final entropy, whatever the version of the optimization. Where that calculation is numerically undefined, the headline is capped_I_median instead, and optimize() warns.

The choice affects the reported value alone: the weights, and so the stacked posterior and its samples, come from the same optimization whatever the headline.

The weights maximize the weighted sum of the components’ expected log-joints, each estimated by its run’s GP, plus the entropy of the mixture. On a noisy target, this optimization favors components whose estimates came out high and adds an upward bias to the raw ELBO. Here bias is measured against the stack’s ELBO computed with the exact log joint, which is itself a lower bound on the log evidence.

The two-level shrinkage estimate targets this selection bias. It first shrinks component expected log-joints toward the mean within each run, using their GP uncertainty, and then shrinks the runs’ levels toward their common mean. It can also remove some bias inherited from the input VBMC runs, but it does not make the headline unbiased: miscalibrated GP uncertainty and different run-level biases can leave an error in either direction. The raw ELBO can likewise fall below the stack’s exact ELBO when the input runs are pessimistic.

On our benchmarks with every component weight optimized, the remaining stacking optimism and the inherited bias removed by shrinkage roughly balanced, so noisy-stack headlines were about as biased as the input runs. With version="ns", stacking adds no weight-selection bias; shrinkage still applies and can move the headline below the raw value.

stacked.elbo_details contains the complete report:

Key	Meaning
`raw`	ELBO from the fresh evaluation at the selected weights; the headline of a noiseless stack.
`capped_I_median`	ELBO with the expected log-joint capped at the component median; the headline of a noisy stack whose shrinkage is unavailable.
`capped_E_median`	ELBO with the expected log-joint capped at the retained-run median.
`naive`	Baseline that gives every retained run equal total weight and keeps that run’s original internal component weights. It inherits errors in the input runs and is not a bound on the optimized ELBO.
`shrunk_two_level`	Empirical-Bayes estimate that first shrinks component expected log-joints within each retained run using their full GP covariance, then shifts each run’s components by the change that shrinkage makes to the run’s level; the headline of a noisy stack. It uses the selected stacking weights and the same final entropy as `raw`. `None` means finite input statistics led to a numerically undefined calculation; `optimize()` emits a `RuntimeWarning` and leaves every other report value available.
`shrinkage_noise_share`	Diagnostic ratio of estimated noise to component spread, averaged over runs using their selected mixture masses. It can exceed one and is `None` when shrinkage is unavailable or no run with positive spread has positive selected mass.
`headline_method`	The key of the headline: `"raw"` for a noiseless stack, `"shrunk_two_level"` for a noisy stack, or `"capped_I_median"` for a noisy stack whose shrinkage is unavailable.
`cap_amount`	Reduction that the median cap applied to the headline, `raw` less `capped_I_median`. It is zero unless the headline is `capped_I_median`, even when a diagnostic cap would bind.
`entropy_sd`	Monte Carlo contribution to the uncertainty of `raw`.
`gp_sd`	GP quadrature contribution, treating retained runs as independent.
`raw_sd`	Combined raw-estimate uncertainty, equal to `stacked.elbo_sd`.
`noisy`	Whether the stack is treated as noisy.
`noise_status_source`	Per-retained-run provenance: `"recorded"`, `"inferred"` or `"override"`.

stacked.elbo_sd combines the stratified entropy estimator’s Monte Carlo variance with the GP quadrature uncertainty at the selected weights. It describes the raw value. It excludes selection bias and does not propagate the shrinkage or the median cap, so it must not be interpreted as defining a calibrated confidence interval for the headline of a noisy stack.

A composite posterior#

Each VBMC run works in its own transformed parameter space (bounded parameters, and possibly a rotation and rescaling learned during the run). The components of different runs therefore live in incompatible spaces and are only comparable after mapping back to the original space. The stacked posterior keeps the input posteriors as they are and draws from each through its own transform. Do not combine the component means and covariances of different runs directly, and do not read them as a single mixture: use sample() for every estimate and plot.

Helpers#

pyvbmc.svbmc.utils holds two helpers from the S-VBMC workflow: overlay_corner_plot draws several sample sets on one corner plot, to compare the individual runs with the stacked posterior, and find_init_bounds returns the box to draw VBMC starting points from given the bounds and plausible bounds of a problem.

Differences from the standalone `svbmc` package#

The implementation derives from the standalone S-VBMC package (acerbilab/svbmc, version 0.1.1), which is deprecated in its favor: use pyvbmc.SVBMC for new and existing code. The package’s BSD 3-Clause license notice ships with the subpackage, pyvbmc.svbmc. The method and the optimizer of the weights are the same, and on our benchmark the two reach the same weights. What differs:

It is faster. The work that the components of one run share, its parameter transform above all, is done once per run rather than once per component.
The reported ELBO is more precise. The corrections that bring the expected log-joint of each run to the original parameter space are computed once, when the object is created, exactly or by numerical quadrature, where the standalone package estimated them by Monte Carlo at every step. After the optimization the ELBO is evaluated again with more draws (optimize(n_samples_final=...)), and it comes with an uncertainty, elbo_sd.
elbo is a number, the headline described under ELBO reporting, which for a noisy stack is a shrinkage estimate that the standalone package does not compute. The entries of the standalone package’s elbo dictionary are in elbo_details, under the names in the table below.
seed takes the place of testing: every random draw comes from the generator of the SVBMC object, and the input posteriors are neither modified nor advanced.
sample(n) returns exactly n draws. The standalone package rounded the share of each run separately and could return a few more or fewer.
Results are float64, problems with one variable work, and malformed inputs raise errors.

Code written for the standalone package’s dictionary-valued elbo can be migrated as follows:

Standalone `svbmc`	`pyvbmc.svbmc`
`stacked.elbo["estimated"]`	`stacked.elbo_details["raw"]`
`stacked.elbo["debiased_I_median"]`	`stacked.elbo_details["capped_I_median"]`
`stacked.elbo["debiased_E_median"]`	`stacked.elbo_details["capped_E_median"]`

Citation#

Silvestrin, F., Li, C., & Acerbi, L. (2025). Stacking Variational Bayesian Monte Carlo. Transactions on Machine Learning Research. arXiv:2504.05004, TMLR. Please cite it together with the VBMC and PyVBMC papers when you use S-VBMC.

class pyvbmc.svbmc.SVBMC(vp_list: list, s_max: float = np.float64(2.23606797749979), M_min: float | int = 0.6666666666666666, seed=None, *, show_tips: bool = True, noisy: bool | None = None)[source]#

Stacking Variational Bayesian Monte Carlo (S-VBMC).

Combines the variational posteriors of several completed VBMC runs on the same problem into one stacked mixture by optimizing the stacked ELBO with respect to the mixture weights only. Component locations and scales, and the expected log-joint statistics VBMC stored with each posterior, stay fixed, so no further evaluations of the model are needed.

Construct the object with the fitted posteriors, then call optimize(). Construction discards runs that VBMC did not mark as converged (stats["stable"]) and runs whose expected log-joint carries excessive uncertainty (s_max), and raises when fewer than M_min runs remain.

The stacked posterior is a composite over the input posteriors: each VBMC run has its own parameter transform, so components of different runs live in different transformed spaces and are only comparable in the original space. Use sample() for estimates and plots; do not combine the component means and covariances of different runs directly.

Requires the optional torch extra (pip install --upgrade "pyvbmc[torch]"); construction raises ImportError without it. An object saved with save() can be loaded, sampled and plotted without it.

Parameters:

vp_listlist of VariationalPosterior: Posteriors returned by completed VBMC runs on the same model, data and original parameter space. Construction verifies the dimension, the mixture shapes, the statistics stacking reads, and that every run carries the same original-space hard bounds; the runs are free to map that space differently. That they were run on the same model and data does not follow from the posteriors and is not verified.
s_maxfloat, optional: Tolerance on the standard deviation of the expected log-joint of any single component, sqrt(max(stats["J_sjk"])). A run at or above it is discarded as poorly converged. Default sqrt(5).
M_minint or float, optional: Minimum number of runs that must survive the filters: a count when greater than 1, a proportion of len(vp_list) when at most 1 (so exactly 1 means every run). Default 2/3.
seedNone, int, array_like[int], SeedSequence, BitGenerator or Generator, optional: Seed of the object’s random generator, which drives the Monte Carlo entropy estimate and sample(). None derives the generator from NumPy’s global state, so np.random.seed before construction still fixes a run. Default None.
show_tipsbool, optional: Show occasional guidance when INFO logging is enabled. Set to False to suppress tips while keeping the messages on progress and on the headline ELBO. Default True.
noisybool or None, optional: Override the noise status of the whole stack. With None (the default), read each retained run’s uncertainty_handling_level statistic; when absent, infer noise from elbo_sd > 0.1. This inference can misclassify old posteriors. Any noisy retained run makes the headline elbo the two-level shrinkage estimate, where it is defined.

Attributes:

vp_listlist of VariationalPosterior: The retained posteriors, in input order. They are read, never modified: sampling works on copies that use this object’s generator.
Dint: Dimension of the problem.
Mint: Number of retained runs.
Klist of int: Number of components of each retained run.
wnp.ndarray, shape (1, K_total): Weights of all components, concatenated in run order and normalized. optimize() replaces them with the optimized ones.
elbofloat or None: After optimize(): the headline estimate of the stacked ELBO at the returned weights, evaluated with fresh draws. For a noiseless stack it is the raw ELBO; for a noisy stack, the two-level shrinkage estimate elbo_details["shrunk_two_level"], or the component-median capped ELBO where the shrinkage is numerically undefined. Choosing the weights on noisy estimates adds an upward bias to the raw ELBO of a noisy stack, which the shrinkage targets. The shrinkage can also remove part of the bias that each run’s ELBO carries from VBMC, and residual bias can remain. The raw ELBO is not an upper bound on the stack’s true ELBO.
elbo_sdfloat or None: Estimated uncertainty of the raw evaluation at the returned weights: entropy Monte Carlo error and GP quadrature uncertainty, treating runs as independent. It excludes selection bias and the uncertainty of the shrinkage or of the cap; it does not define a calibrated confidence interval for the headline of a noisy stack.
elbo_detailsdict or None: After optimize(): raw, capped_I_median, capped_E_median, naive and shrunk_two_level ELBO estimates; shrinkage_noise_share; headline_method (raw, shrunk_two_level or capped_I_median); cap_amount (the reduction the median cap applied to the headline, zero unless the headline is capped_I_median); entropy_sd, gp_sd and raw_sd (equal to elbo_sd); noisy and noise_status_source. Shrinkage adjusts component estimates within each run and then shifts run levels, at the selected weights and with the same final entropy as raw. shrunk_two_level is None when its numerical calculation is undefined, in which case optimization still succeeds with a RuntimeWarning. shrinkage_noise_share is the estimated noise-to-spread ratio averaged over runs using their selected masses; it can exceed one and can be None independently. Naive stacking equally weights retained runs with their original internal weights. Its estimate inherits the runs’ errors and is a diagnostic baseline, not a bound on the optimized ELBO.
noisybool: Whether the stack is treated as noisy.
noise_status_sourcetuple of str: How each retained run’s noise status was obtained: recorded, inferred or override, in retained-run order.
entropyfloat or None: After optimize(): the entropy estimate of the stacked posterior at the optimized weights.
Inp.ndarray, shape (1, K_total): Expected log-joint of every component (the mean over the run’s GP hyperparameter samples of stats["I_sk"]), in the run’s own transformed space.
individual_elboslist of float: The ELBO of each retained run, used to initialize the weights.
I_corrected, E_correctednp.ndarray: The expected log-joints corrected to the original space, per component (shape (1, K_total)) and per run (shape (M,)), computed deterministically at construction. The capped ELBO variants limit the expected log-joint to their medians.
rngnp.random.Generator: The object’s random generator.
loggerlogging.Logger: Progress goes to the "SVBMC" logger at INFO level; raise its level (to WARNING, say) to quieten a run.

classmethod load(file)[source]#

Load a stacked posterior from a file written by save().

Loading needs no torch: sample() and plot() of the loaded object work without it, while optimize() and the methods that estimate the stacked ELBO or its entropy raise ImportError as construction does. The logging level is not saved: progress goes to the "SVBMC" logger at INFO level unless a level has been set on it, as for a newly constructed object.

Parameters:

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

Returns:

stackedSVBMC: The loaded object.

Raises:

TypeError: If the file holds another kind of object, such as a VariationalPosterior.
OSError: If the file cannot be found, or cannot be opened for other reasons.

maximize_ELBO(n_samples: int = 20, lr: float = 0.1, max_steps: int = 500, version: str = 'all-weights')[source]#

Maximize the stacked ELBO over the weights with Adam.

The weights are the softmax of unconstrained logits, initialized from the runs’ own weights and ELBOs so that better runs start heavier. That starting point belongs to the stack, so every call departs from it whatever w currently holds. Every step draws fresh entropy samples. The optimization stops when the ELBO, rounded to five decimals, has not improved for five consecutive steps, or after max_steps; the best iterate is returned.

Parameters:

n_samplesint, optional: Draws per component for the entropy estimate. Default 20.
lrfloat, optional: Adam learning rate. Default 0.1.
max_stepsint, optional: Maximum number of steps. Default 500.
version{“all-weights”, “posterior-only”, “ns”}, optional: "all-weights" optimizes the weight of every component; "posterior-only" optimizes one weight per run, keeping the relative weights within a run; "ns" (naive stacking) performs no optimization and returns the normalized input weights. Default "all-weights".

Returns:

w_finaltorch.Tensor, shape (K_total,): The optimized weights (float64).
elbo_besttorch.Tensor: Optimization-time Monte Carlo ELBO from the selected iteration (float64 scalar). optimize() re-evaluates the returned weights with fresh draws for the reported headline.
entropy_besttorch.Tensor: The entropy estimate at the returned weights (float64 scalar).

optimize(n_samples: int = 20, lr: float = 0.1, max_steps: int = 500, version: str = 'all-weights', *, n_samples_final: int = 100)[source]#

Optimize the stacked posterior in place.

Runs maximize_ELBO(), stores the optimized weights in w (NumPy float64), the entropy estimate in entropy and a fresh final evaluation in elbo, elbo_sd and elbo_details. The report includes a deterministic two-level shrinkage estimate at those weights, computed with the same final entropy and without random draws; where it is defined, it is the headline of a noisy stack. Returns None.

Parameters:

n_samplesint, optional: Draws per component for the entropy estimate. Default 20.
lrfloat, optional: Adam learning rate. Default 0.1.
max_stepsint, optional: Maximum number of steps. Default 500.
version{“all-weights”, “posterior-only”, “ns”}, optional: Optimization mode; see maximize_ELBO(). Default "all-weights".
n_samples_finalint, optional: Fresh entropy draws per component at the returned weights and at the naive-stacking weights. At least 2, to estimate sampling variance. Default 100. Naive mode reuses one final evaluation.

plot(n_samples: int = 10000, color: str = 'black', figsize: tuple | None = None, smooth: float | None = None, **corner_kwargs)[source]#

Draw a corner plot of the stacked posterior from fresh samples.

Parameters:

n_samplesint, optional: Draws taken with sample(). Default 10000.
colorstr, optional: Matplotlib color of the contours and histograms. Default "black".
figsize(float, float) or None, optional: Figure size in inches; None uses 2.5 inches per dimension.
smoothfloat or None, optional: Gaussian kernel smoothing applied by corner.
**corner_kwargs: Passed to corner.corner().

Returns:

matplotlib.figure.Figure: The figure.

sample(n_samples: int, balance_flag: bool = False)[source]#

Draw samples from the stacked posterior in the original space.

The stacked posterior is a mixture over the retained runs, so a draw picks a run with probability equal to that run’s total weight and then a point from that run’s reweighted posterior through its own transform. By default the run counts are multinomial, the component choice within a run is random (see VariationalPosterior.sample()) and the rows are shuffled, so the result is an independent sample of exactly n_samples rows. With balance_flag=True the draws are stratified instead: the number taken from a run is within one draw of that run’s exact share of the weights, and within a run they are spread over its components by the balanced allocation of VariationalPosterior.sample(), which follows the component weights up to the random placement of its remainder and so holds no exact bound per component. Stratification lowers the variance of expectations; it does not give an independent sample. All randomness comes from this object’s generator; the input posteriors are not touched.

Parameters:

n_samplesint: Number of draws.
balance_flagbool, optional: Stratify the draws by run and component. Default False.

Returns:

Xnp.ndarray, shape (n_samples, D): Draws in the original parameter space.

save(file, overwrite=False)[source]#

Save the stacked posterior to a file.

The file holds the whole object: the retained posteriors, the weights, the ELBO report and the state of the generator, so the object that load() returns draws what this one would draw next. It holds no Python bytecode, so it can be loaded, used and saved again under another minor version of Python than the one that wrote it.

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.).

stacked_ELBO(w, n_samples: int = 20)[source]#

Estimate the stacked ELBO at given weights.

The expected log-joint is the weighted sum of the components’ stored expected log-joints (from the runs’ GP surrogates), corrected to the original space with the Jacobian terms of stacked_entropy(); the entropy is estimated by Monte Carlo. The expected log-joints are corrected deterministically during construction.

Parameters:

wtorch.Tensor, np.ndarray or list: Weights of the stacked posterior (K_total values; normalized internally). NumPy or list input is converted to a float64 tensor.
n_samplesint, optional: Draws per component for the entropy estimate. Default 20.

Returns:

ELBOtorch.Tensor: The estimated stacked ELBO (float64 scalar).
Htorch.Tensor: The estimated entropy (float64 scalar).

stacked_entropy(w, n_samples: int = 20)[source]#

Monte Carlo estimate of the entropy of the stacked posterior.

Differentiable with respect to the mixture weights. Every component contributes n_samples draws, taken in its own run’s transformed space and mapped to the original space; the log density of the stacked mixture at every draw is then evaluated by mapping the draws into each run’s space with the matching Jacobian correction. The draws come from this object’s generator.

Parameters:

wtorch.Tensor: Weights of the stacked posterior (any shape holding K_total values; normalized internally). Converted to float64.
n_samplesint, optional: Draws per component. Default 20.

Returns:

Htorch.Tensor: The entropy estimate (a float64 scalar).
J_correctionsnp.ndarray, shape (K_total,): Per-component expected log-Jacobian, computed deterministically at construction. Subtracting it from the stored expected log-joints expresses them in the original space, which I_corrected already holds and stacked_ELBO() reads; this copy is for inspection.

Helpers for comparing VBMC runs and choosing their starting points.

pyvbmc.svbmc.utils.find_init_bounds(LB: ndarray = None, UB: ndarray = None, PLB: ndarray = None, PUB: ndarray = None)[source]#

Find the bounds to sample uniformly from when choosing starting points for VBMC.

If plausible lower bounds (PLB) are specified, it uses those as lower sampling bounds (sample_LB), otherwise it uses lower bounds (LB). Upper sampling bounds (sample_UB) are determined in the same way with plausible upper bounds (PUB) or upper bounds (UB).

At least one of LB, UB, PLB and PUB must be specified as an array with the same dimensionality of the inference problem. Each bound may be a scalar (broadcast to every dimension) or an array holding D values in any layout, including the (1, D) rows PyVBMC uses; the returned bounds are flat (D,) arrays.

Parameters:

LB, UBnp.ndarray, optional: Inputs for VBMC. They represent lower (LB) and upper (UB) 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 [d] = -inf and UB [d] = inf if the d-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). If PLB and PUB are not specified (see below), the lower and upper sampling bounds (sample_LB and sample_UB, respectively) will be determined by LB and UB. Both are by default None.
PLB, PUBnp.ndarray, optional: Inputs for VBMC. They represent a set of plausible lower (PLB) and upper (PUB) bounds 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. Among other things, the plausible box is used by VBMC to draw initial samples and to set priors over hyperparameters of the algorithm. If PLB and PUB are specified, they determine the lower and upper sampling bounds (sample_LB and sample_UB, respectively). If they are not, the sampling bounds will be determined by LB and UB. Both are by default None.

Returns:

sample_LB, sample_UBnp.ndarray: Lower (sample_LB) and upper (sample_UB) bounds to sample uniformly from when initializing VBMC.

pyvbmc.svbmc.utils.overlay_corner_plot(samples: List[ndarray], labels: List[str] | None = None, colors: List[str] | None = None, figsize: Tuple[float, float] | None = None, smooth: float | None = None, base: float = 2.5, **corner_kwargs)[source]#

Overlay a corner plot for multiple sample sets.

Parameters:

sampleslist[np.ndarray]: List of arrays, each of shape (N_i,`D`) where D is the dimensionality.
labelslist[str], optional: Legend labels corresponding to each item in samples. Defaults to “Run 1”, “Run 2”, …
colorslist[str], optional: Matplotlib colors for each sample set. Defaults to the rcParams color cycle.
figsize(float, float) or None, optional: Figure size in inches. If None, uses base`*`D on each side.
smoothfloat or None, optional: Gaussian kernel smoothing (px) applied by corner.
basefloat, default 2.5: Inches per variable when auto-sizing.
**corner_kwargs: Extra keywords forwarded to corner.corner() (e.g. bins=40, levels=(0.68, 0.95)).

Returns:

matplotlib.figure.Figure: The resulting corner-plot figure.

SVBMC

Contents

SVBMC#

Basic usage#

ELBO reporting#

A composite posterior#

Helpers#

Differences from the standalone svbmc package#

Citation#

`SVBMC`#

Differences from the standalone `svbmc` package#