“Parallel-Tempering”-Assisted Hybrid Monte Carlo Algorithm for Bayesian Inference in Dynamical Systems

Bayesian model updating provides a rigorous framework to account for uncertainty induced by lack of knowledge about engineering systems in their respective mathematical models through updates of the joint probability density function (PDF), the so-called posterior PDF, of the unknown model parameters. The Markov chain Monte Carlo (MCMC) methods are currently the most popular approaches for generating samples from the posterior PDF. However, these methods often found wanting when sampling from difficult distributions (e.g., high-dimensional PDFs, PDFs with flat manifolds, multimodal PDFs, and very peaked PDFs). This paper introduces a new multi-level sampling approach for Bayesian model updating, called Sequential Gauss-Newton algorithm, which is inspired by the Transitional Markov chain Monte Carlo (TMCMC) algorithm. The Sequential Gauss-Newton algorithm improves two aspects of TMCMC to make an efficient and effective MCMC algorithm for drawing samples from difficult posterior PDFs. First, the statistical efficiency of the algorithm is enhanced by use of the systematic resampling scheme. Second, a new MCMC algorithm, called Gauss-Newton MCMC algorithm, is proposed which is essentially an M-H algorithm with a Gaussian proposal PDF tailored to the posterior PDF using the gradient and Hessian information of the negative log posterior. The effectiveness of the proposed algorithm for solving the Bayesian model updating problem is illustrated using three examples with irregularly shaped posterior PDFs.

We apply the hybrid Monte Carlo (HMC) algorithm to the financial time sires analysis of the stochastic volatility (SV) model for the first time. The HMC algorithm is used for the Markov chain Monte Carlo (MCMC) update of volatility variables of the SV model in the Bayesian inference. We compute parameters of the SV model from the artificial financial data and compare the results from the HMC algorithm with those from the Metropolis algorithm. We find that the HMC decorrelates the volatility variables faster than the Metropolis algorithm. We also make an empirical analysis based on the Yen/Dollar exchange rates.

<strong class="journal-contentHeaderColor">Abstract.</strong> 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. We find that our 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.

The Hamiltonian Monte Carlo (HMC) algorithm is a Markov Chain Monte Carlo (MCMC) technique, which combines the advantages of Hamiltonian dynamics methods and Metropolis Monte Carlo approach, to sample from complex distributions. The HMC algorithm incorporates gradient information in the dynamic trajectories and thus suppresses the random walk nature in traditional Markov chain simulation methods. This ensures rapid mixing, faster convergence, and improved efficiency of the Markov chain. The leapfrog method is generally used in discrete simulation of the dynamic transitions. In this paper, we refer to this as the leapfrog–HMC. The primary goal of this paper is to present the HMC algorithm as a tool for rapid sampling of high dimensional and complex distributions, and demonstrate its advantages over the classical Metropolis Monte Carlo technique. We demonstrate that the use of an adaptive–step discretization scheme in simulating the dynamic transitions results in an algorithm which significantly outperforms the leapfrog–HMC algorithm. Relevance to reservoir parameter estimation and uncertainty quantification will be discussed.

Sampling across tree space is one of the major challenges in Bayesian phylogenetic inference using Markov chain Monte Carlo (MCMC) algorithms. Standard MCMC tree moves consider small random perturbations of the topology, and select from candidate trees at random or based on the distance between the old and new topologies. MCMC algorithms using such moves tend to get trapped in tree space, making them slow in finding the globally most probable trees (known as “convergence”) and in estimating the correct proportions of the different types of them (known as “mixing”). Here, we introduce a new class of moves, which propose trees based on their parsimony scores. The proposal distribution derived from the parsimony scores is a quickly computable albeit rough approximation of the conditional posterior distribution over candidate trees. We demonstrate with simulations that parsimony-guided moves correctly sample the uniform distribution of topologies from the prior. We then evaluate their performance against standard moves using six challenging empirical data sets, for which we were able to obtain accurate reference estimates of the posterior using long MCMC runs, a mix of topology proposals, and Metropolis coupling. On these data sets, ranging in size from 357 to 934 taxa and from 1740 to 5681 sites, we find that single chains using parsimony-guided moves usually converge an order of magnitude faster than chains using standard moves. They also exhibit better mixing, that is, they cover the most probable trees more quickly. Our results show that tree moves based on quick and dirty estimates of the posterior probability can significantly outperform standard moves. Future research will have to show to what extent the performance of such moves can be improved further by finding better ways of approximating the posterior probability, taking the trade-off between accuracy and speed into account. [Bayesian phylogenetic inference; MCMC; parsimony; tree proposal.]

One of the most common types of models that helps us to understand neuron behavior is based on the Hodgkin–Huxley ion channel formulation (HH model). A major challenge with inferring parameters in HH models is non-uniqueness: many different sets of ion channel parameter values produce similar outputs for the same input stimulus. Such phenomena result in an objective function that exhibits multiple modes (i.e., multiple local minima). This non-uniqueness of local optimality poses challenges for parameter estimation with many algorithmic optimization techniques. HH models additionally have severe non-linearities resulting in further challenges for inferring parameters in an algorithmic fashion. To address these challenges with a tractable method in high-dimensional parameter spaces, we propose using a particular Markov chain Monte Carlo (MCMC) algorithm, which has the advantage of inferring parameters in a Bayesian framework. The Bayesian approach is designed to be suitable for multimodal solutions to inverse problems. We introduce and demonstrate the method using a three-channel HH model. We then focus on the inference of nine parameters in an eight-channel HH model, which we analyze in detail. We explore how the MCMC algorithm can uncover complex relationships between inferred parameters using five injected current levels. The MCMC method provides as a result a nine-dimensional posterior distribution, which we analyze visually with solution maps or landscapes of the possible parameter sets. The visualized solution maps show new complex structures of the multimodal posteriors, and they allow for selection of locally and globally optimal value sets, and they visually expose parameter sensitivities and regions of higher model robustness. We envision these solution maps as enabling experimentalists to improve the design of future experiments, increase scientific productivity and improve on model structure and ideation when the MCMC algorithm is applied to experimental data.

A new micromagnetic method using the hybrid Monte Carlo (HMC) algorithm has been developed, and the magnetic properties at finite temperature can be simulated. In this paper, the HMC micromagnetics are applied on the polycrystalline FePt–C granular films for heat-assisted magnetic recording, where the significant parameters such as the temperature-dependent saturation magnetization $M_{s}$ (T) are measured by experiments. The Voronoi polycrystalline microstructure is built upon a regular mesh of micromagnetic cells. The magnetization (M–H) loops below the Curie temperature are calculated; they agree well with experiments at 300, 400, 500, and 600 K, respectively, which have confirmed the validity of the HMC micromagnetic method. Besides, the M–H loop at 700 K, which currently cannot be measured by experiments, can also be simulated using the HMC algorithm, which shows superparamagnetic properties.

In an effort to build seismic models that are the most consistent with multiple data sets we have applied a new probabilistic inverse technique. This method uses a Markov chain Monte Carlo (MCMC) algorithm to sample models from a prior distribution and test them against multiple data types to generate a posterior distribution. While computationally expensive, this approach has several advantages over deterministic models, notably the seamless reconciliation of different data types that constrain the model, the proper handling of both data and model uncertainties, and the ability to easily incorporate a variety of prior information, all in a straightforward, natural fashion. A real advantage of the technique is that it provides a more complete picture of the solution space. By mapping out the posterior probability density function, we can avoid simplistic assumptions about the model space and allow alternative solutions to be identified, compared, and ranked. Here we use this method to determine the crust and upper mantle structure of the Yellow Sea and Korean Peninsula region. The model is parameterized as a series of seven layers in a regular latitude‐longitude grid, each of which is characterized by thickness and seismic parameters (Vp, Vs, and density). We use surface wave dispersion and body wave traveltime data to drive the model. We find that when properly tuned (i.e., the Markov chains have had adequate time to fully sample the model space and the inversion has converged), the technique behaves as expected. The posterior model reflects the prior information at the edge of the model where there is little or no data to constrain adjustments, but the range of acceptable models is significantly reduced in data‐rich regions, producing values of sediment thickness, crustal thickness, and upper mantle velocities consistent with expectations based on knowledge of the regional tectonic setting.

We explore a general framework in Markov chain Monte Carlo (MCMC) sampling where sequential proposals are tried as a candidate for the next state of the Markov chain. This sequential-proposal framework can be applied to various existing MCMC methods, including Metropolis–Hastings algorithms using random proposals and methods that use deterministic proposals such as Hamiltonian Monte Carlo (HMC) or the bouncy particle sampler. Sequential-proposal MCMC methods construct the same Markov chains as those constructed by the delayed rejection method under certain circumstances. In the context of HMC, the sequential-proposal approach has been proposed as extra chance generalized hybrid Monte Carlo (XCGHMC). We develop two novel methods in which the trajectories leading to proposals in HMC are automatically tuned to avoid doubling back, as in the No-U-Turn sampler (NUTS). The numerical efficiency of these new methods compare favorably to the NUTS. We additionally show that the sequential-proposal bouncy particle sampler enables the constructed Markov chain to pass through regions of low target density and thus facilitates better mixing of the chain when the target density is multimodal.

A reliable assessment of the posterior uncertainties is a crucial aspect of any amplitude versus angle (AVA) inversion due to the severe ill-conditioning of this inverse problem. To accomplish this task, numerical Markov chain Monte Carlo algorithms are usually used when the forward operator is nonlinear. The downside of these algorithms is the considerable number of samples needed to attain stable posterior estimations especially in high-dimensional spaces. To overcome this issue, we assessed the suitability of Hamiltonian Monte Carlo (HMC) algorithm for nonlinear target- and interval-oriented AVA inversions for the estimation of elastic properties and associated uncertainties from prestack seismic data. The target-oriented approach inverts the AVA responses of the target reflection by adopting the nonlinear Zoeppritz equations, whereas the interval-oriented method inverts the seismic amplitudes along a time interval using a 1D convolutional forward model still based on the Zoeppritz equations. HMC uses an artificial Hamiltonian system in which a model is viewed as a particle moving along a trajectory in an extended space. In this context, the inclusion of the derivative information of the misfit function makes possible long-distance moves with a high probability of acceptance from the current position toward a new independent model. In our application, we adopt a simple Gaussian a priori distribution that allows for an analytical inclusion of geostatistical constraints into the inversion framework, and we also develop a strategy that replaces the numerical computation of the Jacobian with a matrix operator analytically derived from a linearization of the Zoeppritz equations. Synthetic and field data inversions demonstrate that the HMC is a very promising approach for Bayesian AVA inversion that guarantees an efficient sampling of the model space and retrieves reliable estimations and accurate uncertainty quantifications with an affordable computational cost.

Recent progress using geometry in the design of efficient Markov chain Monte Carlo (MCMC) algorithms have shown the effectiveness of the Fisher Riemannian structure. Furthermore, the theory of the underlying geometry of spaces of statistical models has made an important breakthrough by extending the classical theory on exponential families to their closures, the so-called extended exponential families. This paper looks at the underlying geometry of the Fisher information, in particular its limiting behaviour near boundaries, which illuminates the excellent behaviour of the corresponding geometric MCMC algorithms. Further, the paper shows how Fisher geodesics in extended exponential families smoothly attach the boundaries of extended exponential families to their relative interior. We conjecture that this behaviour could be exploited for trans-dimensional MCMC algorithms.

Timely detection of localised COVID-19 surges is essential for targeting limited health resources, yet most routine surveillance algorithms ignore spatial dependence and many Bayesian spatio-temporal models are computationally demanding. Evidence on Hamiltonian Monte Carlo (HMC) performance for outbreak detection in low- and middle-income country (LMIC) settings remains limited. We applied a Bayesian spatio-temporal hidden Markov model (HMM) to South African COVID-19 hospital mortality, comparing data-augmented Markov chain Monte Carlo (MCMC) with dynamic HMC. We conducted a retrospective ecological time-series study of in-hospital COVID- 19 deaths across 52 districts over 28 months (March 2020-June 2022), using national hospital surveillance linked to district-level health-system and population indicators. Death counts were modelled with a Poisson log-linear specification incorporating a smooth temporal trend, cyclic seasonality, spatial smoothing, and selected covariates, offset by expected deaths from admissions and a case-fatality ratio. Outbreaks were represented by a two-state HMM with latent indicators integrated out analytically and estimated via dynamic HMC in Stan. Eight candidate models were ranked using marginal likelihoods; the preferred model was re-fitted with both samplers to compare runtime, ESS, and convergence. A spatial HMM with marginalised outbreak states was strongly favoured over nonoutbreak and threshold-based alternatives. Posterior outbreak probabilities reproduced the four recognised national waves while revealing marked district-level heterogeneity, with early intense outbreaks in Western Cape and Gauteng districts and later peaks inland. Outbreaks were short-lived (mean around three months), and residual spatial risks indicated persistent excess mortality in the Eastern Cape and Free State. Dynamic HMC and MCMC yielded similar outbreak probability surfaces; however, HMC produced substantially larger ESS 1 (approximately 5,000 versus 63) and near-ideal convergence, whereas MCMC showed poor mixing. ESS per second was similar, so HMC's extra computation yielded more informative samples. A Bayesian spatio-temporal HMM fitted with dynamic HMC delivers coherent, spatially resolved outbreak probabilities and captures short-lived district-level mortality surges within broader national waves. Despite greater computational intensity, dynamic HMC offers superior mixing and convergence and is preferable for routine surveillance when adequate computing resources are available. The framework is transferable to other routinely collected surveillance data in South Africa and similar LMIC settings.

An integrated Markov chain Monte Carlo algorithm for upscaling hydrological and geochemical parameters from column to field scale

In this paper, we derive a Markov chain Monte Carlo (MCMC) algorithm supported by a neural network. In particular, we use the neural network to substitute derivative calculations made during a Metropolis adjusted Langevin algorithm (MALA) step with inexpensive neural network evaluations. Using a complex, high-dimensional blood coagulation model and a set of measurements, we define a likelihood function on which we evaluate the new MCMC algorithm. The blood coagulation model is a dynamic model, where derivative calculations are expensive and hence limit the efficiency of derivative-based MCMC algorithms. The MALA adaptation greatly reduces the time per iteration, while only slightly affecting the sample quality. We also test the new algorithm on a 2-dimensional example with a non-convex shape, a case where the MALA algorithm has a clear advantage over other state of the art MCMC algorithms. To assess the impact of the new algorithm, we compare the results to previously generated results of the MALA and the random walk Metropolis Hastings (RWMH).

There is an increasing use of Markov chain Monte Carlo (MCMC) algorithms for fitting statistical models in psychometrics, especially in situations where the traditional estimation techniques are very difficult to apply. One of the disadvantages of using an MCMC algorithm is that it is not straightforward to determine the convergence of the algorithm. Using the output of an MCMC algorithm that has not converged may lead to incorrect inferences on the problem at hand. The convergence is not one to a point, but that of the distribution of a sequence of generated values to another distribution, and hence is not easy to assess; there is no guaranteed diagnostic tool to determine convergence of an MCMC algorithm in general. This article examines the convergence of MCMC algorithms using a number of convergence diagnostics for two real data examples from psychometrics. Findings from this research have the potential to be useful to researchers using the algorithms. For both the examples, the number of iterations required (suggested by the diagnostics) to be reasonably confident that the MCMC algorithm has converged may be larger than what many practitioners consider to be safe.