Metropolis-Hastings Markov Chain Monte Carlo Research Articles

A Bayesian approach to identifying the role of hospital structure and staff interactions in nosocomial transmission of SARS-CoV-2.

Nosocomial infections threaten patient safety, and were widely reported during the COVID-19 pandemic. Effective hospital infection control requires a detailed understanding of the role of different transmission pathways, yet these are poorly quantified. Using patient and staff data from a large UK hospital, we demonstrate a method to infer unobserved epidemiological event times efficiently and disentangle the infectious pressure dynamics by ward. A stochastic individual-level, continuous-time state-transition model was constructed to model transmission of SARS-CoV-2, incorporating a dynamic staff-patient contact network as time-varying parameters. A Metropolis-Hastings Markov chain Monte Carlo (MCMC) algorithm was used to estimate transmission rate parameters associated with each possible source of infection, and the unobserved infection and recovery times. We found that the total infectious pressure exerted on an individual in a ward varied over time, as did the primary source of transmission. There was marked heterogeneity between wards; each ward experienced unique infectious pressure over time. Hospital infection control should consider the role of between-ward movement of staff as a key infectious source of nosocomial infection for SARS-CoV-2. With further development, this method could be implemented routinely for real-time monitoring of nosocomial transmission and to evaluate interventions.

Source term estimation for continuous plume dispersion in Fusion Field Trial-07: Bayesian inference probability adjoint inverse method

In scenarios involving sudden releases of unidentified gases or concealed pollution emergencies, source control emerges as a critical procedure to safeguard residential air quality. Appropriate inverse source tracking methodology depending on diverse measurement data could be utilized to promptly identify pollutant source parameters. In this study, source term estimation (STE) method, i.e., jointly combining probability adjoint method with the Bayesian inference method, has been proposed. General form of the pollutant inverse transport equation was firstly established. Subsequently, the pollution source information, assumed from single continuous point releases during Fusion Field Trials 2007 under an unsteady wind field, was identified using the Bayesian inference probability adjoint inverse method. Metropolis-Hastings Markov Chain Monte Carlo (MH-MCMC) and Differential Evolution Markov Chain Monte Carlo (DE-MCMC) were then compared as sampling methods for Bayesian inference. Results indicated that the DE-MCMC algorithm has superior convergence and could present higher accuracy of pollutant source information than that of MH-MCMC algorithm, particularly for highly nonlinear and multi-modal distribution systems. Furthermore, the integration of Union standard Adjoint Location Probability (UALP) as prior information into the Bayesian inference probability adjoint inverse method effectively narrowed the sampling range, enhancing both the accuracy and robustness of the proposed approach. Finally, the impact of the covariance matrix on the inverse identification accuracy was explored. Overall, this research has provided insights into the future applicability of this Bayesian inference inversion technique for point source identification.

Exoplanet characterization using conditional invertible neural networks

Context. The characterization of the interior of an exoplanet is an inverse problem. The solution requires statistical methods such as Bayesian inference. Current methods employ Markov chain Monte Carlo (MCMC) sampling to infer the posterior probability of the planetary structure parameters for a given exoplanet. These methods are time-consuming because they require the evaluation of a planetary structure model ~105 times. Aims. To speed up the inference process when characterizing an exoplanet, we propose to use conditional invertible neural networks to calculate the posterior probability of the planetary structure parameters. Methods. Conditional invertible neural networks (cINNs) are a special type of neural network that excels at solving inverse problems. We constructed a cINN following the framework for easily invertible architectures (FreIA). This neural network was then trained on a database of 5.6 × 106 internal structure models to recover the inverse mapping between internal structure parameters and observable features (i.e., planetary mass, planetary radius, and elemental composition of the host star). We also show how observational uncertainties can be accounted for. Results. The cINN method was compared to a commonly used Metropolis-Hastings MCMC. To do this, we repeated the characterization of the exoplanet K2-111 b, using both the MCMC method and the trained cINN. We show that the inferred posterior probability distributions of the internal structure parameters from both methods are very similar; the largest differences are seen in the exoplanet water content. Thus, cINNs are a possible alternative to the standard time-consuming sampling methods. cINNs allow infering the composition of an exoplanet that is orders of magnitude faster than what is possible using an MCMC method. The computation of a large database of internal structures to train the neural network is still required, however. Because this database is only computed once, we found that using an invertible neural network is more efficient than an MCMC when more than ten exoplanets are characterized using the same neural network.

Plateau proposal distributions for adaptive component-wise multiple-try metropolis

Markov chain Monte Carlo (MCMC) methods are sampling methods that have become a commonly used tool in statistics, for example to perform Monte Carlo integration. As a consequence of the increase in computational power, many variations of MCMC methods exist for generating samples from arbitrary, possibly complex, target distributions. The performance of an MCMC method, in particular that of a Metropolis–Hastings MCMC method, is predominately governed by the choice of the so-called proposal distribution used. In this paper, we introduce a new type of proposal distribution for the use in Metropolis–Hastings MCMC methods that operates component-wise and with multiple trials per iteration. Specifically, the novel class of proposal distributions, called Plateau distributions, does not overlap, thus ensuring that the multiple trials are drawn from different regions of the state space. Furthermore, the Plateau proposal distributions allow for a bespoke adaptation procedure that lends itself to a Markov chain with efficient problem dependent state space exploration and favourable burn-in properties. Simulation studies show that our novel MCMC algorithm outperforms competitors when sampling from distributions with a complex shape, highly correlated components or multiple modes.

A time-domain multisource Bayesian/Markov chain Monte Carlo formulation of time-lapse seismic waveform inversion

Seismic waveform inversion has been shown to be well-suited for identification and characterization of time-lapse changes in a reservoir. However, the subtle medium variations associated with enhanced oil recovery and/or [Formula: see text] storage problems give uncertainty quantification within such inverse approaches a heightened importance. To analyze both of these features of the time-lapse inverse problem, we formulate a Bayesian full-waveform inversion (FWI) procedure, based on a Markov chain Monte Carlo (MCMC) algorithm. The formulation uses several existing strategies, such as the use of a double-difference time-lapse FWI (DDFWI), incorporation of time-domain multisource data, and application of a local-updating target-oriented inversion. However, it incorporates these within a stochastic framework, involving computation of model covariance with an adaptive Metropolis algorithm, and a method to estimate data error statistics based on the features of time-lapse difference data that is incorporated. A random walk Metropolis-Hastings MCMC is adopted for optimization. In conventional, i.e., deterministic, DDFWI, inversions are carried out for the baseline and monitoring models; in the MCMC approach, a deterministic FWI procedure is carried out for the baseline model, and the MCMC algorithm is applied in the monitoring inversion stage; and the final time-lapse model is the difference between these. A feasibility study is carried out using synthetic 2D acoustic models and data, including time-lapse model estimation and uncertainty quantification. We compare the MCMC approach with conventional deterministic optimization DDFWI and remark on benefits derived, which appear to justify the expanded complexity and cost of a global approach. In addition to the availability of posterior distributions, which are critical for the assessment of the estimations, we observe that the MCMC approach tends to produce monitoring images with clearer edges and fewer coherent errors.

Modeling dragonfly population data with a Bayesian bivariate geometric mixed-effects model

We develop a generalized linear mixed model (GLMM) for bivariate count responses for statistically analyzing dragonfly population data from the Northern Netherlands. The populations of the threatened dragonfly species Aeshna viridis were counted in the years 2015–2018 at 17 different locations (ponds and ditches). Two different widely applied population size measures were used to quantify the population sizes, namely the number of found exoskeletons (‘exuviae’) and the number of spotted egg-laying females were counted. Since both measures (responses) led to many zero counts but also feature very large counts, our GLMM model builds on a zero-inflated bivariate geometric (ZIBGe) distribution, for which we show that it can be easily parameterized in terms of a correlation parameter and its two marginal medians. We model the medians with linear combinations of fixed (environmental covariates) and random (location-specific intercepts) effects. Modeling the medians yields a decreased sensitivity to overly large counts; in particular, in light of growing marginal zero inflation rates. Because of the relatively small sample size (n = 114) we follow a Bayesian modeling approach and use Metropolis-Hastings Markov Chain Monte Carlo (MCMC) simulations for generating posterior samples.

Maximum—A Posteriori Estimation of Linear Time-Invariant State-Space Models via Efficient Monte-Carlo Sampling

Abstract This article addresses maximum-a-posteriori (MAP) estimation of linear time-invariant state-space (LTI-SS) models. The joint posterior distribution of the model matrices and the unknown state sequence is approximated by using Rao-Blackwellized Monte-Carlo sampling algorithms. Specifically, the conditional distribution of the state sequence given the model parameters is derived analytically, while only the marginal posterior distribution of the model matrices is approximated using a Metropolis-Hastings Markov Chain Monte-Carlo sampler. From the joint distribution, MAP estimates of the unknown model matrices as well as the state sequence are computed. The performance of the proposed algorithm is demonstrated on a numerical example and on a real laboratory benchmark dataset of a hair dryer process.

Uncertainty quantification in reservoirs with faults using a sequential approach

Reservoir simulation is critically important for optimally managing petroleum reservoirs. Often, many of the parameters of the model are unknown and cannot be measured directly. These parameters must then be inferred from production data at the wells. This is an inverse problem which can be formulated within a Bayesian framework to integrate prior knowledge with observational data. Markov Chain Monte Carlo (MCMC) methods are commonly used to solve Bayesian inverse problems by generating a set of samples which can be used to characterize the posterior distribution. In this work, we present a novel MCMC algorithm which uses a sequential transition kernel designed to exploit the redundancy which is often present in time series data from reservoirs. This method can be used to efficiently generate samples from the Bayesian posterior for time-dependent models. While this method is general and could be useful for many different models. We consider a Bayesian inverse problem in which we wish to infer fault transmissibilities from measurements of pressure at wells using a two-phase flow model. We demonstrate how the sequential MCMC algorithm presented here can be more efficient than a standard Metropolis-Hastings MCMC approach for this inverse problem. We use integrated autocorrelation times along with mean-squared jump distances to determine the performance of each method for the inverse problem.

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.

Individual-Level Modelling of Infectious Disease Data: EpiILM

In this article, we introduce the R package EpiILM, which provides tools for simulation from, and inference for, discrete-time individual-level models of infectious disease transmission proposed by Deardon et al. (2010). The inference is set in a Bayesian framework and is carried out via Metropolis-Hastings Markov chain Monte Carlo (MCMC). For its fast implementation, key functions are coded in Fortran. Both spatial and contact network models are implemented in the package and can be set in either susceptible-infected (SI) or susceptible-infected-removed (SIR) compartmental frameworks. The use of the package is demonstrated through examples involving both simulated and real data.

Bayesian Approach for Recovering Piecewise Constant Viscoelasticity from MRE Data

This paper deals with an inverse problem for recovering the piecewise constant viscoelasticity of a living body from MRE (Magnetic Resonance Elastography) data. Based on a scalar partial differential equation whose solution can approximately simulate MRE data, our inverse coefficient problem is considered as a statistical inverse problem of reconstructing the posterior distribution of unknown viscoelastic modulus. For sampling this distribution, one usually can use the Metropolis-Hastings Markov chain Monte Carlo (MH-MCMC) algorithm. However, without an appropriate “proposal” distribution given artificially, the MH-MCMC algorithm is hard to draw samples efficiently. To avoid this, a so-called slice sampling algorithm is introduced in this paper and applied for solving our problem. The performance of these statistical inversion algorithms is numerically tested basing on simulated data.

A two-stage Markov chain Monte Carlo method for seismic inversion and uncertainty quantification

Bayesian methods for full-waveform inversion allow quantification of uncertainty in the solution, including determination of interval estimates and posterior distributions of the model unknowns. Markov chain Monte Carlo (MCMC) methods produce posterior distributions subject to fewer assumptions, such as normality, than deterministic Bayesian methods. However, MCMC is computationally a very expensive process that requires repeated solution of the wave equation for different velocity samples. Ultimately, a large proportion of these samples (often 40%–90%) is rejected. We have evaluated a two-stage MCMC algorithm that uses a coarse-grid filter to quickly reject unacceptable velocity proposals, thereby reducing the computational expense of solving the velocity inversion problem and quantifying uncertainty. Our filter stage uses operator upscaling, which provides near-perfect speedup in parallel with essentially no communication between processes and produces data that are highly correlated with those obtained from the full fine-grid solution. Four numerical experiments demonstrate the efficiency and accuracy of the method. The two-stage MCMC algorithm produce the same results (i.e., posterior distributions and uncertainty information, such as medians and highest posterior density intervals) as the Metropolis-Hastings MCMC. Thus, no information needed for uncertainty quantification is compromised when replacing the one-stage MCMC with the more computationally efficient two-stage MCMC. In four representative experiments, the two-stage method reduces the time spent on rejected models by one-third to one-half, which is important because most of models tried during the course of the MCMC algorithm are rejected. Furthermore, the two-stage MCMC algorithm substantially reduced the overall time-per-trial by as much as 40%, while increasing the acceptance rate from 9% to 90%.

Bayesian Detection of Piecewise Linear Trends in Replicated Time-Series with Application to Growth Data Modelling.

We consider the situation where a temporal process is composed of contiguous segments with differing slopes and replicated noise-corrupted time series measurements are observed. The unknown mean of the data generating process is modelled as a piecewise linear function of time with an unknown number of change-points. We develop a Bayesian approach to infer the joint posterior distribution of the number and position of change-points as well as the unknown mean parameters. A-priori, the proposed model uses an overfitting number of mean parameters but, conditionally on a set of change-points, only a subset of them influences the likelihood. An exponentially decreasing prior distribution on the number of change-points gives rise to a posterior distribution concentrating on sparse representations of the underlying sequence. A Metropolis-Hastings Markov chain Monte Carlo (MCMC) sampler is constructed for approximating the posterior distribution. Our method is benchmarked using simulated data and is applied to uncover differences in the dynamics of fungal growth from imaging time course data collected from different strains. The source code is available on CRAN.

Estimation of local heat transfer coefficient from natural convection experiments using liquid crystal thermography and Bayesian method

In this work, an inverse methodology is developed for estimating the local heat transfer coefficients on a vertical plate embedded with the three discrete heat sources, under steady state natural convection, with the temperatures measured at the adiabatic surface without disturbing the fluid flow, using simple conduction/surrogate model and Bayesian inference. Liquid crystal thermography (LCT), an optical measurement method based on the colour-temperature relationship of thermochromic liquid crystal sheet (TLC) is used to determine the temperature field of the adiabatic surface. Bayesian framework with Metropolis Hastings-Markov chain Monte Carlo (MH-MCMC) sampling method is considered for exploring the posterior distribution to estimate the parameters in terms of point estimates like mean, Maximum a posteriori (MAP) and standard deviation. A parity plot between simulated (using retrieved parameters) and measured TLC temperatures shows good agreement.

Orbits for the Impatient: A Bayesian Rejection-sampling Method for Quickly Fitting the Orbits of Long-period Exoplanets

Abstract We describe a Bayesian rejection-sampling algorithm designed to efficiently compute posterior distributions of orbital elements for data covering short fractions of long-period exoplanet orbits. Our implementation of this method, Orbits for the Impatient (OFTI), converges up to several orders of magnitude faster than two implementations of Markov Chain Monte Carlo (MCMC) in this regime. We illustrate the efficiency of our approach by showing that OFTI calculates accurate posteriors for all existing astrometry of the exoplanet 51 Eri b up to 100 times faster than a Metropolis–Hastings MCMC. We demonstrate the accuracy of OFTI by comparing our results for several orbiting systems with those of various MCMC implementations, finding the output posteriors to be identical within shot noise. We also describe how our algorithm was used to successfully predict the location of 51 Eri b six months in the future based on less than three months of astrometry. Finally, we apply OFTI to 10 long-period exoplanets and brown dwarfs, all but one of which have been monitored over less than 3% of their orbits, producing fits to their orbits from astrometric records in the literature.

Stochastic Newton Sampler: The R Package sns

The R package sns implements the stochastic Newton sampler (SNS), a MetropolisHastings Markov chain Monte Carlo (MCMC) algorithm where the proposal density function is a multivariate Gaussian based on a local, second-order Taylor-series expansion of log-density. The mean of the proposal function is the full Newton step in the NewtonRaphson optimization algorithm. Taking advantage of the local, multivariate geometry captured in log-density Hessian allows SNS to be more efficient than univariate samplers, approaching independent sampling as the density function increasingly resembles a multivariate Gaussian. SNS requires the log-density Hessian to be negative-definite everywhere in order to construct a valid proposal function. This property holds, or can be easily checked, for many GLM-like models. When the initial point is far from density peak, running SNS in non-stochastic mode by taking the Newton step - augmented with line search - allows the MCMC chain to converge to high-density areas faster. For high-dimensional problems, partitioning the state space into lower-dimensional subsets, and applying SNS to the subsets within a Gibbs sampling framework can significantly improve the mixing of SNS chains. In addition to the above strategies for improving convergence and mixing, sns offers utilities for diagnostics and visualization, sample-based calculation of Bayesian predictive posterior distributions, numerical differentiation, and log-density validation.

Combining principal component analysis with parameter line-searches to improve the efficacy of Metropolis–Hastings MCMC

When Markov chain Monte Carlo (MCMC) algorithms are used with complex mechanistic models, convergence times are often severely compromised by poor mixing rates and a lack of computational power. Methods such as adaptive algorithms have been developed to improve mixing, but these algorithms are typically highly sophisticated, both mathematically and computationally. Here we present a nonadaptive MCMC algorithm, which we term line-search MCMC, that can be used for efficient tuning of proposal distributions in a highly parallel computing environment, but that nevertheless requires minimal skill in parallel computing to implement. We apply this algorithm to make inferences about dynamical models of the growth of a pathogen (baculovirus) population inside a host (gypsy moth, Lymantria dispar). The line-search MCMC appeal rests on its ease of implementation, and its potential for efficiency improvements over classical MCMC in a highly parallel setting, which makes it especially useful for ecological models.

A Hamiltonian Monte–Carlo method for Bayesian inference of supermassive black hole binaries

We investigate the use of a Hamiltonian Monte–Carlo to map out the posterior density function for supermassive black hole binaries. While previous Markov Chain Monte–Carlo (MCMC) methods, such as Metropolis–Hastings MCMC, have been successfully employed for a number of different gravitational wave sources, these methods are essentially random walk algorithms. The Hamiltonian Monte–Carlo treats the inverse likelihood surface as a ‘gravitational potential’ and by introducing canonical positions and momenta, dynamically evolves the Markov chain by solving Hamiltonʼs equations of motion. This method is not as widely used as other MCMC algorithms due to the necessity of calculating gradients of the log-likelihood, which for most applications results in a bottleneck that makes the algorithm computationally prohibitive. We circumvent this problem by using accepted initial phase-space trajectory points to analytically fit for each of the individual gradients. Eliminating the waveform generation needed for the numerical derivatives reduces the total number of required templates for a iteration chain from to . The result is in an implementation of the Hamiltonian Monte–Carlo that is faster, and more efficient by a factor of approximately the dimension of the parameter space, than a Hessian MCMC.

Hzar: hybrid zone analysis using an R software package

We present a new software package (HZAR) that provides functions for fitting molecular genetic and morphological data from hybrid zones to classic equilibrium cline models using the Metropolis-Hastings Markov chain Monte Carlo (MCMC) algorithm. The software applies likelihood functions appropriate for different types of data, including diploid and haploid genetic markers and quantitative morphological traits. The modular design allows flexibility in fitting cline models of varying complexity. To facilitate hypothesis testing, an autofit function is included that allows automated model selection from a set of nested cline models. Cline parameter values, such as cline centre and cline width, are estimated and may be compared statistically across clines. The package is written in the R language and is available through the Comprehensive R Archive Network (CRAN; http://cran.r-project.org/). Here, we describe HZAR and demonstrate its use with a sample data set from a well-studied hybrid zone in western Panama between white-collared (Manacus candei) and golden-collared manakins (M. vitellinus). Comparisons of our results with previously published results for this hybrid zone validate the hzar software. We extend analysis of this hybrid zone by fitting additional models to molecular data where appropriate.

A Flexible Bayesian Framework for Modeling Haplotype Association with Disease, Allowing for Dominance Effects of the Underlying Causative Variants

Multilocus analysis of single-nucleotide-polymorphism (SNP) haplotypes may provide evidence of association with disease, even when the individual loci themselves do not. Haplotype-based methods are expected to outperform single-SNP analyses because (i) common genetic variation can be structured into haplotypes within blocks of strong linkage disequilibrium and (ii) the functional properties of a protein are determined by the linear sequence of amino acids corresponding to DNA variation on a haplotype. Here, I propose a flexible Bayesian framework for modeling haplotype association with disease in population-based studies of candidate genes or small candidate regions. I employ a Bayesian partition model to describe the correlation between marker-SNP haplotypes and causal variants at the underlying functional polymorphism(s). Under this model, haplotypes are clustered according to their similarity, in terms of marker-SNP allele matches, which is used as a proxy for recent shared ancestry. Haplotypes within a cluster are then assigned the same probability of carrying a causal variant at the functional polymorphism(s). In this way, I can account for the dominance effect of causal variants, here corresponding to any deviation from a multiplicative contribution to disease risk. The results of a detailed simulation study demonstrate that there is minimal cost associated with modeling these dominance effects, with substantial gains in power over haplotype-based methods that do not incorporate clustering and that assume a multiplicative model of disease risks.