modified_harmonic_mean#

cultivars.bayes.evidence.modified_harmonic_mean(draws, log_kernel, *, tau=_MHM_TAU, source='draws')[source]#

Geweke’s (1999) modified harmonic mean on posterior draws.

The harmonic-mean identity \(m(y)^{-1} = \mathbb{E}_{\theta \mid y}\bigl[g(\theta) / (f(y \mid \theta)\, \pi(\theta))\bigr]\) holds for any density \(g\) on the parameter space, and the plain harmonic mean, \(g = \pi\), has infinite variance because the prior’s tails are heavier than the posterior’s. Geweke’s modification takes \(g\) to be a Gaussian fitted to the draws and truncated to its central \(\tau\) ellipsoid, so that \(g\) has lighter tails than the posterior by construction:

\[\log \widehat{m}(y) = -\log \frac{1}{S} \sum_{s=1}^{S} \frac{g_\tau(\theta^{(s)})}{f(y \mid \theta^{(s)})\, \pi(\theta^{(s)})}, \qquad g_\tau(\theta) = \frac{\phi(\theta;\, \hat\mu, \hat\Sigma)}{\tau} \, \mathbb{1}\{(\theta - \hat\mu)^\top \hat\Sigma^{-1} (\theta - \hat\mu) \le \chi^2_{d}(\tau)\}.\]

The sum is formed on the log scale, and the Monte Carlo standard error of the log estimate is the standard deviation of the ratios divided by the square root of their Geyer effective sample size, so autocorrelated draws are charged for their autocorrelation. The share of draws that fall inside the ellipsoid is compared with \(\tau\): a share far from it says the posterior is not the Gaussian the envelope assumes, and the record notes it.

Parameters:

  • draws (ArrayLike) – (S,) or (S, d) posterior draws, finite throughout, with S greater than 2d + 2 so a Gaussian can be fitted. Transform bounded parameters to the real line first (a log for a variance, a logit for a probability) and include the Jacobian in log_kernel; a Gaussian envelope that leaves the support biases the estimate.

  • log_kernel (ArrayLike) – (S,) unnormalized log posterior, log likelihood plus log prior, evaluated at each draw in the same parameterization as draws.

  • tau (float) – Probability mass of the Gaussian envelope retained, in (0, 1). Draws outside the tau ellipsoid receive zero weight; a smaller value is more robust to non-Gaussian tails and uses fewer draws.

  • source (str) – Label for the record’s source field, which is the row label in compare().

Return type:

_MarginalLikelihoodSelection

Returns:

The MarginalLikelihoodSelection with method set to "modified harmonic mean", the log estimate, its Monte Carlo standard error (nan when fewer than eight draws fall inside the ellipsoid), n_draws the number of draws supplied, nobs zero since the function cannot know the sample, and a note on the envelope, plus a second note when the inside share differs from tau by more than a tenth.

Raises:

  • SpecificationError – If tau is not inside (0, 1), if there are too few draws to fit the envelope, or if no draw falls inside the truncation ellipsoid.

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

Note

The estimator is biased upward when the posterior’s tails are heavier than the envelope’s inside the ellipsoid, and the coverage check catches only gross failures; on a skewed posterior the share can sit near tau while the estimate is off by a tenth of a nat. For a roughly elliptical posterior bridge_sampling() has lower variance and no such bias, at the price of a callable kernel. Records from this function carry nobs = 0, so compare() cannot check that they score the same sample as the records they are ranked against; that is the caller’s responsibility.

See also

  • bridge_sampling() – the lower-variance estimator when the kernel can be evaluated at new points.

  • marginal_likelihood() – the same estimator applied by the particle-chain DSGE result to its own draws.

  • compare() – ranks the records.

References

Geweke, J. (1999). Using simulation methods for Bayesian econometric models: Inference, development, and communication. Econometric Reviews, 18(1), 1-73.

Newton, M. A., & Raftery, A. E. (1994). Approximate Bayesian inference with the weighted likelihood bootstrap. Journal of the Royal Statistical Society B, 56(1), 3-48.

Example

Draws from a standard normal in three dimensions with the unnormalized kernel \(-\tfrac{1}{2}\|z\|^2\); the evidence is the missing constant \((2\pi)^{3/2}\):

>>> import numpy as np
>>> rng = np.random.default_rng(1)
>>> z = rng.standard_normal((3000, 3))
>>> record = modified_harmonic_mean(z, -0.5 * (z**2).sum(axis=1))
>>> bool(abs(record.log_value - 1.5 * np.log(2 * np.pi)) < 0.05)
True
>>> record.method, record.n_draws
('modified harmonic mean', 3000)

A bounded parameter left on its own scale biases the estimate; the log transform with its Jacobian removes the bias. Here the draws are exponential, the kernel is \(-x\), and the true log evidence is zero:

>>> rng = np.random.default_rng(3)
>>> x = rng.exponential(size=4000)
>>> raw = modified_harmonic_mean(x, -x)
>>> u = np.log(x)
>>> transformed = modified_harmonic_mean(u, -np.exp(u) + u)
>>> round(raw.log_value, 2), round(transformed.log_value, 2)
(0.12, 0.01)