VariationalPosterior

Note

The main output of a PyVBMC run is a VariationalPosterior object vp.

The variational posterior can be queried and manipulated via several methods, such as:

• moments: compute mean and covariance matrix of the variational posterior;

• pdf: evaluate the variational posterior density at a point;

• plot: plot the variational posterior as a corner plot (1D and 2D marginals);

• sample: draw random samples from the variational posterior.

• to_torch(): export an independent torch distribution;

• to_arviz(): export independent draws in the current ArviZ DataTree format.

There are also methods to compare variational posteriors:

• kl_div: compute the Kullback-Leibler divergence between two posteriors;

• mtv: compute the marginal total variation distance between two posteriors.

See below for the full VariationalPosterior class methods and interface.

Density gradients from pdf and log_pdf are available in transformed coordinates (orig_flag=False). Requests for gradients in the original parameter space raise NotImplementedError.

Performance settings

A fitted posterior retains the fixed CalibrationProfile used by its VBMC run. Standalone posteriors use historical defaults; pass calibration=profile or calibration="cached" to opt in. A cache lookup never initiates a campaign. See calibrate for measuring settings and the reproducibility contract.

class pyvbmc.variational_posterior.VariationalPosterior(D: int, K: int = 2, x0=None, parameter_transformer=None, rng=None, calibration=None)

The variational posterior class used in PyVBMC.

The variational posterior represents the approximate posterior as returned by the Variational Bayesian Monte Carlo (VBMC) algorithm.

In VBMC, the variational posterior is a weighted mixture of multivariate normal distributions (see Notes below for details).

Parameters:

Dint

The number of dimensions (i.e., parameters) of the posterior.

Kint, optional

The number of mixture components, default 2.

x0np.ndarray, optional

The starting points for the mixture component means, one per row. A single point of D elements, given as a flat array, a row or a column, starts every component; a 1-by-1 array starts every coordinate of every component; an n0-by-D matrix starts the components at its rows, repeated in order where n0 is below K and cut after the K-th row where it is above. By default np.zeros.

parameter_transformerParameterTransformer, optional

The ParameterTransformer object specifying the transformation of the input space that leads to the current representation used by the variational posterior, by default uses an identity transform.

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

The random generator, or a seed for one, used by every method that draws random numbers (e.g. sample, kl_div, mtv). A Generator is used as is, so it can be shared with a VBMC instance. By default a generator is derived from NumPy’s global random state, so that np.random.seed still gives reproducible results.

calibrationCalibrationProfile, {“cached”, “off”}, or None, optional

Fixed numerical chunk settings. None and "off" use the historical defaults. "cached" defers a machine-local cache lookup until the first density or Monte Carlo entropy evaluation.

Attributes:

wnp.ndarray

The weights of the VP mixture components, shape (1, K).

etanp.ndarray

The unbounded (softmax) parametrization of the VP mixture components, shape (1, K).

munp.ndarray

The means of the VP mixture components, shape (D, K).

sigmanp.ndarray

The per-component scale of the VP mixture components. Shape (1, K).

lambdnp.ndarray

The per-dimension scale of the VP mixture components. Shape (D, 1).

optimize_weightsbool

Whether to optimize the weights.

optimize_mubool

Whether to optimize the means.

optimize_sigmabool

Whether to optimize sigma.

optimize_lambdbool

Whether to optimize lambd.

parameter_transformerParameterTransformer

The parameter transformer implementing transformations to/from unbounded space.

boundsdict

A dictionary containing the soft bounds for each variable to be optimized, as get_bounds last computed them, or None before the first call.

statsdict

A dictionary of statistics and other relevant info computed during optimization.

calibration_profileCalibrationProfile or None

The fixed profile, or None before a cached request resolves.

Notes

In VBMC, the variational posterior is defined as a mixture of multivariate normal distributions as follows:

\[q({\theta}) = \sum_{k = 1}^K w_k N(\theta, \mu_k, \sigma^2_k \Lambda)\]

where \(w_k\) are the mixture weights, \(\mu_k\) the component means, \(\sigma_k\) scaling factors, and \(\Lambda\) is a diagonal matrix common to all components with elements \(\lambda^2_d\) on the diagonal, for \(1 \le d \le D\).

Note that \(q({\theta})\) is defined in an unconstrained space. Constrained variables in the posterior are mapped to a transformed, unconstrained space via a nonlinear mapping (represented by a ParameterTransformer object). The transformation is handled automatically.

In practice, you would almost never create a new VariationalPosterior object, simply use the variational posteriors returned by VBMC.

property calibration_profile: CalibrationProfile | None

The fixed profile, or None before a cached request resolves.

get_bounds(X: ndarray, options, K: int = None)

Compute soft bounds for variational posterior parameters.

These bounds are used during the variational optimization in VBMC. They are a function of the training inputs given here: the component means are bounded by the bounding box of those inputs, and their log scale (the product of sigma and lambd) by the log of the width of that box, down to a fraction tol_length of it. The bounds on the weight parameters span the weights above half the weight tolerance. The box is also stored on the posterior as bounds, replacing any box of an earlier call.

Parameters:

Xndarray, shape (N, D)

Training inputs.

optionsOptions

Program options.

Kint, optional

The number of mixture components. By default we use the number provided at class instantiation.

Returns:

theta_bnddict

A dictionary of soft bounds with the following elements:

lbnp.ndarray

Lower bounds.

ubnp.ndarray

Upper bounds.

tol_confloat

Fractional tolerance for constraint violation of variational parameters.

weight_thresholdfloat, optional

Threshold below which weights are penalized.

weight_penaltyfloat, optional

The penalty for weight below the threshold.

get_parameters(raw_flag=True)

Get variational posterior parameters as single array.

Return all the active VariationalPosterior parameters flattened as a 1D (numpy) array, possibly transformed.

The parameters are first normalized in place: lambd is divided by its root mean square and sigma multiplied by it, and the weights are divided by their sum when they are being optimized. The two scales enter the density through their product alone, so rescaling does not change the distribution. Normalizing the weights changes it only if they did not already sum to one. A mode cached by a previous call to mode() is discarded because the parameters may have changed.

Parameters:

raw_flagbool, optional

Specifies whether the scales sigma and lambd and the weights w are returned as raw (unconstrained) parameters, their logarithms, or as they are, by default True.

Returns:

thetanp.ndarray

The variational posterior parameters flattened as a 1D array.

kl_div(vp2: VariationalPosterior = None, samples: ndarray = None, N: int = 100000, gauss_flag: bool = False)

Kullback-Leibler divergence between two variational posteriors.

Compute the forward and reverse Kullback-Leibler (KL) divergence between two posteriors. The other variational posterior can be specified as vp2 (an instance of the class VariationalPosterior) or with samples. One of the two must be specified.

Parameters:

vp2VariationalPosterior, optional

The other VariationalPosterior, by default None.

samplesnp.ndarray, optional

An N-by-D matrix of samples from the other variational posterior, by default None.

Nint, optional

The number of random samples to estimate the KL divergence, by default int(1e5).

gauss_flagbool, optional

If True, returns a “Gaussianized” KL-divergence, that is the KL divergence between two multivariate normal distributions with the same moments as the variational posteriors given as inputs. By default False.

Returns:

kl_div: np.ndarray

A two-element vector containing the forward and reverse Kullback-Leibler divergence between the two posteriors.

Raises:

ValueError

Raised if neither vp2 nor samples are specified.

ValueError

Raised if vp2 is not provided but gauss_flag = False.

ValueError

Raised if gauss_flag is True and N is 0: the moments are estimated from N samples in the original space.

Notes

Since the KL divergence is not symmetric, the method returns both the forward and the reverse KL divergence, that is KL(vp1 || vp2) and KL(vp2 || vp1).

classmethod load(file, *, calibration=None)

Load a VP from a file.

Parameters:

filepath-like

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

calibrationCalibrationProfile, {“cached”, “off”}, or None, optional

For a posterior saved while its calibration request ("cached") was still pending, the settings that replace the request; None leaves it pending. For a posterior with resolved settings, None, a profile with the same settings, and "off" where the saved settings are the historical defaults keep the saved profile and its provenance; any other value raises.

Returns:

vpVariationalPosterior

The loaded VP instance.

Raises:

OSError

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

ValueError

If calibration asks for other settings than the resolved ones the file holds ("cached", "off" where the saved settings are not the historical defaults, or a profile with other settings), or is not one of the accepted values.

log_pdf(*args, **kwargs)

log-probability density function of the variational posterior.

Compute the log density of the variational posterior at one or multiple input points. The parameters are the same as for vp.pdf() with log_flag=True. These parameters are described again here for reference.

Parameters:

xnp.ndarray

x is a matrix of inputs to evaluate the pdf at. The rows of the N-by-D matrix x correspond to observations or points, and columns correspond to variables or coordinates. x is assumed to be in the original space by default. A one-dimensional array is one point, so for D = 1 several points go in a column.

orig_flagbool, optional

Controls if the value of the posterior density should be evaluated in the original parameter space for orig_flag is True, or in the transformed space if orig_flag is False, by default True. Accordingly, x should be in the original space if orig_flag is True and be in the transformed space if orig_flag is False.

grad_flagbool, optional

If True the gradient of the log-pdf is returned as a second output, by default False. Gradients are available only in transformed space (orig_flag=False).

dffloat, optional

Compute the log-pdf of a heavy-tailed version of the variational posterior, in which the multivariate normal components have been replaced by multivariate t-distributions with df degrees of freedom. The default is df = np.inf, limit in which the t-distribution becomes a multivariate normal. A negative df replaces them instead by products of D univariate t-distributions with abs(df) degrees of freedom; df = 0 also gives the normal components.

Returns:

log_pdf: np.ndarray

The log probability density of the variational posterior evaluated at each row of x.

gradient: np.ndarray

If grad_flag is True, the function returns the gradient as well.

Raises:

ValueError

Raised if the rows of x are not D wide.

NotImplementedError

Raised if df is non-zero and finite and grad_flag = True (Gradient of heavy-tailed pdf not supported yet).

NotImplementedError

Raised if orig_flag = True and grad_flag = True (gradient computation in original space is not supported).

mode(orig_flag=True, n_opts: int | None = None)

Find the mode of the variational posterior.

Parameters:

orig_flagbool, optional

If True find the mode of the variational posterior in the original parameter space, otherwise in the transformed parameter space. By default True.

n_optsint, optional

Maximum number of optimization runs from different starting points to find the mode. By default n_opts is the square root of the number of mixture components K, that is \(n\_opts = \lceil \sqrt{K} \rceil\). A call in the original space that leaves n_opts out returns the mode that an earlier such call found, if the posterior still carries one, and otherwise stores the mode it finds. A call in the transformed space, or one that gives n_opts, runs the search and neither reads nor replaces that stored mode.

Returns:

modenp.ndarray

The mode of the variational posterior.

Notes

Mode estimation (e.g., for the purpose of maximum-a-posteriori estimation) is not recommended with VBMC, since due to the underlying representation (mixture of Gaussians) the mode of the variational posterior is a brittle and potentially unreliable estimator of the mode of the target posterior, especially if it lies close to the boundaries of the space.

The mode is not invariant to nonlinear reparameterizations of the input space, so the mode in the original space and the mode in the transformed (unconstrained) space will generally be in different locations (even after applying the appropriate transformations).

The starting points of the optimization runs are drawn from a copy of self.rng, so a call leaves the posterior’s random stream where it found it, and two calls with nothing drawn from that stream in between give the same answer.

moments(N: int = 1000000, orig_flag=True, cov_flag=False)

Compute mean and covariance matrix of variational posterior.

In the original space, the moments are estimated via Monte Carlo sampling. If requested in the transformed (unconstrained) space used internally by VBMC, the moments can be computed analytically.

Parameters:

Nint, optional

Number of samples used to estimate the moments in the original space, as any scalar that holds a whole number (see sample). By default int(1e6).

orig_flagbool, optional

If True, compute moments in the original parameter space, otherwise in the transformed VBMC space. By default True.

cov_flagbool, optional

If True, return the covariance matrix as a second return value, by default False.

Returns:

mean: np.ndarray

The mean of the variational posterior, of shape (1, D).

cov: np.ndarray

If cov_flag is True, returns the covariance matrix as well, of shape (D, D).

Raises:

ValueError

Raised in the original space if N is not a scalar holding a whole number of samples.

mtv(vp2: VariationalPosterior = None, samples: ndarray = None, N: int = 100000)

Marginal total variation distances between two variational posteriors.

Compute the total variation distance between the variational posterior and a second posterior, separately for each dimension (hence “marginal” total variation distance, MTV). The second posterior can be specified either as a VariationalPosterior or as a set of samples.

Parameters:

vp2VariationalPosterior, optional

The other VariationalPosterior, by default None.

samplesnp.ndarray, optional

An N-by-D matrix of samples from the other variational posterior, by default None.

Nint, optional

The number of random draws to estimate the MTV, by default int(1e5).

Returns:

mtvnp.ndarray

An array of shape (1, D) holding, for each coordinate, the total variation distance between the marginal distributions of this posterior and of vp2 or samples.

Raises:

ValueError

Raised if neither vp2 nor samples are specified.

Notes

The total variation distance between two densities \(p_1\) and \(p_2\) is:

\[TV(p_1, p_2) = \frac{1}{2} \int | p_1(x) - p_2(x) | dx.\]

pdf(x: ndarray, orig_flag: bool = True, log_flag: bool = False, grad_flag: bool = False, df: float = inf)

Probability density function of the variational posterior.

Compute the probability density function (pdf) of the variational posterior at one or multiple input points.

Parameters:

xnp.ndarray

x is a matrix of inputs to evaluate the pdf at. The rows of the N-by-D matrix x correspond to observations or points, and columns correspond to variables or coordinates. x is assumed to be in the original space by default. A one-dimensional array is one point, so for D = 1 several points go in a column.

orig_flagbool, optional

Controls if the value of the posterior density should be evaluated in the original parameter space for orig_flag is True, or in the transformed space if orig_flag is False, by default True. Accordingly, x should be in the original space if orig_flag is True and be in the transformed space if orig_flag is False.

log_flagbool, optional

If log_flag is True return the logarithm of the pdf, by default False. Where the density of the mixture of normal components is below the smallest normal double, its logarithm and the gradient of the logarithm are computed by log-sum-exp over the components, so that they stay finite in the far tails.

grad_flagbool, optional

If True the gradient of the pdf is returned as a second output, by default False. Gradients are available only in transformed space (orig_flag=False).

dffloat, optional

Compute the pdf of a heavy-tailed version of the variational posterior, in which the multivariate normal components have been replaced by multivariate t-distributions with df degrees of freedom. The default is df = np.inf, limit in which the t-distribution becomes a multivariate normal. A negative df replaces them instead by products of D univariate t-distributions with abs(df) degrees of freedom; df = 0 also gives the normal components.

Returns:

pdf: np.ndarray

The probability density of the variational posterior evaluated at each row of x.

gradient: np.ndarray

If grad_flag is True, the function returns the gradient as well.

Raises:

ValueError

Raised if the rows of x are not D wide.

NotImplementedError

Raised if df is non-zero and finite and grad_flag = True (Gradient of heavy-tailed pdf not supported yet).

NotImplementedError

Raised if orig_flag = True and grad_flag = True (gradient computation in the original space is not supported).

plot(n_samples: int = 100000, title: str = None, plot_data: bool = False, highlight_data: np.ndarray = None, plot_vp_centres: bool = False, plot_style: dict = None, gp: GaussianProcess = None)

Plot the variational posterior.

plot displays the variational posterior as a cornerplot showing the 1D and 2D marginals, estimated from samples. It uses the corner package.

plot also optionally displays the centres of the variational mixture components and the datapoints of the underlying Gaussian process (GP) used by VBMC. The plot can be enhanced by custom styles and specific datapoints of the GP can be highlighted.

Parameters:

n_samplesint, optional

The number of posterior samples used to create the plot, by default int(1e5).

titlestr, optional

The title of the plot, by default None.

plot_databool, optional

Whether to plot the datapoints of the GP, by default False.

highlight_datanp.ndarray, optional

Indices of the GP datapoints that should be plotted in a different way than the other datapoints, by default None.

plot_vp_centresbool, optional

Whether to plot the centres of the vp components, by default False.

plot_styledict, optional

A dictionary of plot styling options. The possible options are:

cornerdict, optional

Styling options directly passed to the corner function. By default: {"fig": plt.figure(figsize=(6, 6)), "labels": labels}. See the documentation of corner.

datadict, optional

Styling options used to plot the GP data. By default: {"s":15, "color":'blue', "facecolors": "none"}.

highlight_datadict, optional

Styling options used to plot the highlighted GP data. By default: {"s":15, "color":"orange"}.

vp_centredict, optional

Styling options used to plot the vp centres. By default: {"marker":"x", "color":"red"}.

gpgpyreg.GaussianProcess, optional

The GP to use for plotting actively sampled points.

Returns:

figmatplotlib.figure.Figure

The resulting matplotlib figure of the plot.

property rng: Generator

The random generator used by every method of this object that draws random numbers (sample, kl_div, mtv, plot, …). Copies made with copy.deepcopy share the same generator, so that all copies of a variational posterior draw from a single random stream.

sample(N: int, orig_flag: bool = True, balance_flag: bool = False, df: float = inf)

Draw random samples from the variational posterior.

Parameters:

Nint

Number of samples to draw, as any scalar that holds a whole number: a float such as 1e5, a NumPy scalar or a 0-D array is taken as that number of samples.

orig_flagbool, optional

If orig_flag is True, the random vectors are returned in the original parameter space. If False, they are returned in the transformed, unconstrained space used internally by VBMC. By default True.

balance_flagbool, optional

If balance_flag is True, the generating process is balanced such that the random samples come from each mixture component exactly proportionally (or as close as possible) to the variational mixture weights. If False, the generating mixture component for each sample is determined randomly, according to the mixture weights. By default False.

dffloat, optional

Generate the samples from a heavy-tailed version of the variational posterior, in which the multivariate normal components have been replaced by multivariate t-distributions with df degrees of freedom. The default is np.inf, limit in which the t-distribution becomes a multivariate normal. A finite negative df is refused: pdf reads it as the product of D univariate t densities, a family this method cannot draw from.

Returns:

Xnp.ndarray

X is an N-by-D matrix of random vectors drawn from the variational posterior.

Inp.ndarray

I is an N-element array of integers such that the i-th element of I indicates the index of the variational mixture component from which the i-th row of X has been generated.

Raises:

ValueError

Raised if N is not a scalar holding a whole number of samples.

ValueError

Raised if df is finite and negative (the product of univariate t densities has no sampler here).

Notes

Random draws use self.rng.

save(file, overwrite=False)

Save the VP to a file.

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

set_parameters(theta: ndarray, raw_flag=True)

Set variational posterior parameters from a single array.

Takes as input a numpy array and assigns it to the variational posterior parameters. eta, the unbounded (softmax) parametrization of the mixture weights, is set alongside w, so that softmax(eta) is the resulting w.

When both sigma and lambd are optimized, lambd is then divided by its root mean square and sigma multiplied by it. The two scales enter the density through their product alone, so that rescaling leaves the distribution as it was. A scale that is not optimized keeps its value and the other takes the value in theta, so that the parameters set are a function of theta alone. The weights are divided by their sum when they are being optimized.

Parameters:

thetanp.ndarray

The array with the parameters that should be assigned.

raw_flagbool, optional

Specifies whether the scales sigma and lambd and the weights are passed as raw (unconstrained) parameters or as they are, by default True. The raw scales are their logarithms, and the raw weights are the parameters of a softmax: the weights are their exponentials, normalized to sum to one.

Raises:

ValueError

Raised if raw_flag is False and a scale or a weight in theta is negative.

Export methods

The export methods have dedicated API pages:

• Export to a torch distribution

• Export to an ArviZ DataTree