Markov Chain Monte Carlo Sampling Method Research Articles

Identifiability and convergence behavior for Markov chain Monte Carlo using multivariate probit models

Multivariate probit models have been popularly utilized to analysis multivariate ordinal data. However, the identifiable multivariate probit models entail the covariance matrix for the underlying multivariate normal variables to be a correlation matrix, which brings a rigorous task to conduct efficient statistical analysis. Parameter expansion to make the identifiable model to be non-identifiable has been inevitably explored. However, the effect of the expanded parameters on the convergence of Markov chain Monte Carlo (MCMC) is seldomly investigated; in addition, the comparison of MCMC developed based on the identifiable model and that based on the non-identifiable model is hardly ever explored, especially for data with large sample sizes. In this paper, we conduct a thorough investigation to illustrate the effect of the expanded parameters on the convergence of MCMC and compare the behavior of MCMC between the identifiable and non-identifiable models. Our investigation provides a practical guide regarding the construction of non-identifiable models and development of corresponding MCMC sampling methods. We conduct our investigation using simulation studies and present an application using data from the Russia Longitudinal Monitoring Survey-Higher School of Economics (RLMS-HSE) study.

Bayesian analysis of multivariate mixed longitudinal ordinal and continuous data

SummaryMultivariate longitudinal ordinal and continuous data exist in many scientific fields. However, it is a rigorous task to jointly analyse them due to the complicated correlated structures of those mixed data and the lack of a multivariate distribution. The multivariate probit model, assuming there is a multivariate normal latent variable for each multivariate ordinal data, becomes a natural modeling choice for longitudinal ordinal data especially for jointly analysing with longitudinal continuous data. However, the identifiable multivariate probit model requires the variances of the latent normal variables to be fixed at 1, thus the joint covariance matrix of the latent variables and the continuous multivariate normal variables is restricted at some of the diagonal elements. This constrains to develop both the classical and Bayesian methods to analyse mixed ordinal and continuous data. In this investigation, we proposed three Markov chain Monte Carlo (MCMC) methods: Metropolis–Hastings within Gibbs algorithm based on the identifiable model, and a Gibbs sampling algorithm and parameter‐expanded data augmentation based on the constructed non‐identifiable model. Through simulation studies and a real data application, we illustrated the performance of these three methods and provided an observation of using non‐identifiable model to develop MCMC sampling methods.

A Bayesian level set method for identifying subsurface geometries and rheological properties in Stokes flow

SUMMARY We aim to simultaneously infer the shape of subsurface structures and material properties such as density or viscosity from surface observations. Modelling mantle flow using incompressible instantaneous Stokes equations, the problem is formulated as an infinite-dimensional Bayesian inverse problem. Subsurface structures are described as level sets of a smooth auxiliary function, allowing for geometric flexibility. As inverting for subsurface structures from surface observations is inherently challenging, knowledge of plate geometries from seismic images is incorporated into the prior probability distributions. The posterior distribution is approximated using a dimension-robust Markov-chain Monte Carlo sampling method, allowing quantification of uncertainties in inferred parameters and shapes. The effectiveness of the method is demonstrated in two numerical examples with synthetic data. In a model with two higher-density sinkers, their shape and location are inferred with moderate uncertainty, but a trade-off between sinker size and density is found. The uncertainty in the inferred is significantly reduced by combining horizontal surface velocities and normal traction data. For a more realistic subduction problem, we construct tailored level-set priors, representing “seismic” knowledge and infer subducting plate geometry with their uncertainty. A trade-off between thickness and viscosity of the plate in the hinge zone is found, consistent with earlier work.

Scalable and portable computational framework enabling online probabilistic remaining useful life (RUL) estimation

This work demonstrates a framework that enables online prognostics in adhesive joints by estimating the real-time probabilistic remaining useful life (RUL) using ANNs based hybrid physics models and vectorized sequential Monte Carlo (SMC) simulations. The framework is developed by integrating the physics-based damage degradation model and uncertainty quantification (UQ) techniques to estimate both probabilistic fatigue failure life and RUL. The fatigue damage growth (FDG) simulator, a hybrid surrogate model that simulates real-time fatigue degradation in adhesive joints, is used. In the initial set of results, the generalized framework is validated by estimating the probabilistic fatigue failure life using two UQ methods: Markov Chain Monte Carlo (MCMC) and SMC method. The computational results are successfully compared against experimental data. The conventional MCMC sampling methods are inherently serial, which limits the exploitation of the computational speed-up provided by the FDG simulator and hinders the real-time life predictions. The SMC method quantifies the uncertainties by parallelizing the sampling process, significantly reducing computational time and enabling real-time prediction. Next, the generalized framework is used to estimate probabilistic RUL from the fatigue crack propagation data. The parallel SMC method showed very good speedup compared to the MCMC method. To further enhance the computational speed-up with SMC method, vectorized FDG simulations are introduced into the framework and good scalability is achieved. Finally, the portability of the framework is demonstrated by deploying it on the portable Raspberry Pi cluster.

Bayesian Analysis of Multivariate Longitudinal Ordinal Data Using Multiple Multivariate Probit Models

Multivariate longitudinal ordinal data are often involved in longitudinal studies with each individual having more than one longitudinal ordinal measure. However, due to complicated correlation structures within each individual and no explicit likelihood functions, analyzing multivariate longitudinal ordinal data is quite challenging. In this paper, Markov chain Monte Carlo (MCMC) sampling methods are developed to analyze multivariate longitudinal ordinal data by extending multivariate probit (MVP) models for univariate longitudinal ordinal data to multiple multivariate probit models (MMVP) for multivariate longitudinal ordinal data. The identifiable MVP models require the covariance matrix of the latent multivariate normal variables underlying the longitudinal ordinal variables to be a correlation matrix, thus a Metropolis-Hastings (MH) algorithm is usually necessitated, which brings a rigorous task to develop efficient MCMC sampling methods. In contrast to the identifiable MVP models, the non-identifiable MVP models can be constructed to circumvent a MH algorithm to sample a correlation matrix by a Gibbs sampling to sample a covariance matrix, and hence improve the mixing and convergence of the MCMC components. Therefore, both the identifiable MMVP models and the non-identifiable MMVP models for multivariate longitudinal ordinal data are presented, and their corresponding MCMC sampling methods are developed. The performances of these methods are illustrated through simulation studies and an application using data from the Russia Longitudinal Monitoring Survey-Higher School of Economics (RLMS-HSE).

Research on an Optimized Evaluation Method of the Bearing Capacity of Reinforced Concrete Beam Based on the Bayesian Theory

An optimized evaluation method of the bearing capacity of reinforced concrete beam based on the Bayesian theory was proposed in this paper. This evaluation method optimized the traditional Markov Chain-Monte Carlo (MCMC) sampling method, and proposed an improved Metropolis-Hastings (MH) sampling method and a transitive MCMC (TMCMC) sampling method based on the MCMC theory. These two derived sampling methods solved the problem that the traditional MCMC algorithm makes it difficult to achieve convergence when the number of modified parameters is large. Therefore, on the basis of obtaining the measured sample information and the prior information of uncertain parameters, this paper first used multiple "model components" to form a model sample, then carried out a sensitivity analysis based on the relevant response indicators and selected the key parameters that had a great impact on the bearing capacity, carried out static load tests, and extracted and analyzed the experimental data. Then, based on a large amount of analysis data, the improved MH sampling method and TMCMC sampling method were used to establish a posterior probability distribution database. Finally, multiple posterior probability distributions were used to identify and predict the bearing capacity. The results showed that the method was feasible and effective for the evaluation of the bearing capacity of reinforced concrete beam.

A data-driven optimization method for coarse-graining gene regulatory networks

One major challenge in systems biology is to understand how various genes in a gene regulatory network (GRN) collectively perform their functions and control network dynamics. This task becomes extremely hard to tackle in the case of large networks with hundreds of genes and edges, many of which have redundant regulatory roles and functions. The existing methods for model reduction usually require the detailed mathematical description of dynamical systems and their corresponding kinetic parameters, which are often not available. Here, we present a data-driven method for coarse-graining large GRNs, named SacoGraci, using ensemble-based mathematical modeling, dimensionality reduction, and gene circuit optimization by Markov Chain Monte Carlo methods. SacoGraci requires network topology as the only input and is robust against errors in GRNs. We benchmark and demonstrate its usage with synthetic, literature-based, and bioinformatics-derived GRNs. We hope SacoGraci will enhance our ability to model the gene regulation of complex biological systems.

Bayesian Graph Convolutional Neural Networks via Tempered MCMC

Deep learning models, such as convolutional neural networks, have long been applied to image and multi-media tasks, particularly those with structured data. More recently, there has been more attention to unstructured data that can be represented via graphs. These types of data are often found in health and medicine, social networks, and research data repositories. Graph convolutional neural networks have recently gained attention in the field of deep learning that takes advantage of graph-based data representation with automatic feature extraction via convolutions. Given the popularity of these methods in a wide range of applications, robust uncertainty quantification is vital. This remains a challenge for large models and unstructured datasets. Bayesian inference provides a principled approach to uncertainty quantification of model parameters for deep learning models. Although Bayesian inference has been used extensively elsewhere, its application to deep learning remains limited due to the computational requirements of the Markov Chain Monte Carlo (MCMC) methods. Recent advances in parallel computing and advanced proposal schemes in MCMC sampling methods has opened the path for Bayesian deep learning. In this paper, we present Bayesian graph convolutional neural networks that employ tempered MCMC sampling with Langevin-gradient proposal distribution implemented via parallel computing. Our results show that the proposed method can provide accuracy similar to advanced optimisers while providing uncertainty quantification for key benchmark problems.

Evidence accumulation models with R: A practical guide to hierarchical Bayesian methods

Evidence accumulation models are a useful tool to allow researchers to investigate the latent cognitive variables that underlie response time and response accuracy. However, applying evidence accumulation models can be difficult because they lack easily computable forms. Numerical methods are required to determine the parameters of evidence accumulation that best correspond to the fitted data. When applied to complex cognitive models, such numerical methods can require substantial computational power which can lead to infeasibly long compute times. In this paper, we provide efficient, practical software and a step-by-step guide to fit evidence accumulation models with Bayesian methods. The software, written in C++, is provided in an R package: 'ggdmc'. The software incorporates three important ingredients of Bayesian computation, (1) the likelihood functions of two common response time models, (2) the Markov chain Monte Carlo (MCMC) algorithm (3) a population-based MCMC sampling method. The software has gone through stringent checks to be hosted on the Comprehensive R Archive Network (CRAN) and is free to download. We illustrate its basic use and an example of fitting complex hierarchical Wiener diffusion models to four shooting-decision data sets.

Direct parameter inference from global EoR signal with Bayesian statistics

ABSTRACT In the observation of sky-averaged $\mathrm{H\, \small{I}}$ signal from Epoch of Reionization (EoR), model parameter inference can be a computation-intensive work, which makes it hard to perform a direct one-stage model parameter inference by using Markov Chain Monte Carlo (MCMC) sampling method in Bayesian framework. Instead, a two-stage inference is usually used, i.e. the parameters of some characteristic points on the EoR spectrum model are first estimated, which are then used as the input to estimate physical model parameters further. However, some previous works had noticed that this kind of method could bias results, and it could be meaningful to answer the question of whether it is feasible to perform direct one-stage MCMC sampling and obtain unbiased physical model parameter estimations. In this work, we studied this problem and confirmed the feasibility. We find that unbiased estimations to physical model parameters can be obtained with a one-stage direct MCMC sampling method. We also study the influence of some factors that should be considered in practical observations to model parameter inference. We find that a very tiny amplifier gain calibration error (10−5 relative error) with complex spectral structures can significantly bias the parameter estimation; the frequency-dependent antenna beam and geographical position can also influence the results, so that should be carefully handled.

Neural Langevin Dynamical Sampling

Sampling technique is one of the asymptotically unbiased estimation approaches for inference in Bayesian probabilistic models. Markov chain Monte Carlo (MCMC) is a kind of sampling methods, which is widely used in the inference of complex probabilistic models. However, current MCMC methods can incur high autocorrelation of samples, which means that the samples generated by MCMC samplers are far from independent. In this paper, we propose the neural networks Langevin Monte Carlo (NNLMC) which makes full use of the flexibility of neural networks and the high efficiency of the Langevin dynamics sampling to construct a new MCMC sampling method. We propose the new update function to generate samples and employ appropriate loss functions to improve the performance of NNLMC during the process of sampling. We evaluate our method on a large diversity of challenging distributions and real datasets. Our results show that NNLMC is able to sample from the target distribution with low autocorrelation and rapid convergence, and outperforms the state-of-the-art MCMC samplers.

Inference of unexploded ordnance (UXO) by probabilistic inversion of magnetic data

SUMMARY Magnetic modelling of unexploded ordnance (UXO) is a well-documented method used to interpret magnetic anomalies occurring in UXO excavation surveys. By treating UXO as a ferrous spheroidal object, the induced dipole moment can be estimated by approximation of UXO characteristics such as shape, size and orientation. Inversion of magnetic data with respect to UXO requires one to solve the equation for the induced dipole moment, while also determining the location and orientation of the object. This is a highly nonlinear, non-unique problem, where many solutions often are present, which make it difficult for standard inversion methods, such as linearized approaches and maximum likelihood estimators, in assessing uncertainties and correlations in estimated model parameters that often result in an incomplete solution. In this study, we treat the problem concerning magnetic UXO inversion by a probabilistic approach using Markov chain Monte Carlo (McMC) sampling. To deal with the potential multimodality, a combination of two well-known McMC sampling methods is employed in a single-chain approach: the extended Metropolis algorithm is used for efficient local sampling, and the Gibbs sampler is used to help exploring the possible multimodal density of the posterior. By adding a Gibbs step we significantly increase the efficiency of the extended Metropolis, making it viable to use as a single-chain sampler for this problem. We refer to the algorithm as the Gibbs-within-Metropolis algorithm. We test and compare the proposed algorithm to a multichain McMC parallel tempering setup using the extended Metropolis algorithm. We then present synthetic cases and a real data case, where the algorithm is used. We demonstrate how the probabilistic approach allows full inference of the parameters describing a UXO while also including the possible presence of remanent magnetization.

A transport-based multifidelity preconditioner for Markov chain Monte Carlo

Markov chain Monte Carlo (MCMC) sampling of posterior distributions arising in Bayesian inverse problems is challenging when evaluations of the forward model are computationally expensive. Replacing the forward model with a low-cost, low-fidelity model often significantly reduces computational cost; however, employing a low-fidelity model alone means that the stationary distribution of the MCMC chain is the posterior distribution corresponding to the low-fidelity model, rather than the original posterior distribution corresponding to the high-fidelity model. We propose a multifidelity approach that combines, rather than replaces, the high-fidelity model with a low-fidelity model. First, the low-fidelity model is used to construct a transport map that deterministically couples a reference Gaussian distribution with an approximation of the low-fidelity posterior. Then, the high-fidelity posterior distribution is explored using a non-Gaussian proposal distribution derived from the transport map. This multifidelity “preconditioned” MCMC approach seeks efficient sampling via a proposal that is explicitly tailored to the posterior at hand and that is constructed efficiently with the low-fidelity model. By relying on the low-fidelity model only to construct the proposal distribution, our approach guarantees that the stationary distribution of the MCMC chain is the high-fidelity posterior. In our numerical examples, our multifidelity approach achieves significant speedups compared with single-fidelity MCMC sampling methods.

Basis-constrained Bayesian Markov-chain Monte Carlo difference inversion for geoelectrical monitoring of hydrogeologic processes

Bayesian Markov-chain Monte Carlo (McMC) techniques are increasingly being used in geophysical estimation of hydrogeologic processes due to their ability to produce multiple estimates that enable comprehensive assessment of uncertainty. Standard McMC sampling methods can, however, become computationally intractable for spatially distributed, high-dimensional problems. We have developed a novel basis-constrained Bayesian McMC difference inversion framework for time-lapse geophysical imaging. The strategy parameterizes the Bayesian inversion model space in terms of sparse, hydrologic-process-tuned bases, leading to dimensionality reduction while accounting for the physics of the target hydrologic process. We evaluate the algorithm on cross-borehole electrical resistivity tomography (ERT) field data acquired during a heat-tracer experiment. We validate the ERT-estimated temperatures with direct temperature measurements at two locations on the ERT plane. We also perform the inversions using the conventional smoothness-constrained inversion (SCI). Our approach estimates the heat plumes without excessive smoothing in contrast with the SCI thermograms. We capture most of the validation temperatures within the 90% confidence interval of the mean. Accounting for the physics of the target process allows the detection of small temperature changes that are undetectable by the SCI. Performing the inversion in the reduced-dimensional model space results in significant gains in computational cost.

Langevin-gradient parallel tempering for Bayesian neural learning

Bayesian inference provides a rigorous approach for neural learning with knowledge representation via the posterior distribution that accounts for uncertainty quantification. Markov Chain Monte Carlo (MCMC) methods typically implement Bayesian inference by sampling from the posterior distribution. This not only provides point estimates of the weights, but the ability to propagate and quantify uncertainty in decision making. However, these techniques face challenges in convergence and scalability, particularly in settings with large datasets and neural network architectures. This paper addresses these challenges in two ways. First, parallel tempering MCMC sampling method is used to explore multiple modes of the posterior distribution and implemented in multi-core computing architecture. Second, we make within-chain sampling scheme more efficient by using Langevin gradient information for creating Metropolis–Hastings proposal distributions. We demonstrate the techniques using time series prediction and pattern classification applications. The results show that the method not only improves the computational time, but provides better decision making capabilities when compared to related methods.

Bayesian inverse problems for recovering coefficients of two scale elliptic equations

We consider the Bayesian inverse homogenization problem of recovering the locally periodic two scale coefficient of a two scale elliptic equation, given limited noisy information on the solution. We consider both the uniform and the Gaussian prior probability measures. We use the two scale homogenized equation whose solution contains the solution of the homogenized equation which describes the macroscopic behaviour, and the corrector which encodes the microscopic behaviour. We approximate the posterior probability by a probability measure determined by the solution of the two scale homogenized equation. We show that the Hellinger distance of these measures converges to zero when the microscale converges to zero, and establish an explicit convergence rate when the solution of the two scale homogenized equation is sufficiently regular. Sampling the posterior measure by Markov Chain Monte Carlo (MCMC) method, instead of solving the two scale equation using fine mesh for each proposal with extremely high cost, we can solve the macroscopic two scale homogenized equation. Although this equation is posed in a high dimensional tensorized domain, it can be solved with essentially optimal complexity by the sparse tensor product finite element method, which reduces the computational complexity of the MCMC sampling method substantially. We show numerically that observations on the macrosopic behaviour alone are not sufficient to infer the microstructure. We need also observations on the corrector. Solving the two scale homogenized equation, we get both the solution to the homogenized equation and the corrector. Thus our method is particularly suitable for sampling the posterior measure of two scale coefficients.

Performance assessment of a commonly used “accumulation and wash-off” model from long-term continuous road runoff turbidity measurements

The suitability of a commonly used accumulation and wash-off model for continuous modelling of urban runoff contamination was evaluated based on 11-month turbidity and flow-rate records from an urban street. Calibration and uncertainty analysis were performed using a Markov Chain Monte-Carlo sampling method for both suspended solids loads (discharge rates) and concentration modelling. Selected models failed at replicating suspended solids concentration over the complete monitoring period. The studied dataset indeed suggests that the accumulation process is rather unpredictable and cannot be satisfactorily represented with usual accumulation models unless short periods are considered. Regarding suspended solid loads modelling, noticeably better performance was achieved, but similar results could as well be obtained with much simpler constant concentration models. Unless providing very accurate estimates of concentrations in runoff, accounting for their temporal variability during rain events may therefore not always be necessary for pollutant loads modelling, as loads are in fact mostly explained by runoff volumes.

Bayesian localization of an unknown number of ocean acoustic sources

This paper considers localizing an unknown number of ocean acoustic sources when properties of the environment are poorly known. A Bayesian formulation is developed in which environmental parameters, noise statistics, and the number, locations, and complex spectra (amplitudes and phases) of multiple sources are considered unknown random variables constrained by acoustic data and prior information. The number of sources is determined during a burn-in stage by minimizing the Bayesian information criterion using hybrid optimization with an efficient source birth/death scheme. Optimal estimates and marginal posterior probability distributions for source locations are computed employing a variety of sampling approaches. Environmental properties and source locations and are treated as explicit parameters and marginalized using Markov-chain Monte Carlo sampling methods. In particular, environmental parameters are treated using Metropolis-Hastings sampling applied efficiently in a principal-component space, and source locations are treated using Gibbs sampling since the corresponding conditional probability distributions can be computed efficiently using normal-mode methods. Source and noise spectra are sampled implicitly by applying analytic maximum-likelihood solutions expressed in terms of the explicit parameters. This represents an empirical Bayesian approximation within a hierarchical formulation, and significantly reduces the dimensionality and improves sampling efficiency in the inversion.

Emcee: The MCMC Hammer

ABSTRACTWe introduce a stable, well tested Python implementation of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) proposed by Goodman & Weare (). The code is open source and has already been used in several published projects in the astrophysics literature. The algorithm behind emcee has several advantages over traditional MCMC sampling methods and it has excellent performance as measured by the autocorrelation time (or function calls per independent sample). One major advantage of the algorithm is that it requires hand-tuning of only 1 or 2 parameters compared to ∼N2 for a traditional algorithm in an N-dimensional parameter space. In this document, we describe the algorithm and the details of our implementation. Exploiting the parallelism of the ensemble method, emcee permits any user to take advantage of multiple CPU cores without extra effort. The code is available online at http://dan.iel.fm/emcee under the GNU General Public License v2.

Markov Chain Monte Carlo versus Importance Sampling in Bayesian Inference of the GARCH Model

Usually, the Bayesian inference of the GARCH model is preferably performed by the Markov Chain Monte Carlo (MCMC) method. In this study, we also take an alternative approach to the Bayesian inference by the importance sampling. Using a multivariate Student's t-distribution that approximates the posterior density of the Bayesian inference, we compare the perfor- mance of the MCMC and importance sampling methods. The overall performance can be measured in terms of statistical errors obtained for the same size of Monte Carlo data. The Bayesian inference of the GARCH model is performed by the MCMC method implemented by the Metropolis-Hastings algorithm and the importance sampling method for artificial return data and stock return data. We find that the statistical errors of the GARCH parameters from the importance sampling are smaller than or comparable to those obtained from the MCMC method. Therefore we conclude that the importance sampling method can also be applied effectively for the Bayesian inference of the GARCH model as an alternative method to the MCMC method.