Adaptive Markov Chain Monte Carlo Research Articles

Sequentially Guided MCMC Proposals for Synthetic Likelihoods and Correlated Synthetic Likelihoods

Synthetic likelihood (SL) is a strategy for parameter inference when the likelihood function is analytically or computationally intractable. In SL, the likelihood function of the data is replaced by a multivariate Gaussian density over summary statistics of the data. SL requires simulation of many replicate datasets at every parameter value considered by a sampling algorithm, such as Markov chain Monte Carlo (MCMC), making the method computationally-intensive. We propose two strategies to alleviate the computational burden. First, we introduce an algorithm producing a proposal distribution that is sequentially tuned and made conditional to data, thus it rapidly \textit{guides} the proposed parameters towards high posterior density regions. In our experiments, a small number of iterations of our algorithm is enough to rapidly locate high density regions, which we use to initialize one or several chains that make use of off-the-shelf adaptive MCMC methods. Our "guided" approach can also be potentially used with MCMC samplers for approximate Bayesian computation (ABC). Second, we exploit strategies borrowed from the correlated pseudo-marginal MCMC literature, to improve the chains mixing in a SL framework. Moreover, our methods enable inference for challenging case studies, when the posterior is multimodal and when the chain is initialised in low posterior probability regions of the parameter space, where standard samplers failed. To illustrate the advantages stemming from our framework we consider five benchmark examples, including estimation of parameters for a cosmological model and a stochastic model with highly non-Gaussian summary statistics.

Sampling by divergence minimization

We introduce a Markov Chain Monte Carlo (MCMC) method that is designed to sample from target distributions with irregular geometry using an adaptive scheme. In cases where targets exhibit non-Gaussian behavior, we propose that adaptation should be regional rather than global. Our algorithm minimizes the information projection component of the Kullback–Leibler (KL) divergence between the proposal and target distributions to encourage proposals that are distributed similarly to the regional geometry of the target. Unlike traditional adaptive MCMC, this procedure rapidly adapts to the geometry of the target’s current position as it explores the surrounding space without the need for many preexisting samples. The divergence minimization algorithms are tested on target distributions with irregularly shaped modes and we provide results demonstrating the effectiveness of our methods.

Development of a two-phase adaptive MCMC method for efficient Bayesian model updating of complex dynamic systems

The fundamental problem of Bayesian model updating is identifying the posterior probability density function (PDF) of uncertain model parameters. Markov chain Monte Carlo (MCMC) has been recognized to be potential for this problem. However, the main difficulty is that many samples need to be generated first to locate the regions of high probability. This could be time consuming for complex dynamic systems. A novel two-phase adaptive MCMC method is developed in this work to tackle this difficulty. A parameter-space search algorithm is used to do a systematic search to find output-equivalent points, which are in the neighborhood of the regions of high probability. To generate samples of uncertain model parameters in these regions for describing the shape of the posterior PDF, a Markov chain is constructed with these output-equivalent points as seeds. This is done by the proposed weighted MCMC algorithm based on the modified Metropolis-Hasting (MH) algorithm. To enhance the efficiency of the algorithm, a new proposal PDF is developed using the seeds. The covariance matrix of the proposal PDF is adaptively updated with the latest set of generated samples. Furthermore, the optimal prediction-error variance as a function of model parameters is analytical derived and its values at the seeds are adaptively used to estimate the acceptance ratio in the modified MH algorithm. The proposed weighted MCMC algorithm generates samples in the region of high probability adaptively without going through multiple levels, which are computationally demanding. The posterior uncertainties in model updating can be rigorously considered by the generated samples with their weights proportional to their posterior PDF values. To verify the proposed method, numerical simulation of a two-story shear building with the challenging locally identifiable and unidentifiable cases was conducted. Further verification was also conducted by updating a finite element model of a full-scaled coupled building system utilizing measured vibration data from field tests.

An adaptive multiple-try Metropolis algorithm

Markov chain Monte Carlo (MCMC) methods, specifically samplers based on random walks, often have difficulty handling target distributions with complex geometry such as multi-modality. We propose an adaptive multiple-try Metropolis algorithm designed to tackle such problems by combining the flexibility of multiple-proposal samplers with the user-friendliness and optimality of adaptive algorithms. We prove the ergodicity of the resulting Markov chain with respect to the target distribution using common techniques in the adaptive MCMC literature. In a Bayesian model for loss of heterozygosity in cancer cells, we find that our method outperforms traditional adaptive samplers, non-adaptive multiple-try Metropolis samplers, and various more sophisticated competing methods.

The Barker proposal: Combining robustness and efficiency in gradient-based MCMC.

There is a tension between robustness and efficiency when designing Markov chain Monte Carlo (MCMC) sampling algorithms. Here we focus on robustness with respect to tuning parameters, showing that more sophisticated algorithms tend to be more sensitive to the choice of step‐size parameter and less robust to heterogeneity of the distribution of interest. We characterise this phenomenon by studying the behaviour of spectral gaps as an increasingly poor step‐size is chosen for the algorithm. Motivated by these considerations, we propose a novel and simple gradient‐based MCMC algorithm, inspired by the classical Barker accept‐reject rule, with improved robustness properties. Extensive theoretical results, dealing with robustness to tuning, geometric ergodicity and scaling with dimension, suggest that the novel scheme combines the robustness of simple schemes with the efficiency of gradient‐based ones. We show numerically that this type of robustness is particularly beneficial in the context of adaptive MCMC, giving examples where our proposed scheme significantly outperforms state‐of‐the‐art alternatives.

Convergence assessment for Bayesian phylogenetic analysis using MCMC simulation

Abstract Posterior distributions are commonly approximated by samples produced from a Markov chain Monte Carlo (MCMC) simulation. Every MCMC simulation has to be checked for convergence, that is, that sufficiently many samples have been obtained and that these samples indeed represent the true posterior distribution. Here we develop and test different approaches for convergence assessment in phylogenetics. We analytically derive a threshold for a minimum effective sample size (ESS) of 625. We observe that only the initial sequence estimator provides robust ESS estimates for common types of MCMC simulations (autocorrelated samples, adaptive MCMC, Metropolis‐coupled MCMC). We show that standard ESS computation can be applied to phylogenetic trees if the tree samples are converted into traces of absence/presence of splits. Convergence in distribution between replicated MCMC runs can be assessed with the Kolmogorov–Smirnov test. The commonly used potential scale reduction factor (PSRF) is biased when applied to skewed posterior distribution. Additionally, we provide how the distribution of differences in split frequencies can be computed exactly akin to standard exact tests and show that it depends on the true frequency of a split. Hence, the average standard deviation of split frequencies is too simplistic and the expected difference based on the 95% quantile should be used instead to check for convergence in split frequencies. We implemented the methods described here in the open‐source R package Convenience (https://github.com/lfabreti/convenience), which allows users to easily test for convergence using output from standard phylogenetic inference software.

Parameter Estimation of Partially Observed Turbulent Systems Using Conditional Gaussian Path-Wise Sampler

Parameter estimation of complex nonlinear turbulent dynamical systems using only partially observed time series is a challenging topic. The nonlinearity and partial observations often impede using closed analytic formulae to recover the model parameters. In this paper, an exact path-wise sampling method is developed, which is incorporated into a Bayesian Markov chain Monte Carlo (MCMC) algorithm in light of data augmentation to efficiently estimate the parameters in a rich class of nonlinear and non-Gaussian turbulent systems using partial observations. This path-wise sampling method exploits closed analytic formulae to sample the trajectories of the unobserved variables, which avoid the numerical errors in the general sampling approaches and significantly increase the overall parameter estimation efficiency. The unknown parameters and the missing trajectories are estimated in an alternating fashion in an adaptive MCMC iteration algorithm with rapid convergence. It is shown based on the noisy Lorenz 63 model and a stochastically coupled FitzHugh–Nagumo model that the new algorithm is very skillful in estimating the parameters in highly nonlinear turbulent models. The model with the estimated parameters succeeds in recovering the nonlinear and non-Gaussian features of the truth, including capturing the intermittency and extreme events, in both test examples.

Uncertainty Quantification of Water Level Predictions from Radar‐based Areal Rainfall Using an Adaptive MCMC Algorithm

This study proposes an approach for the uncertainty quantification at each stage of a single hydrological process of water level predictions based on different sources of mean areal precipitation (MAP) forecasts by using an adaptive Bayesian Markov chain Monte Carlo (MCMC) approach. The MAP forecasts are derived from the McGill Algorithm for Precipitation Nowcasting by Lagrangian Extrapolation (MAPLE) system and a long short-term memory (LSTM) network. The predicted water levels at two stations in the Gangnam catchment, Seoul, South Korea, are processed with a coupled 1D/2D urban hydrological model (1D/2D-UHM) forced by MAPLE MAP forecasts and LSTM-corrected MAP forecasts of five heavy rainfall events. The proposed Bayesian approach using the delayed rejection and adaptive Metropolis (DRAM) algorithm was compared with the Metropolis-Hastings (MH) algorithm in the uncertainty estimation of Weibull distribution parameters. The uncertainty contributions of the stages and sources in the related process were analyzed, including quantitative precipitation estimation (QPE) inputs, MAP inputs and 1D/2D-UHM. The results indicate that the uncertainty contribution of the MAPLE MAP forecasting is the highest in the 3-hour forecasting time. The uncertainty contribution of the QPE input for MAPLE MAP forecasting is the smallest and that of two sources, including the LSTM-corrected MAP source, and MAP and the coupled model is more significant than that of the QPE input. This research showed that the adaptive Bayesian MCMC method using the DRAM algorithm might be a robust option in quantitative uncertainty analyses of a single hydrological process, especially for urban flood management.

Adaptive sequential Monte Carlo for posterior inference and model selection among complex geological priors

SUMMARY Bayesian model selection enables comparison and ranking of conceptual subsurface models described by spatial prior models, according to the support provided by available geophysical data. Deep generative neural networks can efficiently encode such complex spatial priors, thereby, allowing for a strong model dimensionality reduction that comes at the price of enhanced non-linearity. In this setting, we explore a recent adaptive sequential Monte Carlo (ASMC) approach that builds on annealed importance sampling (AIS); a method that provides both the posterior probability density function (PDF) and the evidence (a central quantity for Bayesian model selection) through a particle approximation. Both techniques are well suited to parallel computation and rely on importance sampling over a sequence of intermediate distributions, linking the prior and the posterior PDF. Each subsequent distribution is approximated by updating the particle weights and states, compared with the previous approximation, using a small pre-defined number of Markov chain Monte Carlo (MCMC) proposal steps. Compared with AIS, the ASMC method adaptively tunes the tempering between neighboring distributions and performs resampling of particles when the variance of the particle weights becomes too large. We evaluate ASMC using two different conceptual models and associated synthetic cross-hole ground penetrating radar tomography data. For the most challenging test case, we find that the ASMC method is faster and more reliable in locating the posterior PDF than state-of-the-art adaptive MCMC. The evidence estimates are found to be robust with respect to the choice of ASMC algorithmic variables and much less sensitive to the model proposal type than MCMC. The variance of the evidence estimates are best estimated by replication of ASMC runs, while approximations based on single runs provide comparable estimates when using a sufficient number of proposal steps in approximating each intermediate distribution.

Path differential-informed stratified MCMC and adaptive forward path sampling

Markov Chain Monte Carlo (MCMC) rendering is extensively studied, yet it remains largely unused in practice. We propose solutions to several practicability issues, opening up path space MCMC to become an adaptive sampling framework around established Monte Carlo (MC) techniques. We address non-uniform image quality by deriving an analytic target function for imagespace sample stratification. The function is based on a novel connection between variance and path differentials, allowing analytic variance estimates for MC samples, with potential uses in other adaptive algorithms outside MCMC. We simplify these estimates down to simple expressions using only quantities known in any MC renderer. We also address the issue that most existing MCMC renderers rely on bi-directional path tracing and reciprocal transport, which can be too costly and/or too complex in practice. Instead, we apply our theoretical framework to optimize an adaptive MCMC algorithm that only uses forward path construction. Notably, we construct our algorithm by adapting (with minimal changes) a full-featured path tracer into a single-path state space Markov Chain, bridging another gap between MCMC and existing MC techniques.

Markov chain Monte Carlo with neural network surrogates: application to contaminant source identification

Subsurface remediation often involves reconstruction of contaminant release history from sparse observations of solute concentration. Markov Chain Monte Carlo (MCMC), the most accurate and general method for this task, is rarely used in practice because of its high computational cost associated with multiple solves of contaminant transport equations. We propose an adaptive MCMC method, in which a transport model is replaced with a fast and accurate surrogate model in the form of a deep convolutional neural network (CNN). The CNN-based surrogate is trained on a small number of the transport model runs based on the prior knowledge of the unknown release history. Thus reduced computational cost allows one to reduce the sampling error associated with construction of the approximate likelihood function. As all MCMC strategies for source identification, our method has an added advantage of quantifying predictive uncertainty and accounting for measurement errors. Our numerical experiments demonstrate the accuracy comparable to that of MCMC with the forward transport model, which is obtained at a fraction of the computational cost of the latter.

Randomized reduced forward models for efficient Metropolis–Hastings MCMC, with application to subsurface fluid flow and capacitance tomography

Bayesian modelling and computational inference by Markov chain Monte Carlo (MCMC) is a principled framework for large-scale uncertainty quantification, though is limited in practice by computational cost when implemented in the simplest form that requires simulating an accurate computer model at each iteration of the MCMC. The delayed acceptance Metropolis–Hastings MCMC leverages a reduced model for the forward map to lower the compute cost per iteration, though necessarily reduces statistical efficiency that can, without care, lead to no reduction in the computational cost of computing estimates to a desired accuracy. Randomizing the reduced model for the forward map can dramatically improve computational efficiency, by maintaining the low cost per iteration but also avoiding appreciable loss of statistical efficiency. Randomized maps are constructed by a posteriori adaptive tuning of a randomized and locally-corrected deterministic reduced model. Equivalently, the approximated posterior distribution may be viewed as induced by a modified likelihood function for use with the reduced map, with parameters tuned to optimize the quality of the approximation to the correct posterior distribution. Conditions for adaptive MCMC algorithms allow practical approximations and algorithms that have guaranteed ergodicity for the target distribution. Good statistical and computational efficiencies are demonstrated in examples of calibration of large-scale numerical models of geothermal reservoirs and electrical capacitance tomography.

Posterior inference for sparse hierarchical non-stationary models

Gaussian processes are valuable tools for non-parametric modelling, where typically an assumption of stationarity is employed. While removing this assumption can improve prediction, fitting such models is challenging. Hierarchical models are constructed based on Gaussian Markov random fields with stochastic spatially varying parameters. Importantly, this allows for non-stationarity while also addressing the computational burden through a sparse banded representation of the precision matrix. In this setting, efficient Markov chain Monte Carlo (MCMC) sampling is challenging due to the strong coupling a posteriori of the parameters and hyperparameters. Three adaptive MCMC schemes are developed and compared making use of banded matrix operations for faster inference. Furthermore, a novel extension to higher dimensional input spaces is proposed through an additive structure that retains the flexibility and scalability of the model, while also inheriting interpretability from the additive approach. A thorough assessment of the efficiency and accuracy of the methods in nonstationary settings is presented for both simulated experiments and a computer emulation problem.

MCMC Techniques for Parameter Estimation of ODE Based Models in Systems Biology

Ordinary differential equation systems (ODEs) are frequently used for dynamical system modelling in many science fields such as economics, physics, engineering, and systems biology. A special challenge in systems biology is that ODE systems typically contain kinetic rate parameters, which are unknown and have to be estimated from data. However, non-linearity of ODE systems together with noise in the data raise severe identifiability issues. Hence, Markov Chain Monte Carlo (MCMC) approaches have been frequently used to estimate posterior distributions of rate parameters. However, designing a good MCMC sampler for high dimensional and multi-modal parameter distributions remains a challenging task. Here we performed a systematic comparison of different MCMC techniques for this purpose using five public domain models. The comparison included Metropolis-Hastings, parallel tempering MCMC, adaptive MCMC and parallel adaptive MCMC. In conclusion, we found specifically parallel adaptive MCMC to produce superior parameter estimates while benefitting from inclusion of our suggested informative Bayesian priors for rate parameters and noise variance.

A Bayesian Approach to Real-Time Dynamic Parameter Estimation Using Phasor Measurement Unit Measurement

In this work, we develop a polynomial-chaos-expa- nsion (PCE)-based approach for decentralized dynamic parameter estimation through Bayesian inference. Using this approach, the non-Gaussian distribution of the inverted parameters is obtained. More specifically, we first represent the decentralized generator model with the PCE-based surrogate. This surrogate allows us to efficiently evaluate the time-consuming dynamic solver at parameter values through Metropolis-Hastings (M-H)-based Markov chain Monte Carlo (MCMC). Then, we propose a two-stage hybrid Markov chain Monte Carlo (MCMC) to recover a posteriori distribution of the decentralized generator model parameters. In the first stage, we use the gradient-enhanced Langevin MCMC algorithm to characterize an intermediate posterior parameter distribution. This algorithm is computationally scalable to the high-dimensional parameter space. Based on the intermediate posterior distribution, during the second stage, we use the adaptive MCMC algorithm to fine-tune the strong correlations between the parameters. Finally, the fully recovered a posterior distribution is obtained in the end. The simulation results show that the proposed PCE-based hybrid MCMC algorithm can accurately and efficiently estimate the high-dimensional generator dynamic model parameters with full probabilistic distribution provided.

A low-complexity algorithm for the joint antenna selection and user scheduling in multi-cell multi-user downlink massive MIMO systems

The massive MIMO (multiple-input multiple-output) technology plays a key role in the next-generation (5G) wireless communication systems, which are equipped with a large number of antennas at the base station (BS) of a network to improve cell capacity for network communication systems. However, activating a large number of BS antennas needs a large number of radio-frequency (RF) chains that introduce the high cost of the hardware and high power consumption. Our objective is to achieve the optimal combination subset of BS antennas and users to approach the maximum cell capacity, simultaneously. However, the optimal solution to this problem can be achieved by using an exhaustive search (ES) algorithm by considering all possible combinations of BS antennas and users, which leads to the exponential growth of the combinatorial complexity with the increasing of the number of BS antennas and active users. Thus, the ES algorithm cannot be used in massive MIMO systems because of its high computational complexity. Hence, considering the trade-off between network performance and computational complexity, we proposed a low-complexity joint antenna selection and user scheduling (JASUS) method based on Adaptive Markov Chain Monte Carlo (AMCMC) algorithm for multi-cell multi-user massive MIMO downlink systems. AMCMC algorithm is helpful for selecting combination subset of antennas and users to approach the maximum cell capacity with consideration of the multi-cell interference. Performance analysis and simulation results show that AMCMC algorithm performs extremely closely to ES-based JASUS algorithm. Compared with other algorithms in our experiments, the higher cell capacity and near-optimal system performance can be obtained by using the AMCMC algorithm. At the same time, the computational complexity is reduced significantly by combining with AMCMC.

Bayesian modelling of nonlinear negative binomial integer-valued GARCHX models

This study focuses on modelling dengue cases in northeastern Thailand through two meteorological covariates: cumulative rainfall and average maximum temperature. We propose two nonlinear integer-valued GARCHX models (Markov switching and threshold specification) with a negative binomial distribution, as they take into account the stylized features of weekly dengue haemorrhagic fever cases, which contain nonlinear dynamics, lagged dependence, overdispersion, consecutive zeros and asymmetric effects of meteorological covariates. We conduct parameter estimation and one-step-ahead forecasting for two proposed models based on Bayesian Markov chain Monte Carlo (MCMC) methods. A simulation study illustrates that the adaptive MCMC sampling scheme performs well. The empirical results offer strong support for the Markov switching integer-valued GARCHX model over its competitors via Bayes factor and deviance information criterion. We also provide one-step-ahead forecasting based on the prediction interval that offers a useful early warning signal of outbreak detection.

Adaptive Markov chain Monte Carlo algorithms for Bayesian inference: recent advances and comparative study

Condition assessments of structures require prediction models such as empirical model and numerical simulation model. Generally, these prediction models have model parameters to be estimated from experimental data. Bayesian inference is the formal statistical framework to estimate the model parameters and their uncertainties. As a result, uncertainties associated with the model and measurement can be accounted for decision making. Markov Chain Monte Carlo (MCMC) algorithms have been widely employed. However, there still remain some implementation issues from the inappropriate selection of the proposal mechanism in Markov chain. Since the posterior density for a given problem is often problem-dependent and unknown, users require a trial-and-error approach to select and tune optimal proposal mechanism. To relieve this difficulty, various adaptive MCMC algorithms have been recently appeared. Users must understand their mechanism and limitations before applying the algorithms to their problems. However, there is no comprehensive work to provide detailed exposition and their performance comparison together. This study aims to bring together different adaptive MCMC algorithms with the goal of providing their mechanisms and evaluating their performances through comparative study. Three algorithms are chosen as the representative proposal mechanism. From comparative studies, the discussions were drawn in terms of performances, simplicity and computational costs for less-experienced users.

Robust and Fast Markov Chain Monte Carlo Sampling of Diffusion MRI Microstructure Models.

In diffusion MRI analysis, advances in biophysical multi-compartment modeling have gained popularity over the conventional Diffusion Tensor Imaging (DTI), because they can obtain a greater specificity in relating the dMRI signal to underlying cellular microstructure. Biophysical multi-compartment models require a parameter estimation, typically performed using either the Maximum Likelihood Estimation (MLE) or the Markov Chain Monte Carlo (MCMC) sampling. Whereas, the MLE provides only a point estimate of the fitted model parameters, the MCMC recovers the entire posterior distribution of the model parameters given in the data, providing additional information such as parameter uncertainty and correlations. MCMC sampling is currently not routinely applied in dMRI microstructure modeling, as it requires adjustment and tuning, specific to each model, particularly in the choice of proposal distributions, burn-in length, thinning, and the number of samples to store. In addition, sampling often takes at least an order of magnitude, more time than non-linear optimization. Here we investigate the performance of the MCMC algorithm variations over multiple popular diffusion microstructure models, to examine whether a single, well performing variation could be applied efficiently and robustly to many models. Using an efficient GPU-based implementation, we showed that run times can be removed as a prohibitive constraint for the sampling of diffusion multi-compartment models. Using this implementation, we investigated the effectiveness of different adaptive MCMC algorithms, burn-in, initialization, and thinning. Finally we applied the theory of the Effective Sample Size, to the diffusion multi-compartment models, as a way of determining a relatively general target for the number of samples needed to characterize parameter distributions for different models and data sets. We conclude that adaptive Metropolis methods increase MCMC performance and select the Adaptive Metropolis-Within-Gibbs (AMWG) algorithm as the primary method. We furthermore advise to initialize the sampling with an MLE point estimate, in which case 100 to 200 samples are sufficient as a burn-in. Finally, we advise against thinning in most use-cases and as a relatively general target for the number of samples, we recommend a multivariate Effective Sample Size of 2,200.

An adaptive Kriging surrogate method for efficient joint estimation of hydraulic and biochemical parameters in reactive transport modeling

Groundwater reactive transport models that consider the coupling of hydraulic and biochemical processes are vital tools for predicting the fate of groundwater contaminants and effective groundwater management. The models involve a large number of parameters whose specification greatly affects the model performance. Thus model parameters calibration is crucial to its successful application. The Bayesian inference framework implemented by Markov chain Monte Carlo (MCMC) sampling provides a comprehensive framework to estimate the model parameters. However, its application is hampered by the large computational requirements caused by repeated evaluations of the model in MCMC sampling. This study develops an adaptive Kriging-based MCMC method to overcome the bottleneck of Bayesian inference by replacing the simulation model with a computationally inexpensive Kriging surrogate model. In the adaptive Kriging-based MCMC method, instead of constructing a globally accurate surrogate of the simulation model, we sequentially build a locally accurate surrogate with an iterative refinement to the high probability regions. The performance of the proposed method is demonstrated using a synthetic groundwater reactive transport model for describing sequential Kinetic degradation of Tetrachloroethene (PCE), whose hydraulic and biochemical parameters are jointly estimated. The results suggest that the adaptive Kriging-based MCMC method is able to achieve an accurate Bayesian inference with a hundredfold reduction in the computational cost compared to the conventional MCMC method.