PaperlayerView exact revision on arXiv

Source-linked AI summary

Markov Chain Monte Carlo Methods for Bayesian Data Analysis in Astronomy

Sanjib Sharma

arXiv:1706.01629v1astro-ph.IMphysics.comp-phstat.CO

TL;DR

Astronomical Bayesian analyses require a framework for expressing inference problems and numerical methods for evaluating generally intractable distributions. This review explains Bayesian foundations, surveys MCMC algorithms, presents applications and advanced models, and distributes software implementing selected methods and examples. It emphasizes efficient sampling and extensions for complex problems while noting that thinning discards information and that full distributions or samples can better communicate results.

  • Problem

    Astronomical Bayesian analyses require suitable ways to formulate problems and efficient MCMC methods, whose performance can suffer from tuning, dimensionality, correlations, degeneracies, and multimodality.

  • Method

    The review explains Bayesian theory, surveys Monte Carlo methods, discusses applications and advanced models, and provides software implementing selected algorithms and examples.

  • Results

    The review presents a broad account of Bayesian-MCMC data analysis in astronomy, including methods for sampling, hierarchical modelling, uncertainty handling, and stellar-parameter estimation.

  • Takeaways & Limitations

    The review offers introductory guidance, recent developments, diverse applications, and advanced Bayesian ideas for addressing complex astronomical data-analysis problems.

  • Takeaways & Limitations

    Thinning is useful only when autocorrelation is known beforehand and large; discarding samples cannot improve estimates over the full chain.

Abstract

from arXiv · show

Markov Chain Monte Carlo based Bayesian data analysis has now become the method of choice for analyzing and interpreting data in almost all disciplines of science. In astronomy, over the last decade, we have also seen a steady increase in the number of papers that employ Monte Carlo based Bayesian analysis. New, efficient Monte Carlo based methods are continuously being developed and explored. In this review, we first explain the basics of Bayesian theory and discuss how to set up data analysis problems within this framework. Next, we provide an overview of various Monte Carlo based methods for performing Bayesian data analysis. Finally, we discuss advanced ideas that enable us to tackle complex problems and thus hold great promise for the future. We also distribute downloadable computer software (available at https://github.com/sanjibs/bmcmc/ ) that implements some of the algorithms and examples discussed here.

1. Introduction

The review explains how Bayesian theory and MCMC combine to analyze astronomical data, emphasizing numerical sampling, efficient algorithms, applications, and advanced models. It also targets both newcomers and experienced readers while providing software implementing selected algorithms and examples.

  • Bayesian theory provides a unified framework for inferring causes or their probabilities from observed outcomes and prior information.
  • MCMC numerically samples parameter distributions such as p(θ|D) when analytical solutions are unavailable.The parameters may represent a set rather than a single quantity.
  • MCMC efficiency remains challenging because proposal tuning worsens with dimensionality, correlations, degeneracies, and multimodal distributions.The review therefore focuses on efficient algorithms, including tempering methods for chains that become stuck near local density maxima.
  • Advanced topics include non-parametric and hierarchical Bayesian models, which the review describes as powerful but not yet mainstream in astronomy.
  • The review explains Bayesian basics, surveys recent developments, and discusses emerging ideas for complex problems.It is aimed at readers ranging from newcomers to those already familiar with Bayesian data analysis.
  • The authors distribute downloadable software implementing some discussed algorithms and examples.The software is available through the bmcmc GitHub repository or the bmcmc package command.

2. Bayesian Data Analysis

This section introduces Bayesian theory for modeling data, incorporating priors and measurement errors, and illustrates inference, outlier modeling, and model comparison. It emphasizes that predictive criteria and Bayesian methods serve different goals.

  • Bayesian foundations: Bayesian analysis updates prior knowledge with data through likelihood, posterior, and evidence quantities.The posterior expresses updated belief, while the evidence normalizes the posterior and supports model comparison.
  • Modeling data: A generative model describes how data arise from parameters, while measurement uncertainty is incorporated by integrating over unknown true values.The framework allows Gaussian measurement errors and explicitly integrates over latent true values.
  • Modeling data: Outlier modeling supplements the background model with an outlier probability and separate parameters.The straight-line example generated 50 points with m = 2.0 and c = 10.0, including 20% outliers.
  • Modeling data: MCMC posterior sampling produced best-fit parameters that resembled the true parameters in the synthetic straight-line example.The example is implemented in the distributed software.
  • Model comparison: WAIC uses Bayesian predictive density, handles singular models better than AIC, and is asymptotically equivalent to Bayesian LOOCV.DIC and WAIC also estimate effective degrees of freedom from likelihood information and posterior samples.
  • Model comparison: When priors are well justified, Bayes-factor methods suit model selection, whereas WAIC is preferable when the goal is predictive accuracy for future data.The section presents Bayesian and predictive methods as having distinct strengths and weaknesses.

3. Monte Carlo methods for Bayesian computations

This section explains MCMC as numerical sampling for Bayesian distributions that generally lack closed-form solutions. It develops Markov-chain foundations and compares Metropolis-Hastings proposals and acceptance behavior.

  • Foundations: MCMC numerically samples parameter distributions when Bayesian posterior calculations lack analytical solutions.Monte Carlo estimates use averages from the output of a Markov chain.
  • Markov chains: A Markov chain makes the future conditionally independent of the past given its present state, and time-homogeneous chains use transition probabilities independent of iteration.The state space may be continuous or discrete.
  • Markov chains: Irreducibility, positive recurrence, and aperiodicity support convergence toward a unique limiting stationary distribution.Under these conditions, the chain’s distribution approaches equilibrium from any initial distribution.
  • Markov chains: Detailed balance is commonly imposed because reversibility guarantees a stationary distribution, although reversibility is not necessary.Most MCMC algorithms are designed to satisfy detailed balance.
  • Metropolis-Hastings: Metropolis-Hastings proposes a new state from q(y|x) and accepts it using an acceptance ratio that preserves the target distribution.The transition kernel is K(x, y) = q(y|x)α(x, y).
  • Metropolis-Hastings: The Metropolis-Hastings acceptance ratio yields the fastest mixing rate among the considered reversible alternatives.The transition kernel satisfies detailed balance with respect to f.
  • Metropolis-Hastings: Proposal choices include symmetric Metropolis, random-walk, independence, and Langevin algorithms.Langevin proposals use gradient information to allow faster chain mixing.
  • Metropolis-Hastings: For random-walk Metropolis-Hastings, decreasing σ raises acceptance but lowers independence, while increasing σ has the opposite trade-off.The mean acceptance ratio is below unity except for a uniform target density.

3.3. Gibbs sampling

Gibbs sampling constructs a Markov chain by updating coordinates from conditional distributions. Its stationary-distribution property and block-update structure make it useful for complicated multivariate distributions.

  • Gibbs sampling: Gibbs sampling splits the transition kernel into coordinate-wise updates using each coordinate’s conditional density given the others.The method targets f(x) on a d-dimensional state space.
  • Stationarity and reversibility: The systematic-scan Gibbs transition kernel has f as its stationary distribution but is not reversible.Reversible variants can be obtained by selecting a component randomly at each iteration.
  • Stationarity and reversibility: Random-scan Gibbs sampling can be viewed as a special case of Metropolis-Hastings.Its acceptance ratio is expressed through the corresponding proposal probabilities.

3.4. Metropolis within Gibbs

Metropolis within Gibbs replaces difficult conditional draws with Metropolis-Hastings steps. The approach can also update subsets of variables jointly through block updates.

  • Metropolis within Gibbs: Metropolis within Gibbs uses an MH step when sampling directly from a Gibbs conditional distribution is difficult.It is presented as an alternative implementation of coordinate-wise Gibbs updates.
  • Metropolis within Gibbs: The method is implemented in the software distributed with the review.The algorithm is identified as Algorithm 3.
  • Block updates: Block updates allow subsets of variables to be updated together instead of advancing every variable separately.Breaking a complicated distribution into smaller, easier samplings is described as Gibbs sampling’s main strength.

3.5. Adaptive Metropolis

Adaptive Metropolis methods update the proposal distribution during sampling to reduce user tuning while preserving convergence to the target distribution. The review presents a general adaptive algorithm and illustrates that its chain can resemble an ideally tuned random-walk chain.

  • Motivation: The random-walk Metropolis algorithm becomes difficult to tune in high-dimensional or correlated spaces because proposal efficiency depends on the supplied covariance.Small proposals yield high acceptance but slow mixing; large proposals yield low acceptance and also slow mixing.
  • Adaptive principle: Adaptive schemes estimate and update proposal parameters using past samples while ensuring the proposal converges toward an optimal choice.The adaptation addresses the need to estimate the target covariance before the actual run.
  • Adaptive principle: Robbins–Monro recursion provides the stochastic updating framework, with step sizes γn = γ/n^β for 0 < β < 1.The review states that the iterates converge to the true solution under the stated conditions.
  • Algorithm: The adaptive symmetric random-walk algorithm initializes a point, mean, covariance, target acceptance, and function, then updates scale, mean, and covariance after each proposal.Each proposal is drawn from N(xi, λiΣi), accepted using the current ratio, and followed by parameter updates.
  • Results: β = 0.6 is satisfactory for most test cases, and the resulting adaptive chain closely resembles the ideal chain in Figure 5b.Figure 5d shows the adaptive chain produced by Algorithm 4.

3.6. Affine invariant sampling

Affine-invariant ensemble samplers use interacting walkers to avoid manually selecting a suitable proposal distribution. The review also explains how to assess convergence and effective sample size, while warning that thinning discards information.

  • Affine invariant sampling: Ensemble samplers run multiple interacting chains whose collective information adapts proposal densities, reducing the need for problem-specific proposal tuning.The affine invariant sampler is presented as a general-purpose ensemble method.
  • Affine invariant sampling: Partial resampling updates one walker conditional on the others while satisfying detailed balance, so each walker samples from the target distribution π(x).The ensemble targets the product distribution Π(X) = π(x1)π(x2)...π(xL).
  • Affine invariant sampling: The stretch move proposes y = xi + (r −1)(xi −xj) = xj + r(xi −xj), using another walker and a scalar stretch factor.The proposal is constructed along the line joining two walkers, requiring an appropriate Jacobian in the acceptance calculation.
  • Limitations: Affine-invariant sampling elegantly addresses proposal selection, but its warm-up computational cost scales linearly with the ensemble size.This is identified as a drawback of the Goodman–Weare algorithm.
  • Convergence diagnostics: Effective sample size is N/(2τint,x), accounting for the variance inflation caused by autocorrelation in MCMC samples.The review suggests using ESS to judge whether enough effectively independent samples have been obtained.
  • Convergence diagnostics: Thinning provides no inferential advantage beyond storage savings because discarding correlated samples cannot improve estimates from the original chain.Choosing the thinning interval also requires examining the full-chain autocorrelation.

3.8. Parallel Tempering

Parallel tempering addresses multimodal target distributions by running chains at different temperatures and exchanging states. Temperature broadening helps explore parameter space, but the ladder must balance swap acceptance against computational cost.

  • Method: Multimodal distributions challenge MCMC algorithms, motivating parallel tempering with multiple chains simulated at different temperatures.The target chain has temperature Tn = 1, while higher-temperature chains explore broadened distributions.
  • Method: Higher temperatures broaden the target distribution, enabling wider exploration of parameter space before information is exchanged through state swaps.The chains form a ladder T1 > T2 > ... > Tn, with Tn set to 1.
  • Limitations: The temperature ladder requires careful spacing: distant neighbors reduce swap acceptance and slow mixing, whereas close neighbors increase the number of chains and computational cost.Trial runs can be used to construct a suitable ladder.
  • Extensions: Evolutionary MCMC extends the parallel-chain framework by exchanging information through mutation and crossover ideas while maintaining detailed balance.The construction operates on the product space defined by the parallel chains.

3.9. Monte Carlo Metropolis Hastings

Monte Carlo Metropolis–Hastings approximates Bayesian sampling when likelihoods or posterior ratios cannot be evaluated analytically. It introduces auxiliary samples to estimate the acceptance ratio, with fresh auxiliary draws improving robustness against stochastic trapping.

  • Algorithm: Monte Carlo Metropolis–Hastings generates proposals and auxiliary samples, estimates the MH ratio, and accepts proposals using the estimated ratio.The algorithm outputs points sampled approximately from p(θ|x).
  • Motivation: The method targets cases where p(x|θ) or the posterior cannot be easily evaluated in analytically tractable form.Auxiliary Monte Carlo calculations provide an approximate evaluation of p(θ|x).
  • Robustness: Refreshing the auxiliary sample each iteration makes the algorithm more robust than alternatives that reuse a previous stochastic posterior estimate.Reusing an estimate can cause the chain to become stuck at a stochastic maximum of the estimated likelihood.
  • Robustness: Smaller auxiliary samples or less accurate Monte Carlo ratio estimates worsen the tendency of the chain to become stuck at stochastic maxima.Using the same auxiliary sample for both terms in the ratio reduces noise in the estimated ratio.
  • Unknown normalization constant: For unknown normalization constants, the MH ratio separates into a normalization-constant ratio and a likelihood-shape ratio, with the former estimated using auxiliary samples.When exact sampling is infeasible, Monte Carlo or importance sampling can estimate the unknown ratio.

3.10. Hamiltonian Monte Carlo

Hamiltonian Monte Carlo addresses slow exploration in high-dimensional posteriors by using Hamiltonian dynamics to make distant moves. Population Monte Carlo offers adaptive importance sampling, while importance-sampling estimates are unbiased under the stated conditions.

  • Hamiltonian Monte Carlo: In high dimensions, posterior typical sets form large-volume shells that are difficult to traverse with large steps, slowing traditional MCMC exploration.The difficulty arises because the chain takes a long time to explore the shell-shaped typical set.
  • Hamiltonian Monte Carlo: HMC introduces momentum variables and uses Hamiltonian dynamics before a Metropolis-Hastings accept/reject step.The dynamics allow proposals to travel far from the current location, enabling faster parameter-space exploration.
  • Hamiltonian Monte Carlo: HMC requires target-density gradients and tuning of both the step size ϵ and the number of evolution steps.These requirements have hindered widespread use, although progress has addressed both obstacles.
  • Population Monte Carlo: Population Monte Carlo iteratively adapts an importance function using previous samples and produces approximately target-distributed samples with importance weights.The weighted samples can estimate quantities integrated over the target distribution.
  • Population Monte Carlo: Importance-sampling estimates have lower variance when the importance function is closer to the target, but selecting a good function is difficult in practice.The method can use a mixture of normal or t-distributions as a simple importance-function choice.
  • Population Monte Carlo: Under the stated construction, importance sampling gives an unbiased estimate of the expectation of any function h(x).The equality is valid for any proposal distribution g.

3.12. Nested Sampling

Nested sampling transforms evidence computation into one-dimensional integration over prior mass ordered by likelihood. Its iterative replacement scheme concentrates samples toward higher likelihood and supports cosmological model selection and parameter estimation.

  • Motivation: Nested sampling is designed to ease the computational challenge of computing Bayesian evidence.The method targets integrals associated with the evidence rather than relying on direct high-dimensional integration.
  • Evidence transformation: The method labels iso-likelihood contours with enclosed prior mass X and rewrites the evidence as an integral involving the inverse likelihood function L(X).For normalized priors, X ranges from 0 at the highest likelihood to 1 at the lowest likelihood.
  • Iterative algorithm: Nested sampling starts with N prior draws, repeatedly removes the lowest-likelihood point, and replaces it with a prior draw satisfying L > L_i.The expected enclosed prior mass after iteration i is X_i = exp(−i/N).
  • Iterative algorithm: The evidence integral should be performed in equal steps of ln(X), because small-prior-mass regions can contribute substantially.This motivates the iterative procedure used to generate the likelihood sequence.
  • Applications: Nested sampling is widely used for cosmological model selection and parameter estimation, with packages including CosmoNest, MultiNest, and DNEST.The review identifies all three as publicly available nested-sampling packages.

4. Bayesian hierarchical modelling (BHM)

Bayesian hierarchical modelling connects group-specific and shared information while accommodating hyperparameters, hidden variables, and missing data. The review also relates EM, data augmentation, Gibbs sampling, and stochastic methods to practical hierarchical inference.

  • Bayesian hierarchical modelling: Hierarchical Bayesian models represent dependencies among parameters, hyperparameters, and hidden variables through joint conditional structures.They provide a framework for marginalizing unknown or missing quantities needed to model the data.
  • Bayesian hierarchical modelling: When groups are small, group-specific means are uncertain, so relying partly on the global mean can be preferable to using each group mean alone.For example, with n_j = 2, the group-mean uncertainty is large.
  • Bayesian hierarchical modelling: Hierarchical modelling provides a middle ground between global and group-specific means by using information from all groups.This pooling makes the estimates systematically closer to the global mean than standard group means.
  • Bayesian hierarchical modelling: In the synthetic group-mean example, hierarchical estimates had smaller error bars and shifted more toward the global mean when group error bars were larger.The model used information available from other groups in estimating each group mean.
  • Inference algorithms: Gibbs sampling is presented as an accessible way to analyze hierarchical models, with data augmentation and Metropolis-within-Gibbs extending inference to missing data and unknown parameters.The sequence of sampled parameters and latent variables forms a Markov chain targeting their joint posterior.
  • Inference algorithms: EM alternates expectation and maximization steps, increasing marginal likelihood at each iteration but not guaranteeing a global maximum.The algorithm can converge to a saddle point or local maximum; stochastic variants instead generate a distribution and are less likely to get stuck at a local maximum.
  • Inference algorithms: In a 100,000-iteration comparison, data augmentation and explicit integration produced identical parameter-estimate pdfs, while data augmentation had slightly higher autocorrelation time.The higher autocorrelation resulted from sampling an extra parameter for each data point.

5. Case studies in astronomy

The review presents diverse astronomical case studies showing how Bayesian models and MCMC sampling are used for inference, population modelling, stellar spectroscopy, and incorporation of prior information. These examples emphasize both the breadth of applications and the practical importance of choosing suitable models and sampling schemes.

  • Case-study scope: Astronomical applications include radial-velocity inference, population modelling of binary systems, stellar-parameter estimation from spectra, and Milky Way dynamical modelling.The cases span exoplanets and binary systems, hierarchical population models, data-driven spectroscopy, and analyses incorporating external constraints.
  • Exoplanets and binary systems: Radial-velocity data can constrain binary-system parameters such as mass ratio, orbital period, eccentricity, and systemic velocity using Bayesian posterior sampling.The review models seven orbital parameters and samples their posterior distribution with MCMC.
  • Exoplanets and binary systems: Different eccentricities produce different radial-velocity curve shapes, demonstrating that the observations are sensitive to orbital eccentricity.The example compares systems sharing other parameters while varying e, and the MCMC analysis recovers posterior distributions for selected parameters.
  • Population modelling: Hierarchical Bayesian models avoid population estimates based only on maximum-a-posteriori orbital parameters, which ignore uncertainty in individual systems.The hierarchical formulation estimates system-level parameters alongside population hyperparameters, with Metropolis-within-Gibbs sampling available for the resulting posterior.
  • Stellar spectroscopy: For stellar spectroscopy, a probabilistic generative model maps stellar labels such as temperature, gravity, and abundances to observed spectra, enabling inference of labels for new stars.The Cannon uses fitted model parameters to estimate stellar parameters from a new spectrum and is illustrated with best-fit spectra for four stars.
  • Prior information: External measurements, including solar angular velocity and local surface mass density, can be incorporated as prior information or additional constraints in Bayesian analyses.The review gives ω = 30.24 ± 0.12 km s−1kpc−1 and Σ1.1,⊙ = 72 ± 6 as examples of such constraints.

6. Concluding remarks

The concluding remarks describe MCMC and Bayesian methods as increasingly capable tools for complex astronomical inference, while identifying persistent challenges in computation, multimodality, communication, and high-dimensional hierarchical exploration. The authors also point to future needs including better samplers, visualization, and non-parametric Bayesian methods.

  • Current capabilities: Robust MCMC algorithms now sample multidimensional and complex probability distributions, with adaptive Metropolis, affine-invariant samplers, and parallel tempering addressing important practical difficulties.These methods reduce proposal-tuning problems or support sampling of distributions with multiple challenging features.
  • Model comparison: Bayesian evidence remains difficult to compute efficiently, although alternative model-comparison criteria can use an already computed MCMC chain.The conclusion distinguishes the availability of comparison criteria from the continuing challenge of efficient evidence calculation.
  • Hierarchical models: Bayesian hierarchical models address missing data, marginalization, and observational uncertainty, and Metropolis-within-Gibbs is presented as suitable for sampling their posteriors.The review also provides software for this class of hierarchical-model analyses.
  • Open computational challenges: Multimodal distributions remain problematic for most MCMC algorithms, while parallel tempering can help at additional computational cost and with careful ladder selection.Very high-dimensional, complex distributions are also difficult to explore efficiently, motivating derivative-based methods such as Hamiltonian Monte Carlo.
  • Communicating results: Reporting full posterior distributions or thinned samples preserves variable correlations and supports feeding results from one MCMC analysis into another.Analytical approximations and Gaussian mixture models are offered as more storage-efficient alternatives, alongside better visualization tools for complex outputs.
  • Future directions: Future astronomical datasets will be larger and more precise, increasing demand for computationally efficient methods for higher-dimensional and multi-level models.The review identifies non-parametric Bayesian methods, including Gaussian processes and Dirichlet process mixtures, as topics it did not address.

DISCLOSURE STATEMENT

The authors report no affiliations, memberships, funding, or financial holdings known to affect the review’s objectivity.

  • No relevant affiliations, memberships, funding, or financial holdings were reported.
Loading 1706.01629v1…