Bayesian Inference

Bayesian inference is a formal framework of statistical deduction and inductive learning in which probability is interpreted as a quantitative measure of epistemic uncertainty about unknown parameters, latent states, or competing scientific hypotheses. Named after Thomas Bayes (whose 1763 essay was edited and published posthumously by Richard Price) and developed into its general continuous analytical form by Pierre-Simon Laplace in 1774, Bayesian inference treats an unknown parameter vector θ ∈ Θ as a random variable governed by a prior probability distribution p(θ) that encodes domain knowledge or physical constraints before new observations are collected. Once empirical data D are observed, the prior is combined with the sampling likelihood function p(D | θ) through Bayes' theorem to yield the full posterior probability distribution p(θ | D).

Unlike frequentist estimation—which treats θ as a fixed constant and evaluates procedures over hypothetical repeated samples from the sample space (as in maximum-likelihood phylogenetics)—Bayesian inference conditions directly on the observed dataset D and propagates parameter uncertainty across all nuisance variables via probability integration. In modern computational biology, Bayesian methods underpin divergence-time dating under relaxed molecular clocks, species-tree estimation under the multispecies coalescent, positive selection scanning via Bayes Empirical Bayes in codon substitution models, and high-dimensional posterior sampling using Markov chain Monte Carlo (MCMC) algorithms.

Mathematical Formulation and the Marginal Likelihood (Model Evidence)

For a continuous parameter vector θ ∈ Θ and observed dataset D under an assumed model structure M, Bayes' theorem expresses the posterior density p(θ | D, M) as the normalized product of the likelihood and the prior:

p(θ | D, M) = [ p(D | θ, M) p(θ | M) ] / p(D | M)

Each component of Bayes' equation plays a distinct mathematical role:

• Prior Distribution p(θ** | M): Specifies the initial probability density over the parameter space before observing D** (for example, an exponential distribution p(t) = λ exp(−λt) over phylogenetic branch lengths t ≥ 0, or a Dirichlet distribution over stationary nucleotide frequencies π).
• Likelihood Function p(D** | θ, M): Evaluates the probability density of generating the observed data D** given parameter values θ, computed in molecular systematics from a sequence alignment via Felsenstein's tree-pruning algorithm.
• Marginal Likelihood (Model Evidence) p(D** | M)**: The normalizing denominator obtained by integrating the joint density over the entire d-dimensional parameter space Θ:

p(D | M) = ∫_Θ p(D | θ, M) p(θ | M) dθ

• Posterior Distribution p(θ** | D, M): The updated probability distribution over θ** given the data D, from which 95% Highest Posterior Density (HPD) credible intervals are extracted as the shortest parameter interval containing 95% of the integrated posterior probability mass.

Conjugate Priors, Hierarchical Models, and Maximum A Posteriori (MAP) Estimation

When the prior distribution p(θ) and the sampling likelihood p(D | θ) belong to matched algebraic families—typically within the exponential family of probability distributions—the resulting posterior distribution p(θ | D) falls into the same functional family as the prior with updated hyperparameters. Such priors are called conjugate priors and permit closed-form analytical posterior updating without numerical integration:

• Beta-Binomial Conjugacy: If the allele frequency θ ∈ [0, 1] has a Beta(α, β) prior density proportional to θ^α−1(1 − θ)^β−1 and k mutant alleles are observed in n independent trials (Binomial likelihood), the posterior is Beta(α + k, β + n − k), where α and β act as prior "pseudocounts."
• Dirichlet-Multinomial Conjugacy: Generalizes Beta-Binomial conjugacy to K-state categorical probabilities π = (π₁, ..., π_K) on the simplex ∑_i=1^K π_i = 1; a Dirichlet(α) prior combined with observed state counts c = (c₁, ..., c_K) yields a Dirichlet(α + c) posterior, widely used in sequence profile estimation and admixture modeling.
• Gamma-Poisson and Normal-Normal Conjugacy: A Gamma(α, β) prior on a Poisson substitution rate λ updates to Gamma(α + ∑ x_i, β + n), while a Gaussian prior on a Normal mean μ with known variance σ² yields a precision-weighted Gaussian posterior.

In hierarchical Bayesian models (and empirical Bayes approximations such as Bayes Empirical Bayes), the hyperparameters φ of the prior p(θ | φ) are themselves assigned hyperpriors p(φ) and estimated jointly across thousands of genes or alignment sites, shrinking noisy site-specific estimates toward the genome-wide mean. When a single point summary of the posterior is required, Maximum A Posteriori (MAP) estimation maximizes the unnormalized log-posterior objective:

θ_MAP = arg max_θ [ ln p(D | θ) + ln p(θ) ]

A zero-mean Gaussian prior on θ makes MAP estimation mathematically equivalent to L₂ (Ridge / Tikhonov) regularization, whereas a Laplace double-exponential prior exp(−λ‖θ‖₁) yields sparse L₁ (LASSO) regularization.

Model Selection via Marginal Likelihoods and Bayes Factors

To compare two competing evolutionary or physical models M₁ and M₂ (for example, a strict molecular clock vs. an uncorrelated lognormal relaxed clock, or a single-rate vs. partitioned substitution model), Bayesian inference evaluates the posterior odds ratio:

[ P(M₁ | D) / P(M₂ | D) ] = [ p(D | M₁) / p(D | M₂) ] × [ P(M₁) / P(M₂) ]

The ratio of marginal likelihoods, BF₁₂ = p(D** | M₁) / p(D | M₂), is the Bayes factor, which quantifies how strongly the observed data update prior model odds into posterior model odds. Because the marginal likelihood p(D** | M_i) = ∫ p(D | θ_i, M_i) p(θ_i | M_i) dθ_i averages the likelihood across the entire prior volume of M_i, Bayes factors incorporate an automatic Occam's razor: an overly complex model with many redundant parameters spreads its prior mass thinly over a vast parameter space, reducing p(D | M) unless those extra parameters substantially improve the fit to the data. Under Harold Jeffreys' logarithmic scale, 2 ln(BF₁₂) ∈ (2, 6) indicates positive evidence for M₁, (6, 10) indicates strong evidence, and 2 ln(BF₁₂) > 10 (BF₁₂ > 148.4) indicates decisive evidence. Accurate marginal likelihood estimation requires thermodynamic integration (path sampling) or stepping-stone sampling (Xie et al., 2011) along a power-posterior path p_β(θ | D) ∝ p(D | θ)^β p(θ) bridging the prior (β = 0) and the posterior (β = 1).

Metropolis-Hastings, Gibbs Sampling, and Parallel Tempering (MC³)

In phylogenetic inference over tree topologies τ and continuous parameters θ, the normalizing denominator p(D) sums over (2N − 5)!! trees and integrates over 2N − 3 branch lengths, making direct calculation impossible. In 1953, Nicholas Metropolis and colleagues—generalized in 1970 by W. Keith Hastings—showed that an ergodic Markov chain can be constructed whose stationary distribution equals the target posterior p(ψ | D) (where ψ = (τ, t, θ)) without ever computing p(D). At iteration t, a candidate state ψ′ is drawn from a proposal density q(ψ′ | ψ) (such as a Subtree Pruning and Regrafting move on τ or a sliding-window multiplier on branch lengths) and accepted with the Metropolis-Hastings acceptance probability:

α(ψ, ψ′) = min { 1, [ p(D | ψ′) p(ψ′) / (p(D | ψ) p(ψ)) ] × [ q(ψ | ψ′) / q(ψ′ | ψ) ] }

Because p(D) appears in both the numerator and denominator of p(ψ′ | D) / p(ψ | D), it cancels out. A special case is the Gibbs sampler (Geman and Geman, 1984), which updates one parameter block ψ_j at a time by sampling directly from its full conditional distribution p(ψ_j | ψ_−j, D), resulting in an acceptance probability of α = 1.

When the posterior landscape contains multiple isolated "tree islands" separated by deep valleys of low probability, a single MCMC chain can remain trapped in a local mode. Software packages such as MrBayes 3.2 and RevBayes overcome this barrier using Metropolis-coupled MCMC (MC³, or parallel tempering), which runs M parallel chains heated to inverse temperatures 1 = β₁ > β₂ > ... > β_M > 0, where chain m samples proportional to [p(ψ | D)]^β_m. Raising the posterior to a fractional power β_m < 1 flattens peaks and raises valleys, allowing heated chains to traverse topological barriers rapidly and propose state swaps with the cold chain (β₁ = 1), which alone is recorded after discarding initial burn-in iterations and verifying convergence via the Effective Sample Size (ESS ≥ 200) and Gelman-Rubin potential scale reduction factor (R̂ < 1.01).

How does a Bayesian 95% credible interval differ from a frequentist 95% confidence interval?

A Bayesian 95% Highest Posterior Density (HPD) credible interval conditions directly on the observed dataset D and means there is a 0.95 posterior probability that the true parameter θ lies within the stated bounds given the model and prior. A frequentist 95% confidence interval treats θ as fixed and means that if the data-generating experiment were repeated infinitely many times, 95% of the intervals constructed by that procedure would cover the fixed parameter value.

Why is the harmonic mean estimator of the marginal likelihood discredited for Bayes factors?

The harmonic mean estimator approximates 1 / p(D | M) by averaging 1 / p(D | θ^(t)) across posterior MCMC samples θ^(t). Because the reciprocal likelihood explodes in the low-likelihood tails of the prior where the posterior MCMC chain rarely visits, the estimator has infinite variance and systematically overestimates the marginal likelihood of parameter-rich models. Stepping-stone sampling (Xie et al., 2011) eliminates this bias by bridging the prior and posterior across β-power rungs.

Why can Bayesian clade posterior probabilities exceed non-parametric bootstrap proportions?

Bayesian posterior probabilities measure the probability that a clade is true conditional on the observed alignment, prior, and substitution model being correctly specified, whereas Felsenstein's non-parametric bootstrap measures sampling repeatability across resampled alignment columns. When the substitution model is moderately misspecified on long alignments, Bayesian posterior probabilities can converge rapidly toward 1.0 for an incorrect topology, making model adequacy testing and comparison with maximum-likelihood phylogenetics UFBoot values essential.

1. Hastings, W. K. (1970). "Monte Carlo sampling methods using Markov chains and their applications." Biometrika, 57(1), 97–109. https://doi.org/10.1093/biomet/57.1.97
2. Huelsenbeck, J. P., Ronquist, F., Nielsen, R., & Bollback, J. P. (2001). "Bayesian inference of phylogeny and its impact on evolutionary biology." Science, 294(5550), 2310–2314. https://doi.org/10.1126/science.1065889
3. Xie, W., Lewis, P. O., Fan, Y., Kuo, L., & Chen, M.-H. (2011). "Improving marginal likelihood estimation for Bayesian phylogenetic model selection." Systematic Biology, 60(2), 150–160. https://doi.org/10.1093/sysbio/syq085
4. Ronquist, F., Teslenko, M., van der Mark, P., Ayres, D. L., Darling, A., Höhna, S., Larget, B., Liu, L., Suchard, M. A., & Huelsenbeck, J. P. (2012). "MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space." Systematic Biology, 61(3), 539–542. https://doi.org/10.1093/sysbio/sys029