bridge_sampling#

cultivars.bayes.evidence.bridge_sampling(draws, log_kernel, kernel, *, seed=None, n_proposal=None, source='draws')[source]#

Meng-Wong bridge sampling with a Gaussian proposal fitted to the draws.

Bridge sampling estimates the ratio of two normalizing constants from draws of both densities. With the unnormalized posterior \(p(\theta) = f(y \mid \theta)\, \pi(\theta)\) and a normalized proposal \(q(\theta)\), for any bridge function \(\alpha\),

\[m(y) = \frac{\mathbb{E}_{q}\bigl[p(\theta)\, \alpha(\theta)\bigr]} {\mathbb{E}_{p/m}\bigl[q(\theta)\, \alpha(\theta)\bigr]},\]

and Meng and Wong’s optimal bridge, \(\alpha \propto 1 / (s_1 p + s_2\, m\, q)\) with \(s_1, s_2\) the proposal and posterior shares of the draws, depends on \(m\) itself, so the estimate is the fixed point of an iteration started at the median log ratio. The draws are split in half: the first half fits the Gaussian proposal, the second half enters the bridge, so that the proposal is not fitted to the draws it is then compared with. Because the proposal’s draws must be scored by the posterior, the kernel has to be a callable and not a vector of values, and it has to return a finite number wherever a Gaussian around the draws can land.

On a roughly elliptical posterior the estimate has lower variance than modified_harmonic_mean() and none of its truncation bias. The reported error is the Frühwirth-Schnatter (2004) approximation, which treats the proposal draws as independent and charges the posterior half for its autocorrelation through its Geyer effective size.

Parameters:

  • draws (ArrayLike) – (S,) or (S, d) posterior draws, finite throughout, with S greater than 2d + 2. Bounded parameters must be transformed to the real line with the Jacobian in the kernel, as for modified_harmonic_mean(); here the requirement is strict rather than a matter of bias, because the kernel is evaluated at proposal points that may fall outside a bounded support.

  • log_kernel (ArrayLike) – (S,) unnormalized log posterior at each draw.

  • kernel (Callable[[NDArray[double]], float]) – The same log posterior as a function of one (d,) point, returning a float. It is called once per proposal draw.

  • seed (int | Generator | None) – Seed or generator for the proposal draws. An integer gives a reproducible estimate; a numpy.random.Generator is used in place and advanced.

  • n_proposal (int | None) – Proposal draws, at least eight. None uses as many as enter the bridge from the posterior side, S // 2.

  • source (str) – Label for the record’s source field.

Return type:

_MarginalLikelihoodSelection

Returns:

The MarginalLikelihoodSelection with method set to "bridge sampling", the log estimate, its approximate Monte Carlo standard error, n_draws the number of posterior draws supplied, nobs zero, and a note recording the proposal and the iterations to convergence.

Raises:

  • SpecificationError – If there are too few draws to fit the proposal, or n_proposal is below eight.

  • DimensionError – If draws is not one- or two-dimensional, or log_kernel does not have one value per draw.

  • NumericalError – If any draw or kernel value is not finite, if the kernel returns a non-finite value at a proposal point, or if the fixed-point iteration does not converge within a thousand steps.

Note

The kernel is evaluated n_proposal times, once per proposal draw; for a model whose likelihood is a Kalman filter that is S // 2 filter passes by default. Records from this function carry nobs = 0, so compare() cannot check that they score the same sample as the records they are ranked against.

See also

References

Meng, X.-L., & Wong, W. H. (1996). Simulating ratios of normalizing constants via a simple identity: A theoretical exploration. Statistica Sinica, 6(4), 831-860.

Frühwirth-Schnatter, S. (2004). Estimating marginal likelihoods for mixture and Markov switching models using bridge sampling techniques. Econometrics Journal, 7(1), 143-167.

Gronau, Q. F., Sarafoglou, A., Matzke, D., Ly, A., Boehm, U., Marsman, M., Leslie, D. S., Forster, J. J., Wagenmakers, E.-J., & Steingroever, H. (2017). A tutorial on bridge sampling. Journal of Mathematical Psychology, 81, 80-97.

Example

A bivariate standard normal, kernel \(-\tfrac{1}{2}\|z\|^2\), whose evidence is \(2\pi\):

>>> import numpy as np
>>> rng = np.random.default_rng(2)
>>> z = rng.standard_normal((3000, 2))
>>> record = bridge_sampling(
...     z, -0.5 * (z**2).sum(axis=1), lambda x: -0.5 * float(x @ x), seed=0
... )
>>> bool(abs(record.log_value - np.log(2 * np.pi)) < 0.05)
True
>>> record.method, record.mcse < 0.01
('bridge sampling', True)
>>> record.notes[0]
'Bridge sampling (Meng & Wong, 1996), ... converged in ... iterations.'
On this page