Getting started

Getting started#

The best way to get started with PyVBMC is via the tutorials and worked examples. In particular, start with PyVBMC Example 1: Basic usage and continue from there.

Optionally, we recommend running pyvbmc.calibrate() once to tune performance for your machine (memory, processor, etc.); it takes tens of seconds and can make PyVBMC run faster. Saved calibration settings are reused automatically by future runs (see the calibration guide).

If you are already familiar with approximate inference methods, you can find a summary usage below.

Summary usage#

The typical usage pipeline of PyVBMC follows four steps:

  1. Define the model, which defines a target log density (i.e., an unnormalized log posterior density);

  2. Set up the parameters (parameter bounds, starting point);

  3. Initialize and run the inference;

  4. Examine and visualize the results.

PyVBMC is not concerned with how you define your model in step 1, as long as you can provide an (unnormalized) target log density. Running the inference in step 3 only involves a couple of lines of code:

from pyvbmc import VBMC
# ...
vbmc = VBMC(target, x0, LB, UB, PLB, PUB)
vp, results = vbmc.optimize()

with input arguments:

  • target: the target (unnormalized) log density — often an unnormalized log posterior. target takes as input a parameter vector and returns the log density at the point. The returned log density must be a finite real value, i.e. neither NaN nor +/-inf. See the FAQ for more details;

  • x0: the starting point of the inference in parameter space;

  • LB and UB: hard lower and upper bounds for the parameters (can be -inf and inf, or bounded);

  • PLB and PUB: plausible lower and upper bounds, that is a box that ideally brackets a region of high density of the target.

You can also pass a log-likelihood as target and supply its prior separately with prior= or log_prior=. For model evidence, include the likelihood’s normalization constants and use a normalized prior; see the FAQ.

The outputs are:

  • vp: a VariationalPosterior object which approximates the true target density;

  • results: a dict with additional information. Important keys are: - "elbo": the estimated lower bound on the log model evidence (log normalization constant); - "elbo_sd": the standard deviation of the estimate of the ELBO (not the error between the ELBO and the true log model evidence, which is generally unknown).

The vp object can be manipulated in various ways, see the VariationalPosterior class documentation.

See the examples for more detailed information. The Basic options may also be useful.

For noisy likelihoods, see PyVBMC Example 6: Noisy log-likelihood evaluations and the FAQ for the target interface and noise options.

PyVBMC occasionally prints a tip or a reminder when a new run starts. Pass options={"show_tips": False} to disable tips and reminders about old releases. To suppress the calibration reminder and ordinary optimization output as well, use options={"display": "off"}.

Reproducible runs#

Pass seed= to control the random generator used by a VBMC instance and its variational posterior. The same integer seed reproduces a run with the same target, setup and numerical environment. Reusing the same integer seed repeats PyVBMC’s random stream. For independent runs, leave seed unset or use different seeds; comparing their results helps assess convergence (see PyVBMC Example 4: Multiple runs as validation).

The seed can also be a NumPy Generator. PyVBMC advances that generator, so reusing it continues its random stream. The posterior shares the run’s generator: drawing samples from it before resuming inference changes the subsequent draws of the run.

If your likelihood function uses random simulations, VBMC(seed=...) does not control the randomness inside that function. To reproduce the whole analysis, also set the simulator’s seed before starting the run. For example, with NumPy, create sim_rng = np.random.default_rng(123) once outside your target function, then use that generator for its random draws (e.g., sim_rng.normal(...) instead of np.random.normal(...)). If you use another simulation library, use its seed or generator argument. Let the simulator draw fresh random numbers on each evaluation; do not reset its seed inside the target function.

The FAQ explains reproducibility and the behavior when seed is omitted.

Bring a torch or JAX model into PyVBMC#

The default target interface calls one point at a time. It accepts a one-dimensional NumPy parameter vector and returns one finite scalar. A model implemented in Torch can be connected with a small user-written wrapper that keeps framework objects outside PyVBMC. The example below uses independent Gaussian observations with unit standard deviation and a Gaussian prior on their common mean. Install the modelling framework used by your model separately:

import numpy as np
import torch
from pyvbmc import VBMC
from scipy.stats import norm

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
observations_t = torch.tensor(
    [0.4, -0.2, 0.7], dtype=torch.float64, device=device
)

def torch_log_likelihood(x):
    theta_t = torch.as_tensor(x, dtype=torch.float64, device=device)
    value_t = -0.5 * torch.sum(
        (observations_t - theta_t[0]) ** 2 + np.log(2 * np.pi)
    )
    return np.asarray(
        value_t.detach().cpu().numpy(), dtype=np.float64
    ).item()

x0 = np.array([0.0], dtype=np.float64)
lb, ub = np.array([-np.inf]), np.array([np.inf])
plb, pub = np.array([-2.0]), np.array([2.0])
vbmc = VBMC(
    torch_log_likelihood,
    x0,
    lb,
    ub,
    plb,
    pub,
    prior=norm(loc=0.0, scale=2.0),
)

detach().cpu().numpy() completes device execution and returns the scalar to the host. This also makes PyVBMC’s target timing cover the framework work. The explicit torch.float64 matches PyVBMC’s numerical precision.

For JAX, enable 64-bit values before importing jax.numpy or constructing the model. jax.device_get completes the computation and moves its result to the host:

from jax import config
config.update("jax_enable_x64", True)

import jax
import jax.numpy as jnp
import numpy as np

observations_j = jnp.array([0.4, -0.2, 0.7], dtype=jnp.float64)

def jax_log_likelihood(x):
    theta_j = jnp.asarray(x, dtype=jnp.float64)
    value_j = -0.5 * jnp.sum(
        (observations_j - theta_j[0]) ** 2 + jnp.log(2 * jnp.pi)
    )
    return np.asarray(
        jax.device_get(value_j), dtype=np.float64
    ).item()

Vectorized targets#

Set options={"vectorized_target": True} when the likelihood can evaluate a batch. PyVBMC collects all initial-design rows whose values were not supplied and evaluates them in one target call. All subsequent target calls still contain one point, with shape (1, D). The input is always a NumPy array in original coordinates and the target must return a finite NumPy array with shape (N,) or (N, 1). PyVBMC does not probe a target to determine whether it supports batches. A scalar return is invalid even when N is one.

The Torch wrapper above becomes:

def torch_vectorized_log_likelihood(x):
    theta_t = torch.as_tensor(x, dtype=torch.float64, device=device)
    values_t = -0.5 * torch.sum(
        (observations_t[None, :] - theta_t[:, :1]) ** 2
        + np.log(2 * np.pi),
        dim=1,
    )
    return np.asarray(
        values_t.detach().cpu().numpy(), dtype=np.float64
    )

vbmc = VBMC(
    torch_vectorized_log_likelihood,
    x0,
    lb,
    ub,
    plb,
    pub,
    prior=norm(loc=0.0, scale=2.0),
    options={"vectorized_target": True},
)

The equivalent JAX wrapper is:

def jax_vectorized_log_likelihood(x):
    theta_j = jnp.asarray(x, dtype=jnp.float64)
    values_j = -0.5 * jnp.sum(
        (observations_j[None, :] - theta_j[:, :1]) ** 2
        + jnp.log(2 * jnp.pi),
        axis=1,
    )
    return np.asarray(jax.device_get(values_j), dtype=np.float64)

A separately supplied prior keeps its scalar interface. PyVBMC evaluates it one point at a time and combines it with the batched likelihood; custom prior functions therefore do not need to accept a batch.

For a target with user-provided noise estimates, set "specify_target_noise": True and return a pair of arrays. Each array must have shape (N,) or (N, 1); the standard deviations must be finite and positive. Do not combine the two outputs into an (N, 2) array. For example:

def torch_noisy_log_likelihood(x):
    values = torch_vectorized_log_likelihood(x)
    sds = np.full(values.shape, 0.1, dtype=np.float64)
    return values, sds

def jax_noisy_log_likelihood(x):
    values = jax_vectorized_log_likelihood(x)
    sds_j = jnp.full(values.shape, 0.1, dtype=jnp.float64)
    sds = np.asarray(jax.device_get(sds_j), dtype=np.float64)
    return values, sds

noisy_options = {
    "vectorized_target": True,
    "specify_target_noise": True,
}

Pass noisy_options as the options argument when constructing VBMC with either noisy wrapper.

To reuse cached target values, pass precomputed_evaluations=(X, y) to VBMC. Each value in y must be the output of the target supplied to VBMC at the corresponding row of X. When a prior is supplied separately, the target and y contain the log-likelihood and PyVBMC adds the prior; when the target is a log joint, y contains that log joint. See Reusing evaluations and charging initialization work for shapes, noisy observations and budget accounting.

For a complete Torch fit, an equivalent short JAX target, Torch posterior predictions and density gradients, and ArviZ summaries, see PyVBMC Example 9: Torch and JAX models and posterior exports.

Bring a PyMC model into PyVBMC#

Install the optional integration on Python 3.12 or newer with python -m pip install --upgrade "pyvbmc[pymc]". A PyMCTarget provides the flat target and setup that PyVBMC needs. This model has a named coefficient vector and a positive observation scale:

import numpy as np
import pymc as pm
from pyvbmc import PyMCTarget, VBMC
rng = np.random.default_rng(7)
x_data = np.linspace(-1.0, 1.0, 30)
design = np.column_stack((np.ones(x_data.size), x_data))
observations = design @ np.array([0.5, 1.0])
observations += rng.normal(0.0, 0.7, size=x_data.size)
coords = {
    "observation": np.arange(x_data.size),
    "coefficient": ["intercept", "slope"],
}
with pm.Model(coords=coords) as model:
    X = pm.Data("X", design, dims=("observation", "coefficient"))
    beta = pm.Normal("beta", 0.0, 2.0, dims="coefficient")
    sigma = pm.HalfNormal("sigma", 1.0)
    mean = pm.Deterministic("mean", X @ beta, dims="observation")
    pm.Normal("y", mean, sigma, observed=observations, dims="observation")
target = PyMCTarget(model, seed=7)
print(target)
vbmc = VBMC(
    target,
    seed=7,
    options={"max_fun_evals": target.setup_cost + 60},
)
vp, results = vbmc.optimize()
posterior_data = target.to_arviz(vp, n_samples=1000)
posterior_data = pm.compute_deterministics(
    posterior_data, model=model, extend_dataset=True
)
predictions = pm.sample_posterior_predictive(
    posterior_data, model=model, var_names=["y"], random_seed=7
)

print(target) lists each model variable, its shape, its flat VBMC coordinates, hard and plausible bounds, and its start. Here beta becomes beta[intercept] and beta[slope]. The positive sigma keeps PyMC’s log transform, so its VBMC coordinate is sigma_log__ and the target includes the corresponding Jacobian.

The adapter chooses a starting point and plausible bounds automatically. Its setup work counts toward max_fun_evals through target.setup_cost; the example allows 60 more evaluations for inference. Use a larger budget and validate repeated runs for your own model. Pass setup_budget= to limit the setup work, or start= and plausible_bounds= to supply your own choices in model coordinates.

The adapter supports continuous float64 variables with recognized support and standard PyMC transforms. It does not support discrete variables. See PyMCTarget for the supported models and setup options.

PyMCTarget keeps the model data used at construction. Build a new target to refit changed data, and export each posterior with the target used to fit it. PyMCTarget(seed=) controls setup draws; VBMC(seed=) controls the inference run. PyVBMC Example 8: Fitting a PyMC model shows the full workflow.

Use a fitted posterior downstream#

Stacking several runs#

Stacking Variational Bayesian Monte Carlo (S-VBMC) (Silvestrin et al., 2025) combines the posteriors of several completed VBMC runs on the same model and data into a stacked posterior, without further model evaluations. Combining runs often improves the approximation to the true posterior. See PyVBMC Example 7: Stacking the posteriors of several runs (S-VBMC) for a worked example. S-VBMC requires the optional torch extra; the installation guide also covers dependencies for the posterior exports below.

Torch distribution#

to_torch() exports the fitted posterior as a torch distribution:

import torch

posterior_t = vp.to_torch()
torch.manual_seed(7)
samples_t = posterior_t.sample((1000,))

point_t = samples_t[0].detach().clone().requires_grad_(True)
gradient_t = torch.autograd.grad(posterior_t.log_prob(point_t), point_t)[0]

The default export uses original parameter coordinates, CPU tensors, and torch.float64, regardless of torch’s global defaults. Pass orig_flag=False for PyVBMC’s unbounded internal coordinates, or specify dtype=torch.float32 and device=... explicitly when needed. The export is an independent copy; sampling it uses torch’s random generator.

For bounded parameters, log_prob accepts finite values strictly inside the hard bounds. Support validation rejects exact-bound and outside values. See the method reference for supported shapes and boundary behavior.

The complete workflow in PyVBMC Example 9: Torch and JAX models and posterior exports uses this distribution for psychometric probability predictions and density gradients.

ArviZ DataTree#

to_arviz() draws independent samples into the current ArviZ DataTree format:

import arviz as az

posterior_data = vp.to_arviz(
    n_samples=2000, var_names=["location"]
)
summary = az.summary(posterior_data, group="posterior", kind="stats")
axes = az.plot_dist(posterior_data, group="posterior")

The result contains one posterior group and one chain. By default, each parameter is a scalar variable named x_0, …, x_{D-1}; the structured form below supports vector and matrix variables. Pass orig_flag=False to export samples in internal coordinates.

For a three-dimensional posterior whose first two columns form a vector, preserve the variable shape and labels with a structured export:

# Here vp has D == 3.
structured_data = vp.to_arviz(
    n_samples=2000,
    variables={"beta": (2,), "sigma": ()},
    dims={"beta": ["coef"]},
    coords={"coef": ["intercept", "slope"]},
)

Here beta has dimensions (chain, draw, coef) and sigma has dimensions (chain, draw). Variables consume consecutive posterior columns in mapping order.

This method advances vp.rng in the same way as sample(). The draws are independent samples from a variational approximation, so MCMC convergence diagnostics computed on them do not measure the quality of that approximation. Use posterior summaries and plots, and assess PyVBMC convergence from its own results and repeated runs.