Hamiltonian Monte Carlo Algorithm Research Articles

Bayesian Inference for Mixed Gaussian GARCH-Type Model by Hamiltonian Monte Carlo Algorithm.

MCMC algorithm is widely used in parameters' estimation of GARCH-type models. However, the existing algorithms are either not easy to implement or not fast to run. In this paper, Hamiltonian Monte Carlo (HMC) algorithm, which is easy to perform and also efficient to draw samples from posterior distributions, is firstly proposed to estimate for the Gaussian mixed GARCH-type models. And then, based on the estimation of HMC algorithm, the forecasting of volatility prediction is investigated. Through the simulation experiments, the HMC algorithm is more efficient and flexible than the Griddy-Gibbs sampler, and the credibility interval of forecasting for volatility prediction is also more accurate. A real application is given to support the usefulness of the proposed HMC algorithm well.

On L 𝔮 convergence of the Hamiltonian Monte Carlo

Abstract We establish L 𝔮 {L^{\mathfrak{q}}} convergence for Hamiltonian Monte Carlo (HMC) algorithms. More specifically, under mild conditions for the associated Hamiltonian motion, we show that the outputs of the algorithms converge (strongly for 2 ≤ 𝔮 < ∞ {2\leq\mathfrak{q}<\infty} and weakly for 1 < 𝔮 < 2 {1<\mathfrak{q}<2} ) to the desired target distribution. In addition, we establish a general convergence rate for an L 𝔮 {L^{\mathfrak{q}}} convergence given a convergence rate at a specific q * {q^{*}} , and apply this result to conclude geometric convergence in the Euclidean space for HMC with uniformly strongly logarithmic concave target and auxiliary distributions. We also present the results of experiments to illustrate convergence in L 𝔮 {L^{\mathfrak{q}}} .

Bayesian inference for fitting cardiac models to experiments: estimating parameter distributions using Hamiltonian Monte Carlo and approximate Bayesian computation.

Customization of cardiac action potential models has become increasingly important with the recognition of patient-specific models and virtual patient cohorts as valuable predictive tools. Nevertheless, developing customized models by fitting parameters to data poses technical and methodological challenges: despite noise and variability associated with real-world datasets, traditional optimization methods produce a single "best-fit" set of parameter values. Bayesian estimation methods seek distributions of parameter values given the data by obtaining samples from the target distribution, but in practice widely known Bayesian algorithms like Markov chain Monte Carlo tend to be computationally inefficient and scale poorly with the dimensionality of parameter space. In this paper, we consider two computationally efficient Bayesian approaches: the Hamiltonian Monte Carlo (HMC) algorithm and the approximate Bayesian computation sequential Monte Carlo (ABC-SMC) algorithm. We find that both methods successfully identify distributions of model parameters for two cardiac action potential models using model-derived synthetic data and an experimental dataset from a zebrafish heart. Although both methods appear to converge to the same distribution family and are computationally efficient, HMC generally finds narrower marginal distributions, while ABC-SMC is less sensitive to the algorithmic settings including the prior distribution.

Hamiltonian Monte Carlo Probabilistic Joint Inversion of 2D (2.75D) Gravity and Magnetic Data

AbstractTwo‐dimensional modeling of gravity and magnetic anomalies in terms of polygonal bodies is a popular approach to infer possible configurations of geological structures in the subsurface. Alternatively to the traditional trial‐and‐error manual fit of measured data, here we illustrate a probabilistic strategy to solve the inverse problem. First we derive a set of formulae for solving a 2.75‐dimensional forward model, where the polygonal bodies have a given finite lateral extent, and then we devise a Hamiltonian Monte Carlo algorithm to jointly invert gravity and magnetic data for the geometry and properties of the polygonal bodies. This probabilistic approach fully addresses the nonlinearity of the forward model and provides uncertainty estimation. The result of the inversion is a collection of models which represent the posterior distribution, analysis of which provides estimates of sought properties and may reveal different scenarios.

Efficient Posterior Sampling for Bayesian Poisson Regression

ABSTRACT Poisson log-linear models are ubiquitous in many applications, and one of the most popular approaches for parametric count regression. In the Bayesian context, however, there are no sufficient specific computational tools for efficient sampling from the posterior distribution of parameters, and standard algorithms, such as random walk Metropolis-Hastings or Hamiltonian Monte Carlo algorithms, are typically used. Herein, we developed an efficient Metropolis-Hastings algorithm and importance sampler to simulate from the posterior distribution of the parameters of Poisson log-linear models under conditional Gaussian priors with superior performance with respect to the state-of-the-art alternatives. The key for both algorithms is the introduction of a proposal density based on a Gaussian approximation of the posterior distribution of parameters. Specifically, our result leverages the negative binomial approximation of the Poisson likelihood and the successful Pólya-gamma data augmentation scheme. Via simulation, we obtained that the time per independent sample of the proposed samplers is competitive with that obtained using the successful Hamiltonian Monte Carlo sampling, with the Metropolis-Hastings showing superior performance in all scenarios considered. Supplementary materials for this article are available online.

Bayesian models for spatially explicit interactions between neighbouring plants

Abstract Interactions between neighbouring plants drive population and community dynamics in terrestrial ecosystems. Understanding these interactions is critical for both fundamental and applied ecology. Spatial approaches to model neighbour interactions are necessary, as interaction strength depends on the distance between neighbouring plants. Recent Bayesian advancements, including the Hamiltonian Monte Carlo algorithm, offer the flexibility and speed to fit models of spatially explicit neighbour interactions. We present a guide for parameterizing these models in the Stan programming language and demonstrate how Bayesian computation can assist ecological inference on plant–plant interactions. Modelling plant neighbour interactions presents several challenges for ecological modelling. First, nonlinear models for distance decay can be prone to identifiability problems, resulting in lack of model convergence. Second, the pairwise data structure of plant–plant interaction matrices often leads to large matrices that demand high computational power. Third, hierarchical structure in plant–plant interaction data is ubiquitous, including repeated measurements within field plots, species and individuals. Hierarchical terms (e.g. ‘random effects’) can result in model convergence problems caused by correlations between coefficients. We explore modelling solutions for these challenges with examples representing spatial data on plant demographic rates: growth, survival and recruitment. We show that ragged matrices reduce computational challenges inherent to pairwise matrices, resulting in higher efficiency across data types. We also demonstrate how metrics for model convergence, including divergent transitions and effective sample size, can help diagnose problems that result from complex nonlinear structures. Finally, we explore when to use different model structures for hierarchical terms, including centred and non‐centred parameterizations. We provide reproducible examples written in Stan to enable ecologists to fit and troubleshoot a broad range of neighbourhood interaction models. Spatially explicit models are increasingly central to many ecological questions. Our work illustrates how novel Bayesian tools can provide flexibility, speed and diagnostic capacity for fitting plant neighbour models to large, complex datasets. The methods we demonstrate are applicable to any dataset that includes a response variable and locations of observations, from forest inventory plots to remotely sensed imagery. Further developments in statistical models for neighbour interactions are likely to improve our understanding of plant population and community ecology across systems and scales.

Convergence of preconditioned Hamiltonian Monte Carlo on Hilbert spaces

Abstract In this article, we consider the preconditioned Hamiltonian Monte Carlo (pHMC) algorithm defined directly on an infinite-dimensional Hilbert space. In this context, and under a condition reminiscent of strong log-concavity of the target measure, we prove convergence bounds for adjusted pHMC in the standard 1-Wasserstein distance. The arguments rely on a synchronous coupling of two copies of pHMC, which is controlled by adapting elements from Bou-Rabee, Eberle and Zimmer (2020).

Adaptive variational Hamiltonian Monte Carlo sampling with application to ocean acoustics

The Hamiltonian Monte Carlo (HMC) algorithm is a variant Markov Chain Monte Carlo (MCMC) excelling at sampling high dimensional target distributions. The HMC explores a distribution as a particle trajectory in a Hamiltonian system by augmenting distribution parameters with auxiliary momentum variables. A particle trajectory is numerically integrated from initial point to some distant low-correlation point. Projecting back to parameter space yields a proposal point with high acceptance probability leading to faster exploration of the target distribution. The Leap-frog integrator is commonly used for its volume preservation (symplecticity) and reversibility traits that maintain the detailed balance condition of MCMC. Although effective, Leap-frog has sensitivity to choice of integration step-size and number of integration steps. Poor choices cause slow exploration or bias when regions of the distribution are not explored. Adaptive step methods could remove step size issues and explore regions of a distribution where Leap-frog fails. Existing adaptive methods self-adjust step-size while maintaining symplecticity and reversibility. In this work, we apply adaptive step methods to simulated target distributions and real world ocean acoustic data and compare performance. [Work supported by Office of Naval Research.]

Fully Bayesian estimation of virtual brain parameters with self-tuning Hamiltonian Monte Carlo

Virtual brain models are data-driven patient-specific brain models integrating individual brain imaging data with neural mass modeling in a single computational framework, capable of autonomously generating brain activity and its associated brain imaging signals. Along the example of epilepsy, we develop an efficient and accurate Bayesian methodology estimating the parameters linked to the extent of the epileptogenic zone. State-of-the-art advances in Bayesian inference using Hamiltonian Monte Carlo (HMC) algorithms have remained elusive for large-scale differential-equations based models due to their slow convergence. We propose appropriate priors and a novel reparameterization to facilitate efficient exploration of the posterior distribution in terms of computational time and convergence diagnostics. The methodology is illustrated for in-silico dataset and then, applied to infer the personalized model parameters based on the empirical stereotactic electroencephalography recordings of retrospective patients. This improved methodology may pave the way to render HMC methods sufficiently easy and efficient to use, thus applicable in personalized medicine.

An efficient probabilistic workflow for estimating induced earthquake parameters in 3D heterogeneous media

Abstract. We present an efficient probabilistic workflow for the estimation of source parameters of induced seismic events in three-dimensional heterogeneous media. Our workflow exploits a linearized variant of the Hamiltonian Monte Carlo (HMC) algorithm. Compared to traditional Markov chain Monte Carlo (MCMC) algorithms, HMC is highly efficient in sampling high-dimensional model spaces. Through a linearization of the forward problem around the prior mean (i.e., the “best” initial model), this efficiency can be further improved. We show, however, that this linearization leads to a performance in which the output of an HMC chain strongly depends on the quality of the prior, in particular because not all (induced) earthquake model parameters have a linear relationship with the recordings observed at the surface. To mitigate the importance of an accurate prior, we integrate the linearized HMC scheme into a workflow that (i) allows for a weak prior through linearization around various (initial) centroid locations, (ii) is able to converge to the mode containing the model with the (global) minimum misfit by means of an iterative HMC approach, and (iii) uses variance reduction as a criterion to include the output of individual Markov chains in the estimation of the posterior probability. Using a three-dimensional heterogeneous subsurface model of the Groningen gas field, we simulate an induced earthquake to test our workflow. We then demonstrate the virtue of our workflow by estimating the event's centroid (three parameters), moment tensor (six parameters), and the earthquake's origin time. Using the synthetic case, we find that our proposed workflow is able to recover the posterior probability of these source parameters rather well, even when the prior model information is inaccurate, imprecise, or both inaccurate and imprecise.

Variational inference for nonlinear inverse problems via neural net kernels: Comparison to Bayesian neural networks, application to topology optimization

Inverse problems and, in particular, inferring unknown or latent parameters from data are ubiquitous in engineering simulations. A predominant viewpoint in identifying unknown parameters is Bayesian inference where both prior information about the parameters and the information from the observations via likelihood evaluations are incorporated into the inference process. In this paper, we adopt a similar viewpoint with a slightly different numerical procedure from standard inference approaches to provide insight about the localized behavior of unknown underlying parameters. We present a variational inference approach which mainly incorporates the observation data in a point-wise manner, i.e. we invert a limited number of observation data leveraging the gradient information of the forward map with respect to parameters, and find true individual samples of the latent parameters when the forward map is noise-free and one-to-one. For statistical calculations (as the ultimate goal in simulations), a large number of samples are generated from a trained neural network which serves as a transport map from the prior to posterior latent parameters. Our neural network machinery, developed as part of the inference framework and referred to as Neural Net Kernels (NNK), is based on hierarchical (deep) kernels which provide greater flexibility for training compared to standard neural networks. We showcase the effectiveness of our inference procedure in identifying bimodal and irregular distributions compared to a number of approaches including a maximum a posteriori probability (MAP)-based approach, Markov Chain Monte Carlo sampling with both Metropolis–Hastings and Hamiltonian Monte Carlo algorithms, and a Bayesian neural network approach, namely the widely-known Bayes by Backprop algorithm via a pedagogical example. We further apply our inference procedure to two and three dimensional topology optimization problems where we identify the latent parameters in the random field elastic modulus, modeled as a Karhunen–Loéve expansion, considering the high dimensional design iterate-displacement pair as training data. As future research, we discuss the application of this inference approach in identifying constitutive models in nonlinear elasticity and development of fast linear algebraic solvers for large scale similar, i.e. finite element-based inverse problems calculations.

Seismic inversion for density using a transdimensional approach

Commercial and low-saturation gas (also called paleoresidual gas [PRG]) show similar strong amplitude signatures on P-wave seismic data. This poses an exploration risk in gas reservoir regions. However, density correlates inversely with gas saturation and can differentiate a zone of full gas saturation from PRG. This can improve the chances of success in terms of predrill prediction of gas saturation. Amplitude-variation-with-offset (AVO) inversion using prestack seismic data is the most commonly used technique that can estimate elastic parameters such as P-wave velocity, S-wave velocity, and density. Out of these three parameters, extracting density from seismic data is the most challenging due to its weak sensitivity to seismic reflection amplitude and the lack of good quality seismic data at far offsets. However, with recent improvements in seismic data acquisition and processing technology, which produces reliable AVO gathers, density estimates have improved. This requires that strong density sensitivity to AVO exists. Note that multiple density models may fit the data equally well. Therefore, quantifying uncertainty is crucial for interpretation and risk assessment. We apply a recently developed stochastic approach based on the Bayesian framework to solve the problem in a transdimensional framework, where the number of model parameters is treated as a variable and estimated along with the elastic properties. We use the reversible jump Hamiltonian Monte Carlo (RJHMC) algorithm to sample models from a variable dimensional model space and obtain a globally optimum model and uncertainty estimates. We use a synthetic and good quality real data set from Columbus Basin in Trinidad, which has a proven gas reservoir, to demonstrate the algorithm. The RJHMC results calibrate well with the logs and show the areal extents of the density anomalies within the 3D volume.

Hamiltonian Monte Carlo with strict convergence criteria reduces run-to-run variability in forensic DNA mixture deconvolution

Motivation: Analysing mixed DNA profiles is a common task in forensic genetics. Due to the complexity of the data, such analysis is often performed using Markov Chain Monte Carlo (MCMC)-based genotyping algorithms. These trade off precision against execution time. When default settings (including default chain lengths) are used, as large as a 10-fold changes in inferred log-likelihood ratios (LR) are observed when the software is run twice on the same case. So far, this uncertainty has been attributed to the stochasticity of MCMC algorithms. Since LRs translate directly to strength of the evidence in a criminal trial, forensic laboratories desire LR with small run-to-run variability. Results:We present the use of a Hamiltonian Monte Carlo (HMC) algorithm that reduces run-to-run variability in forensic DNA mixture deconvolution by around an order of magnitude without increased runtime. We achieve this by enforcing strict convergence criteria. We show that the choice of convergence metric strongly influences precision. We validate our method by reproducing previously published results for benchmark DNA mixtures (MIX05, MIX13, and ProvedIt). We also present a complete software implementation of our algorithm that is able to leverage GPU acceleration for the inference process. In the benchmark mixtures, on consumer-grade hardware, the runtime is less than 7 min for 3 contributors, less than 35 min for 4 contributors, and less than an hour for 5 contributors with one known contributor.

A Bayesian multivariate model using Hamiltonian Monte Carlo inference to estimate total organic carbon content in shale

The prediction of total organic carbon (TOC) content using geophysical logs is one of the key steps in shale reservoir characterization. Various empirical relations have previously been used for the estimation of TOC content from well-logs; however, uncertainty quantification in the model estimation is often ignored while performing TOC estimation in a deterministic framework. We introduce the problem of TOC estimation in a Bayesian setting with the goal of enhancing the TOC content prediction together with the quantification of the uncertainty in the model prediction. To signify the uncertainty, we draw random samples of model parameters from the posterior distribution by realizing multidimensional stochastic processes within the Hamiltonian Monte Carlo algorithm. The posterior model for the variables that influence TOC estimation is conditioned on the available well-log observations and is further defined by a priori and likelihood distributions. We demonstrate examples of applications of this approach to estimate the TOC content on two real field data sets from the well-known Devonian Duvernay shale of Western Canada and the Silurian shale of the Ahnet Basin. The accuracy in the estimation is arbitrated by comparing the prediction results with those obtained using the two most widely used empirical models. Corroborating the results by the laboratory-measured TOC contents demonstrate that the Bayesian approach offers a more reliable and better confidence in predictions when compared with the empirical models, as it provides additional information on the prediction uncertainty. Finally, the implications of the present approach are derived in terms of depositional environments to characterize the high TOC content zone in the studied organic shale formations.

Meta-analysis of dichotomous and ordinal tests with an imperfect gold standard.

Standard methods for the meta‐analysis of medical tests, without assuming a gold standard, are limited to dichotomous data. Multivariate probit models are used to analyse correlated dichotomous data, and can be extended to model ordinal data. Within the context of an imperfect gold standard, they have previously been used for the analysis of dichotomous and ordinal test data from a single study, and for the meta‐analysis of dichotomous tests. However, they have not previously been used for the meta‐analysis of ordinal tests. In this article, we developed a Bayesian multivariate probit latent class model for the simultaneous meta‐analysis of ordinal and dichotomous tests without assuming a gold standard, which also allows one to obtain summary estimates of joint test accuracy. We fitted the models using the software Stan, which uses a state‐of‐the‐art Hamiltonian Monte Carlo algorithm, and we applied the models to a dataset in which studies evaluated the accuracy of tests, and test combinations, for deep vein thrombosis. We demonstrate the issues with dichotomising ordinal test accuracy data in the presence of an imperfect gold standard, before applying and comparing several variations of our proposed model which do not require the data to be dichotomised. The models proposed will allow researchers to more appropriately meta‐analyse ordinal and dichotomous tests without a gold standard, potentially leading to less biased estimates of test accuracy. This may lead to a better understanding of which tests, and test combinations, should be used for any given medical condition.

Crop specific inversion of PROSAIL to retrieve green area index (GAI) from several decametric satellites using a Bayesian framework

The objective of this study is to evaluate the performances of a semi-empirical approach based on the Bayesian theory to retrieve Green Area Index (GAI) from multiple decametric satellites. It is designed to overcome some limitations in existing Radiative Transfer Model (RTM) inversion methods, including the high dimensionality of the inverse problem, the convergence problem due to possible equifinality, and the dependence of some RTM variables on the crop-specific architecture. The PROSAIL model is first inverted in a calibration step using the Hamiltonian Monte Carlo (HMC) algorithm over a global dataset of ground GAI measurements (for maize, wheat, and rice) and the corresponding reflectance observations from Landsat-8, Sentinel-2, and Quickbird to derive crop-specific distributions of PROSAIL input variables. These distributions were then used as prior information to predict GAI over an independent set of reflectance observations.Results show that the full Bayesian approach provides close estimates of GAI to ground truth, with respective Root Mean Square Error (RMSE) of 1.01, 1.33, and 0.97 for maize, wheat, and rice (R2=0.67, 0.76 and 0.63, respectively). The performances are better than those approaches generally reported using radiative transfer models that are non-crop-specific, like the SNAP algorithm for Sentinel-2, but are slightly behind the purely empirical models based on machine learning. However, the proposed approach provides an explicit insight of the joint distribution of PROSAIL variables that are valid for any satellite platform. This constitutes a major advantage against purely empirical models, as it enables to fully exploit large observational datasets from multiple sensors and generalize to other platforms.

Hamiltonian Monte Carlo reconstruction from peculiar velocities

ABSTRACTThe problem of the reconstruction of the large-scale density and velocity fields from peculiar velocity surveys is addressed here within a Bayesian framework by means of Hamiltonian Monte Carlo (HMC) sampling. The HAmiltonian Monte carlo reconstruction of the Local EnvironmenT (hamlet) algorithm is designed to reconstruct the linear large-scale density and velocity fields in conjunction with the undoing of lognormal bias in the derived distances and velocities of peculiar velocity surveys, such as the Cosmicflows (CF) data. The hamlet code has been tested against CF mock catalogues consisting of up to 3 × 104 data points with mock errors akin to those of the Cosmicflows-3 (CF3) data, within the framework of the Lambda cold dark matter standard model of cosmology. The hamlet code outperforms previous applications of Gibbs sampling Markov chain Monte Carlo reconstruction from the CF3 data by two to four orders of magnitude in CPU time. The gain in performance is due to the inherent higher efficiency of the HMC algorithm and due to parallel computing on GPUs rather than CPUs. This gain will enable an increase in the reconstruction of the large-scale structure from the upcoming CF4 data and the setting of constrained initial conditions for cosmological high-resolution simulations.

Couplings for Andersen dynamics

Andersen dynamics is a standard method for molecular simulations, and a precursor of the Hamiltonian Monte Carlo algorithm used in MCMC inference. The stochastic process corresponding to Andersen dynamics is a PDMP (piecewise deterministic Markov process) that iterates between Hamiltonian flows and velocity randomizations of randomly selected particles. Both from the viewpoint of molecular dynamics and MCMC inference, a basic question is to understand the convergence to equilibrium of this PDMP particularly in high dimension. Here we present couplings to obtain sharp convergence bounds in the Wasserstein sense that do not require global convexity of the underlying potential energy.

Bayesian analysis of multiple break-points threshold ARMA model with exogenous inputs

In this paper, we introduce a Bayesian statistical inference approach for multiple break-points threshold autoregressive moving average model with exogenous inputs (MB-TARMAX) which change in state space and time domain. Based on the appropriate prior information of parameters, we give the full conditional posterior distribution of parameters including the thresholds and break-points. In order to obtain the estimates of parameters, we employ the Markov chain Monte Carlo (MCMC) method via Gibbs sampler with Metropolis-Hastings algorithm. Compared with Metropolis-Hastings algorithm, we apply Hamiltonian Monte Carlo algorithm to avoid the slow space exploration from simple random walk and improve the sampling efficiency. As applications, we demonstrate the effectiveness of our method from simulation experiments and a real example.

Manifold Markov Chain Monte Carlo Methods for Bayesian Inference in Diffusion Models

AbstractBayesian inference for nonlinear diffusions, observed at discrete times, is a challenging task that has prompted the development of a number of algorithms, mainly within the computational statistics community. We propose a new direction, and accompanying methodology—borrowing ideas from statistical physics and computational chemistry—for inferring the posterior distribution of latent diffusion paths and model parameters, given observations of the process. Joint configurations of the underlying process noise and of parameters, mapping onto diffusion paths consistent with observations, form an implicitly defined manifold. Then, by making use of a constrained Hamiltonian Monte Carlo algorithm on the embedded manifold, we are able to perform computationally efficient inference for a class of discretely observed diffusion models. Critically, in contrast with other approaches proposed in the literature, our methodology is highly automated, requiring minimal user intervention and applying alike in a range of settings, including: elliptic or hypo-elliptic systems; observations with or without noise; linear or non-linear observation operators. Exploiting Markovianity, we propose a variant of the method with complexity that scales linearly in the resolution of path discretisation and the number of observation times. Python code reproducing the results is available at http://doi.org/10.5281/zenodo.5796148.