Marginal likelihood
In Bayesian statistics, the marginal likelihood, also known as the Bayesian evidence or model evidence, is the probability of observing the data given a model, obtained by integrating the joint likelihood of the data and model parameters with respect to the prior distribution over the parameters.[1] Mathematically, for data D and model M_k with parameters h_k, it is expressed as p(D \mid M_k) = \int p(D \mid h_k, M_k) \, p(h_k \mid M_k) \, dh_k.[2] This integral marginalizes out the parameters, yielding a quantity that depends only on the data and the model structure, including the prior.[1] The marginal likelihood plays a central role in Bayesian model comparison and selection, as it quantifies how well a model predicts the data while automatically incorporating a complexity penalty through the prior's influence on the integration.[2] Specifically, the Bayes factor, which is the ratio of the marginal likelihoods of two competing models, serves as a measure of evidence in favor of one model over the other, updating prior odds to posterior odds without requiring nested models.[3] It is essential for Bayesian model averaging, where posterior model probabilities are computed as p(M_k \mid D) = \frac{p(D \mid M_k) p(M_k)}{\sum_l p(D \mid M_l) p(M_l)}, allowing inference that accounts for model uncertainty.[1] Applications span fields like hydrology,[1] phylogenetics,[4] and cosmology,[5] where it aids in hypothesis testing and hyperparameter optimization. Computing the marginal likelihood exactly is often intractable for complex models due to the high-dimensional integral, necessitating approximations such as Laplace's method, thermodynamic integration, or Markov chain Monte Carlo (MCMC) estimators like the harmonic mean.[1] These methods balance accuracy and computational feasibility, with thermodynamic integration noted for its consistency in environmental modeling contexts.[1] The sensitivity to prior choices underscores its interpretive challenges, as it encodes Occam's razor by favoring parsimonious models that fit the data well.[2]Definition and Basics
Formal Definition
The marginal likelihood, denoted as p(y), is formally defined as the marginal probability of the observed data y, obtained by integrating the joint distribution of the data and the model parameters \theta over the parameter space: p(y) = \int p(y, \theta) \, d\theta = \int p(y \mid \theta) p(\theta) \, d\theta, where p(y \mid \theta) is the likelihood function and p(\theta) is the prior distribution over the parameters \theta.[6] This integral represents the prior predictive distribution of the data, averaging the likelihood across all possible parameter values weighted by the prior.[6] The process of marginalization in this context involves integrating out the parameters \theta, which are treated as nuisance parameters, to yield a probability for the data y that is independent of any specific parameter values.[6] This eliminates dependence on \theta by averaging over its uncertainty as specified by the prior, resulting in a quantity that summarizes the model's predictive content for the observed data without conditioning on particular estimates of \theta.[6] In contrast to the joint probability p(y, \theta) = p(y \mid \theta) p(\theta), which depends on both the data and parameters, the marginal likelihood p(y) marginalizes away \theta and thus provides the unconditional probability of the data under the model.[6] This distinction underscores how marginalization transforms the parameter-dependent joint into a data-only marginal.[6] The term "marginal likelihood" is used interchangeably with "marginal probability" or "evidence" to refer to p(y), particularly in contexts emphasizing its role as a model-level summary.[7]Bayesian Interpretation
In Bayesian statistics, the marginal likelihood represents the predictive probability of the observed data under a given model, obtained by averaging the likelihood over all possible parameter values weighted by their prior distribution. This integration marginalizes out the parameters, providing a coherent measure of how well the model explains the data while incorporating prior beliefs about parameter uncertainty. As such, it serves as the Bayesian evidence for the model, distinct from conditional measures that fix parameters at specific values.[8] The marginal likelihood functions as a key indicator of model plausibility, where a higher value suggests the model offers a better overall fit to the data after accounting for the full range of parameter uncertainty encoded in the prior. Unlike point estimates, this averaging process naturally balances goodness-of-fit with the model's capacity to generalize, implicitly favoring parsimonious models that do not overfit by spreading probability mass too thinly across parameters. In model comparison, it enables direct assessment of competing hypotheses by quantifying the relative support each model receives from the data alone.[8] In contrast to maximum likelihood estimation, which maximizes the likelihood at a single point estimate of the parameters and thus relies on plug-in predictions that ignore uncertainty, the marginal likelihood integrates over the entire parameter space. This integration implicitly penalizes model complexity through the prior's influence, as overly flexible models tend to dilute the evidence by assigning low density to the data under broad parameter ranges, whereas maximum likelihood can favor intricate models that fit noise without such restraint. Consequently, marginal likelihood promotes more robust inference in finite samples by embedding Occam's razor directly into the evaluation.[2] Historically, the marginal likelihood emerged as a central component in Bayesian model assessment through its role in computing Bayes factors, which compare the evidence for rival models and update prior odds to posterior odds in a coherent manner. This framework, formalized in seminal work on Bayes factors, provided a principled alternative to frequentist testing for evaluating scientific theories.[8]Mathematical Framework
Expression in Parametric Models
In parametric statistical models, the marginal likelihood integrates out the model parameters to obtain the probability of the observed data under the model specification. Consider a parametric model M indexed by a parameter vector \theta \in \Theta, where \Theta is the parameter space. The marginal likelihood is then expressed as p(y \mid M) = \int_{\Theta} p(y \mid \theta, M) \, \pi(\theta \mid M) \, d\theta, where p(y \mid \theta, M) denotes the likelihood function and \pi(\theta \mid M) is the prior distribution over the parameters given the model.[8] This formulation arises naturally in Bayesian inference as the normalizing constant for the posterior distribution.[8] When the parameter space \Theta is discrete, the continuous integral is replaced by a summation: p(y \mid M) = \sum_{\theta \in \Theta} p(y \mid \theta, M) \, \pi(\theta \mid M). This discrete form is applicable in scenarios such as finite mixture models with a discrete number of components or models with categorical parameters.[8] The dimensionality of the integral or sum corresponds directly to the dimension of \theta, reflecting the number of parameters being marginalized out; higher dimensionality increases computational demands but maintains the core structure of the expression.[8] A concrete example illustrates this in a simple univariate normal model, where the observed data y follows y \sim \mathcal{N}(\mu, \sigma^2) with known \sigma^2, and the prior on the mean is \mu \sim \mathcal{N}(0, 1). The marginal likelihood involves the integral p(y \mid M) = \int_{-\infty}^{\infty} \mathcal{N}(y \mid \mu, \sigma^2) \, \mathcal{N}(\mu \mid 0, 1) \, d\mu. Completing the square in the exponent yields a closed-form normal distribution for the marginal: p(y \mid M) = \mathcal{N}(y \mid 0, \sigma^2 + 1). However, in the more general univariate normal model with unknown variance and a conjugate normal-inverse-gamma prior on (\mu, \sigma^2), the full marginal likelihood over both parameters results in a closed-form Student's t-distribution for y, highlighting how prior choices enable analytical tractability.Relation to Posterior and Likelihood
In Bayesian inference, the marginal likelihood plays a central role in Bayes' theorem, which expresses the posterior distribution asp(\theta \mid y) = \frac{p(y \mid \theta) p(\theta)}{p(y)},
where p(y \mid \theta) is the likelihood function, p(\theta) is the prior distribution, and p(y) is the marginal likelihood of the data y. This formulation updates prior beliefs about the parameters \theta with observed data to yield the posterior p(\theta \mid y).[6] The marginal likelihood p(y) serves as the normalizing constant—or evidence—that ensures the posterior integrates to 1 over the parameter space, transforming the unnormalized product p(y \mid \theta) p(\theta) into a proper probability density. Computationally, it is obtained by integrating the joint density of data and parameters: p(y) = \int p(y \mid \theta) p(\theta) \, d\theta. This integration averages the likelihood across all possible parameter values weighted by the prior, providing a data-dependent measure of model plausibility independent of specific \theta.[6] Unlike the conditional likelihood p(y \mid \theta), which conditions on fixed parameters and evaluates model fit for particular \theta, the marginal likelihood marginalizes the full likelihood over the prior distribution, incorporating parameter uncertainty from the outset. This contrasts with the profiled likelihood in frequentist settings, where nuisance parameters are eliminated by maximization rather than integration, leading to a point-estimate-based adjustment without prior weighting.[6][9] The marginal likelihood thus quantifies the total predictive uncertainty in the data under the model, encompassing both the variability in the likelihood due to unknown parameters and the prior's influence on that variability, whereas the likelihood alone fixes \theta and ignores such averaging. This broader uncertainty measure supports coherent probabilistic reasoning in Bayesian frameworks.[6]
Computation Methods
Analytical Approaches
Analytical approaches to computing the marginal likelihood rely on cases where the integral over the parameter space can be evaluated exactly in closed form, which occurs primarily under the use of conjugate priors.[10] Conjugate priors are distributions from the same family as the posterior, allowing the marginal likelihood to be derived as a normalizing constant without numerical integration.[11] This solvability holds when the likelihood and prior combine to yield a posterior in the same parametric family, simplifying the integration to known functions like the beta or gamma integrals.[10] A classic example is the beta-binomial model, where the binomial likelihood for coin flips is paired with a beta prior on the success probability \theta. The marginal likelihood is then given by p(y \mid n, \alpha, \beta) = \binom{n}{y} \frac{B(\alpha + y, \beta + n - y)}{B(\alpha, \beta)}, where B denotes the beta function, representing the integral over \theta.[12] This closed-form expression arises directly from the conjugacy, enabling exact Bayesian inference for binary data.[13] In Gaussian models, conjugate priors such as the normal-inverse-Wishart facilitate analytical marginal likelihoods, often resulting in a multivariate Student's t-distribution for the data. For a multivariate normal likelihood with unknown mean \mu and precision \Lambda, the marginal distribution of the data is p(D) = \pi^{-nd/2} \frac{\Gamma_d(\nu_n/2)}{\Gamma_d(\nu_0/2)} \frac{|\Lambda_0|^{\nu_0/2}}{|\Lambda_n|^{\nu_n/2}} \left( \frac{\kappa_0}{\kappa_n} \right)^{d/2}, where \nu_n = \nu_0 + n, \kappa_n = \kappa_0 + n, \Lambda_n = \Lambda_0 + S + \frac{\kappa_0 n}{\kappa_n} (\bar{x} - \mu_0)(\bar{x} - \mu_0)^T , S = \sum_{i=1}^n (x_i - \bar{x})(x_i - \bar{x})^T , d is the dimensionality, and n is the number of observations.[10] This highlights the tractability in linear Gaussian settings. However, such analytical solutions are rare, particularly in high-dimensional settings or with non-conjugate priors, where the parameter integral becomes intractable and necessitates numerical methods.[14] These limitations stem from the exponential growth in integration complexity as dimensionality increases, restricting exact computations to low-dimensional or specially structured models.[15]Numerical and Approximation Techniques
When exact analytical computation of the marginal likelihood is infeasible, such as in non-conjugate models with high-dimensional parameter spaces, numerical and approximation techniques become essential for estimation. These methods leverage sampling or asymptotic expansions to approximate the integral \int p(y \mid \theta) p(\theta) \, d\theta, balancing computational feasibility with accuracy. Monte Carlo-based approaches, in particular, provide unbiased or consistent estimators but often require careful tuning to manage variance, while deterministic methods offer faster but potentially biased approximations suitable for large-scale applications. Monte Carlo methods, including importance sampling, form a foundational class of estimators for the marginal likelihood. The importance sampling estimator approximates the integral as \hat{p}(y) \approx \frac{1}{N} \sum_{i=1}^N \frac{p(y \mid \theta_i) p(\theta_i)}{q(\theta_i)}, where \{\theta_i\}_{i=1}^N are samples drawn from a proposal distribution q(\theta) chosen to approximate the posterior p(\theta \mid y). This self-normalized form ensures consistency under mild conditions on q, though high variance arises if q poorly covers the posterior support, necessitating techniques like adaptive proposals or multiple importance sampling to improve efficiency in complex models.[16] Markov Chain Monte Carlo (MCMC) methods extend these ideas by generating dependent samples from the posterior to estimate the marginal likelihood without direct prior sampling. Bridge sampling, for instance, uses samples from both the prior and posterior to compute ratios via an optimal bridge function, yielding \hat{p}(y) = \frac{1}{N} \sum_{i=1}^N \frac{p(y \mid \theta_i^\text{prior}) p(\theta_i^\text{prior})}{g(\theta_i^\text{prior} \mid y)} \cdot \frac{1}{M} \sum_{j=1}^M \frac{g(\theta_j^\text{post} \mid y)}{p(y \mid \theta_j^\text{post}) p(\theta_j^\text{post})}, where g minimizes estimation variance; this approach is particularly robust for multimodal posteriors. The harmonic mean estimator, another posterior-based method, approximates p(y) as the reciprocal of the average inverse likelihood over posterior samples, \hat{p}(y) \approx \left( \frac{1}{N} \sum_{i=1}^N \frac{1}{p(y \mid \theta_i)} \right)^{-1}, but suffers from infinite variance in heavy-tailed cases, prompting stabilized variants.[17][18] Deterministic approximations, such as the Laplace method, provide closed-form estimates by exploiting local behavior around the posterior mode. This technique approximates the log-posterior as quadratic via its Hessian matrix H at the mode \theta^*, leading to p(y) \approx p(y \mid \theta^*) p(\theta^*) (2\pi)^{d/2} |H|^{-1/2}, where d is the parameter dimension; the approximation improves asymptotically as n \to \infty for large samples but can bias results in small-data or highly nonlinear settings. It is computationally efficient, requiring only optimization and Hessian evaluation, and serves as a building block for higher-order corrections in moderate dimensions.[19] Variational inference offers a scalable lower bound on the log-marginal likelihood through optimization, framing estimation as minimizing the Kullback-Leibler divergence between a tractable variational distribution q(\theta) and the true posterior. The evidence lower bound (ELBO) states \log p(y) \geq \mathbb{E}_q [\log p(y, \theta)] - \mathbb{E}_q [\log q(\theta)], maximized by adjusting q (often mean-field or structured forms) via stochastic gradients; this bound is tight when q \approx p(\theta \mid y) and enables fast inference in massive datasets, though it underestimates the true value and requires careful family selection to avoid loose bounds.[20] Recent advancements, including those post-2020, have enhanced these techniques for scalability in deep learning models, where parameter spaces exceed millions of dimensions. Annealed importance sampling (AIS) refines importance sampling by introducing intermediate distributions bridging prior and posterior, with differentiable variants enabling end-to-end optimization of annealing schedules for tighter estimates in generative models.[21][22] Sequential Monte Carlo (SMC) methods, propagating particles through annealed sequences with resampling, have been adapted for deep Bayesian networks, achieving unbiased marginal likelihoods via thermodynamic integration while integrating neural proposals for efficiency in high-dimensional tasks like variational autoencoders.[23] More recent methods, such as generalized stepping-stone sampling (2024), improve efficiency in specific domains like pulsar timing analysis.[24] These developments prioritize variance reduction and GPU acceleration, facilitating model comparison in neural architectures.Applications in Statistics
Model Comparison and Selection
One key application of the marginal likelihood in Bayesian statistics is model comparison, where it serves as the basis for the Bayes factor, a measure of the relative evidence provided by the data for two competing models. The Bayes factor BF_{12} comparing model M_1 to model M_2 is defined as the ratio of their marginal likelihoods:BF_{12} = \frac{p(y \mid M_1)}{p(y \mid M_2)},
where y denotes the observed data.[25] This ratio quantifies how much more likely the data are under M_1 than under M_2, after integrating out model parameters via their priors, thereby providing a coherent framework for hypothesis testing and model selection without relying on arbitrary significance thresholds.[25] For nested models, where M_1 is a special case of M_2 (e.g., by imposing parameter restrictions), the Savage-Dickey density ratio offers a convenient way to compute the Bayes factor using posterior and prior densities evaluated at the boundary values of the restricted parameters. Specifically, the Bayes factor is given by
BF_{12} = \frac{p(\theta \mid y, M_2)}{\pi(\theta \mid M_2)},
evaluated at the values of \theta (the restricted parameter) that define the nesting boundary, under the encompassing model M_2, assuming the prior under M_2 matches the prior under M_1 for the common parameters. This approach simplifies computation by avoiding full marginal likelihood estimation for both models, though it requires careful prior specification to ensure validity.[26] The marginal likelihood inherently implements Occam's razor in model selection by favoring parsimonious models that adequately fit the data, as more complex models must allocate prior probability mass across larger parameter spaces, effectively penalizing overparameterization unless the data strongly support the added complexity.[25] In practice, Bayes factors are interpreted using guidelines such as those proposed by Kass and Raftery, where values between 1 and 3 provide barely worth mentioning evidence, 3 to 20 indicate positive evidence, 20 to 150 strong evidence, and greater than 150 very strong evidence in favor of the numerator model (e.g., BF_{12} > 10 suggests substantial support for M_1).[25] However, Bayes factors exhibit sensitivity to the choice of priors, which can lead to varying conclusions if priors are not chosen judiciously, and their computation becomes prohibitively expensive for high-dimensional or large-scale models, often necessitating approximations like those from MCMC methods.[25]