cultivars.bayes.evidence#
Marginal likelihoods, Bayes factors, and posterior model probabilities.
The marginal likelihood
is the probability a model assigned to the data before seeing them, and
it is the Bayesian answer to “which model”: the ratio \(m_1(y) /
m_2(y)\) is the Bayes factor, and with prior model probabilities it gives
posterior model probabilities. Its logarithm is the one number a
posterior can report in place of a log likelihood, and unlike a
likelihood it carries its own penalty for complexity, since a model that
spreads prior mass over parameters the data do not want pays for it in
\(m(y)\). Three routes to it live in the package, each on the results
that admit it. The conjugate
BVAR integrates its
Normal-inverse-Wishart posterior in closed form, so its number is exact.
The GibbsBVAR under a
static prior uses Chib’s (1995) identity, \(\log m(y) = \log f(y
\mid \theta^*) + \log \pi(\theta^*) - \log \pi(\theta^* \mid y)\), with
the posterior ordinate at a high-density point estimated from the
two-block Gibbs output and no extra runs. The particle-chain
PerturbationDSGE
posterior uses Geweke’s (1999) modified harmonic mean, and since its
likelihood is itself a particle-filter estimate, the number inherits that
noise.
Two commitments shape the surface. First, every route returns the same
record with its Monte Carlo standard error, zero when exact, so that
compare() can say when two models are not separated by the
estimates at hand, and so that
bayesian_model_average() can warn that
a weight ratio is only as precise as the difference of two noisy
estimates. Second, what is not implemented is refused with the reason
rather than approximated silently: the plain harmonic mean, whose
variance is infinite because the prior’s tails are heavier than the
posterior’s; Chib’s estimator for the Gibbs BVAR under an adaptive prior,
which needs reduced runs; and Chib’s estimator with dummy observations,
where the augmented sample makes the ordinate’s meaning subtle. Sampled
results without a route – the stochastic-volatility family, the
Student-t and stochastic-volatility BVARs, the time-varying-parameter
models – say so by name and point to out-of-sample predictive score.
Layout. marginal_likelihood() is the dispatcher: it calls a result’s
own marginal_likelihood() and returns the record, which the three
results above define on themselves and which delegates to _core –
the closed-form update for the conjugate model, _core._estimators.
_chib_independent_normal_wishart for the Gibbs model, and
_modified_harmonic_mean for the DSGE posterior.
modified_harmonic_mean() and bridge_sampling() expose the two
generic estimators on draws that came from elsewhere, another sampler or
a kernel of the user’s own, through _core._estimators.
_modified_harmonic_mean and _bridge_sampling, both of which fit
their Gaussian envelope with _gaussian_envelope and validate their
input with _validate_posterior_draws. compare() ranks any set
of records by delegating to the record’s own compare method in
_internals._selections, which builds the table and reads the
Kass-Raftery evidence scale through _core._converters._evidence_label.
The record itself,
MarginalLikelihoodSelection, is the
_internals dataclass every route constructs.
References
Chib, S. (1995). Marginal likelihood from the Gibbs output. Journal of the American Statistical Association, 90(432), 1313-1321.
Geweke, J. (1999). Using simulation methods for Bayesian econometric models: Inference, development, and communication. Econometric Reviews, 18(1), 1-73.
Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773-795.
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.
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
The evidence of a bivariate standard normal from its unnormalized kernel \(-\tfrac{1}{2}\|z\|^2\) is the missing constant \(2\pi\):
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> z = rng.standard_normal((4000, 2))
>>> record = modified_harmonic_mean(z, -0.5 * (z**2).sum(axis=1), source="N(0, I)")
>>> round(record.log_value, 1) # log(2 pi)
1.8
>>> record.method
'modified harmonic mean'
A package result carries its own route, and the records rank:
>>> from cultivars.multivariate.large_dim.bayesian import BVAR
>>> y = np.zeros((120, 2))
>>> for t in range(1, 120):
... y[t] = 0.7 * y[t - 1] + rng.standard_normal(2)
>>> one = marginal_likelihood(BVAR(y, order=1).fit(n_draws=100, seed=0))
>>> one.method, one.mcse
('analytic', 0.0)