Markov Chain Monte Carlo Output Research Articles

An extended colloid filtration theory for modeling Escherichia coli transport in 3-D fracture networks

Microbial transport in fractured carbonate rock using enhanced solutions is a significant and neglected research topic in the literature. We propose an extended colloid filtration theory (CFT) combined with a particle-tracking following streamlines (PTFS) model for the rapid prediction of breakthrough curves (BTCs) and plumes of pathogens in three-dimensional (3-D) discrete fracture networks (DFNs). We adapted CFT in porous media to pathogen transport in fractures containing Terra Rossa (soil) deposits. As an example of the model capability, a simulation was used to predict the 3-D motion field and Escherichia coli count in groundwater originating from the Forcatella managed aquifer recharge (MAR) Facility (Brindisi, Italy) using a DFN composed of 3,900 fractures. In arid regions, MAR facilities are significant for sustaining basic human needs, such as freshwater supply for drinking and crop production. The Markov chain Monte Carlo (MCMC) technique was applied to E. coli counts in the collected water samples to increase data representativeness. The pathogen transport coefficients were further supported by batch filtration tests carried out in the CNR/IRSA Laboratory (Bari, Italy). The mean E. coli attachment rate coefficient of 0.15 × 10−8 m2 d–1 (sticking efficiency = 1.1 × 10−8 m) resulted in a 2.1 log10 removal in 600 m of reclaimed water filtration. The simulation output visualized the E. coli 3-D pathways in groundwater and the positions of contaminated groundwater spring outflows on Forcatella Beach. The simulation results agreed with the mean MCMC output of E. coli concentrations in bathing water under unperturbed geochemical and environmental flow and transport conditions. However, results indicate that concentrations of pathogenic strains, parasites, and enteric viruses may enter the marine environment of MAR sites during flood periods.

A New Bayesian Person‐Fit Analysis Method Using Pivotal Discrepancy Measures

AbstractA common method of checking person‐fit in Bayesian item response theory (IRT) is the posterior‐predictive (PP) method. In recent years, more powerful approaches have been proposed that are based on resampling methods using the popular statistic. There has also been proposed a new Bayesian model checking method based on pivotal discrepancy measures (PDMs). A PDM T is a discrepancy measure that is a pivotal quantity with a known reference distribution. A posterior sample of T can be generated using standard Markov chain Monte Carlo output, and a p‐value is obtained from probability bounds computed on order statistics of the sample. In this paper, we propose a general procedure to apply this PDM method to person‐fit checking in IRT models. We illustrate this using the and measures. Simulation studies are done comparing these with the PP method and one of the more recent resampling methods. The results show that the PDM method is more powerful than the PP method. Under certain conditions, it is more powerful than the resampling method, while in others, it is less. The PDM method is also applied to a real data set.

Semi-exact control functionals from Sard’s method

Summary A novel control variate technique is proposed for the post-processing of Markov chain Monte Carlo output, based on both Stein’s method and an approach to numerical integration due to Sard. The resulting estimators of posterior expected quantities of interest are proven to be polynomially exact in the Gaussian context, while empirical results suggest that the estimators approximate a Gaussian cubature method near the Bernstein–von Mises limit. The main theoretical result establishes a bias-correction property in settings where the Markov chain does not leave the posterior invariant. Empirical results across a selection of Bayesian inference tasks are presented.

Fast Compression of MCMC Output.

We propose cube thinning, a novel method for compressing the output of an MCMC (Markov chain Monte Carlo) algorithm when control variates are available. It allows resampling of the initial MCMC sample (according to weights derived from control variates), while imposing equality constraints on the averages of these control variates, using the cube method (an approach that originates from survey sampling). The main advantage of cube thinning is that its complexity does not depend on the size of the compressed sample. This compares favourably to previous methods, such as Stein thinning, the complexity of which is quadratic in that quantity.

Scheffé Style Simultaneous Credible Bands for Regression Surfaces with Application to Ache Honey Gathering

We propose two simple, easy-to-implement methods for obtaining simultaneous credible bands in hierarchical models from standard Markov chain Monte Carlo output. The methods generalize Scheffe's (1953) approach to this problem, but in a Bayesian context. A small simulation study is followed by an application of the methods to a seasonal model for Ache honey gathering.

A General Three-Parameter Logistic Model With Time Effect.

Within the framework of item response theory, a new and flexible general three-parameter logistic model with response time (G3PLT) is proposed. The advantage of this model is that it can combine time effect, ability, and item difficulty to influence the correct-response probability. In contrast to the traditional response time models used in educational psychology, the new model incorporates the influence of the time effect on the correct-response probability directly, rather than linking them through a hierarchical method via latent and speed parameters as in van der Linden's model. In addition, the Metropolis–Hastings within Gibbs sampling algorithm is employed to estimate the model parameters. Based on Markov chain Monte Carlo output, two Bayesian model assessment methods are used to assess the goodness of fit between models. Finally, two simulation studies and a real data analysis are performed to further illustrate the advantages of the new model over the traditional three-parameter logistic model.

R package rjmcmc: reversible jump MCMC using post‐processing

SummaryThe rjmcmc package for R implements the post‐processing reversible jump Markov chain Monte Carlo (MCMC) algorithm of Barker & Link. MCMC output from each of the models is used to estimate posterior model probabilities and Bayes factors. Automatic differentiation is used to simplify implementation. The package is demonstrated on two examples.

Modeling Population Structure Under Hierarchical Dirichlet Processes

We propose a Bayesian nonparametric model to infer population admixture, extending the hierarchical Dirichlet process to allow for correlation between loci due to linkage disequilibrium. Given multilocus genotype data from a sample of individuals, the proposed model allows inferring and classifying individuals as unadmixed or admixed, inferring the number of subpopulations ancestral to an admixed population and the population of origin of chromosomal regions. Our model does not assume any specific mutation process, and can be applied to most of the commonly used genetic markers. We present a Markov chain Monte Carlo (MCMC) algorithm to perform posterior inference from the model and we discuss some methods to summarize the MCMC output for the analysis of population admixture. Finally, we demonstrate the performance of the proposed model in a real application, using genetic data from the ectodysplasin-A receptor (EDAR) gene, which is considered to be ancestry-informative due to well-known variations in allele frequency as well as phenotypic effects across ancestry. The structure analysis of this dataset leads to the identification of a rare haplotype in Europeans. We also conduct a simulated experiment and show that our algorithm outperforms parametric methods.

Partition Weighted Approach For Estimating the Marginal Posterior Density With Applications

The computation of marginal posterior density in Bayesian analysis is essential in that it can provide complete information about parameters of interest. Furthermore, the marginal posterior density can be used for computing Bayes factors, posterior model probabilities, and diagnostic measures. The conditional marginal density estimator (CMDE) is theoretically the best for marginal density estimation but requires the closed-form expression of the conditional posterior density, which is often not available in many applications. We develop the partition weighted marginal density estimator (PWMDE) to realize the CMDE. This unbiased estimator requires only a single Markov chain Monte Carlo output from the joint posterior distribution and the known unnormalized posterior density. The theoretical properties and various applications of the PWMDE are examined in detail. The PWMDE method is also extended to the estimation of conditional posterior densities. We carry out simulation studies to investigate the empirical performance of the PWMDE and further demonstrate the desirable features of the proposed method with two real data sets from a study of dissociative identity disorder patients and a prostate cancer study, respectively. Supplementary materials for this article are available online.

Automated Sensitivity Computations for Bayesian Markov Chain Monte Carlo Inference: A New Approach for Prior Robustness and Convergence Analysis

Infinitesimal perturbation analysis is a widely used approach to assess the input sensitivities of stochas- tic dynamic systems in the classical simulation context. In this paper, we introduce an efficient nu- merical approach to undertake Infinitesimal perturbation analysis in the context of dependent Markov Chain Monte Carlo simulations widely used in Bayesian inference. Building on recent developments, we develop a scheme based on automatic differentiation to compute posterior sensitivities based on exact (up to computer floating point error) first-order derivatives of draws (Jacobians) with respect to prior input parameters alongside the estimation algorithm. Assessing such local robustness of posterior inference poses a challenge for existing methods such as finite differencing, symbolic differentiation, and likelihood ratio methods due to the complex stochastic dependence structure and computational intensity of dependent sample draws. The proposed methods allow for a comprehensive and exact local sensitivity analysis of Markov Chain Monte Carlo output with respect to all input parameters, i.e., prior hyper-parameters (prior robustness) and chain starting values (convergence). Our imple- mentation focuses on methods for Gibbs-based algorithms widely used for inference and forecasting in complex high-dimensional settings. We illustrate how the methods can help practitioners assess convergence and prior robustness in an application of Bayesian Vector Autoregression with shrinkage priors for US macroeconomic time series data and forecasting.

Quantifying uncertainty in transdimensional Markov chain Monte Carlo using discrete Markov models

Bayesian analysis often concerns an evaluation of models with different dimensionality as is necessary in, for example, model selection or mixture models. To facilitate this evaluation, transdimensional Markov chain Monte Carlo (MCMC) relies on sampling a discrete indexing variable to estimate the posterior model probabilities. However, little attention has been paid to the precision of these estimates. If only few switches occur between the models in the transdimensional MCMC output, precision may be low and assessment based on the assumption of independent samples misleading. Here, we propose a new method to estimate the precision based on the observed transition matrix of the model-indexing variable. Assuming a first-order Markov model, the method samples from the posterior of the stationary distribution. This allows assessment of the uncertainty in the estimated posterior model probabilities, model ranks, and Bayes factors. Moreover, the method provides an estimate for the effective sample size of the MCMC output. In two model selection examples, we show that the proposed approach provides a good assessment of the uncertainty associated with the estimated posterior model probabilities.

Control variates for stochastic gradient MCMC

It is well known that Markov chain Monte Carlo (MCMC) methods scale poorly with dataset size. A popular class of methods for solving this issue is stochastic gradient MCMC (SGMCMC). These methods use a noisy estimate of the gradient of the log-posterior, which reduces the per iteration computational cost of the algorithm. Despite this, there are a number of results suggesting that stochastic gradient Langevin dynamics (SGLD), probably the most popular of these methods, still has computational cost proportional to the dataset size. We suggest an alternative log-posterior gradient estimate for stochastic gradient MCMC which uses control variates to reduce the variance. We analyse SGLD using this gradient estimate, and show that, under log-concavity assumptions on the target distribution, the computational cost required for a given level of accuracy is independent of the dataset size. Next, we show that a different control-variate technique, known as zero variance control variates, can be applied to SGMCMC algorithms for free. This postprocessing step improves the inference of the algorithm by reducing the variance of the MCMC output. Zero variance control variates rely on the gradient of the log-posterior; we explore how the variance reduction is affected by replacing this with the noisy gradient estimate calculated by SGMCMC.

Comparing hierarchical models via the marginalized deviance information criterion.

Hierarchical models are extensively used in pharmacokinetics and longitudinal studies. When the estimation is performed from a Bayesian approach, model comparison is often based on the deviance information criterion (DIC). In hierarchical models with latent variables, there are several versions of this statistic: the conditional DIC (cDIC) that incorporates the latent variables in the focus of the analysis and the marginalized DIC (mDIC) that integrates them out. Regardless of the asymptotic and coherency difficulties of cDIC, this alternative is usually used in Markov chain Monte Carlo (MCMC) methods for hierarchical models because of practical convenience. The mDIC criterion is more appropriate in most cases but requires integration of the likelihood, which is computationally demanding and not implemented in Bayesian software. Therefore, we consider a method to compute mDIC by generating replicate samples of the latent variables that need to be integrated out. This alternative can be easily conducted from the MCMC output of Bayesian packages and is widely applicable to hierarchical models in general. Additionally, we propose some approximations in order to reduce the computational complexity for large-sample situations. The method is illustrated with simulated data sets and 2 medical studies, evidencing that cDIC may be misleading whilst mDIC appears pertinent.

Post-processing of Markov chain Monte Carlo output in Bayesian latent variable models with application to multidimensional scaling

In Bayesian analysis of multidimensional scaling model with MCMC algorithm, we encounter the indeterminacy of rotation, reflection and translation of the parameter matrix of interest. This type of indeterminacy may be seen in other multivariate latent variable models as well. In this paper, we propose to address this indeterminacy problem with a novel, offline post-processing method that is easily implemented using easy-to-use Markov chain Monte Carlo (MCMC) software. Specifically, we propose a post-processing method based on the generalized extended Procrustes analysis to address this problem. The proposed method is compared with four existing methods to deal with indeterminacy thorough analyses of artificial as well as real datasets. The proposed method achieved at least as good a performance as the best existing method. The benefit of the offline processing approach in the era of easy-to-use MCMC software is discussed.

Modified Cholesky Riemann Manifold Hamiltonian Monte Carlo: exploiting sparsity for fast sampling of high-dimensional targets

Riemann manifold Hamiltonian Monte Carlo (RMHMC) has the potential to produce high-quality Markov chain Monte Carlo output even for very challenging target distributions. To this end, a symmetric positive definite scaling matrix for RMHMC is proposed. The scaling matrix is obtained by applying a modified Cholesky factorization to the potentially indefinite negative Hessian of the target log-density. The methodology is able to exploit the sparsity of the Hessian, stemming from conditional independence modeling assumptions, and thus admit fast implementation of RMHMC even for high-dimensional target distributions. Moreover, the methodology can exploit log-concave conditional target densities, often encountered in Bayesian hierarchical models, for faster sampling and more straightforward tuning. The proposed methodology is compared to alternatives for some challenging targets and is illustrated by applying a state-space model to real data.

Automated Sensitivity Analysis for Bayesian Inference via Markov Chain Monte Carlo: Applications to Gibbs Sampling

Bayesian inference relies heavily on numerical Markov chain Monte carlo (MCMC) methods for the estimation of the typically intractable high-dimensional posterior distributions and requires specific inputs. In this paper we introduce a new general and efficient numerical approach to address important robustness concerns of MCMC analysis with respect to prior input assumptions, a major obstacle to wider acceptance of Bayesian inference, including MCMC algorithm performance (convergence) reflected in the dependence on the chain starting values. The approach builds on recent developments in sensitivity analysis of high-dimensional numerical integrals for classical simulation methods using automatic numerical differentiation methods to compute first order derivatives of algorithmic output with respect to all inputs. We introduce a range of new robustness measures based on Jacobian matrices of MCMC output w.r.t. to the two sets of input parameters, prior parameters and chain starting values, to enable researchers to routinely undertake a comprehensive sensitivity analysis of their MCMC results. The methods are implemented for a range of Gibbs samplers and illustrated using both simulated and real data examples.

A Bayesian heteroscedastic GLM with application to fMRI data with motion spikes

We propose a voxel-wise general linear model with autoregressive noise and heteroscedastic noise innovations (GLMH) for analyzing functional magnetic resonance imaging (fMRI) data. The model is analyzed from a Bayesian perspective and has the benefit of automatically down-weighting time points close to motion spikes in a data-driven manner. We develop a highly efficient Markov Chain Monte Carlo (MCMC) algorithm that allows for Bayesian variable selection among the regressors to model both the mean (i.e., the design matrix) and variance. This makes it possible to include a broad range of explanatory variables in both the mean and variance (e.g., time trends, activation stimuli, head motion parameters and their temporal derivatives), and to compute the posterior probability of inclusion from the MCMC output. Variable selection is also applied to the lags in the autoregressive noise process, making it possible to infer the lag order from the data simultaneously with all other model parameters. We use both simulated data and real fMRI data from OpenfMRI to illustrate the importance of proper modeling of heteroscedasticity in fMRI data analysis. Our results show that the GLMH tends to detect more brain activity, compared to its homoscedastic counterpart, by allowing the variance to change over time depending on the degree of head motion.

The Number of MCMC Draws Needed to Compute Bayesian Credible Bounds

ABSTRACTIn the past 20 years, there has been a staggering increase in the use of Bayesian statistical inference, based on Markov chain Monte Carlo (MCMC) methods, to estimate model parameters and other quantities of interest. This trend exists in virtually all areas of engineering and science. In a typical application, researchers will report estimates of parametric functions (e.g., quantiles, probabilities, or predictions of future outcomes) and corresponding intervals from MCMC methods. One difficulty with the use of inferential methods based on Monte Carlo (MC) is that reported results may be inaccurate due to MC error. MC error, however, can be made arbitrarily small by increasing the number of MC draws. Most users of MCMC methods seem to use indirect diagnostics, trial-and-error, or guess-work to decide how long to run a MCMC algorithm and accuracy of MCMC output results is rarely reported. Unless careful analysis is done, reported numerical results may contain digits that are completely meaningless. In this article, we describe an algorithm to provide direct guidance on the number of MCMC draws needed to achieve a desired amount of precision (i.e., a specified number of accurate significant digits) for Bayesian credible interval endpoints.

Markov chain Monte Carlo methods: an introductory example

When the Guide to the Expression of Uncertainty in Measurement (GUM) and methods from its supplements are not applicable, the Bayesian approach may be a valid and welcome alternative. Evaluating the posterior distribution, estimates or uncertainties involved in Bayesian inferences often requires numerical methods to avoid high-dimensional integrations. Markov chain Monte Carlo (MCMC) sampling is such a method—powerful, flexible and widely applied. Here, a concise introduction is given, illustrated by a simple, typical example from metrology.The Metropolis–Hastings algorithm is the most basic and yet flexible MCMC method. Its underlying concepts are explained and the algorithm is given step by step. The few lines of software code required for its implementation invite interested readers to get started. Diagnostics to evaluate the performance and common algorithmic choices are illustrated to calibrate the Metropolis–Hastings algorithm for efficiency. Routine application of MCMC algorithms may be hindered currently by the difficulty to assess the convergence of MCMC output and thus to assure the validity of results. An example points to the importance of convergence and initiates discussion about advantages as well as areas of research. Available software tools are mentioned throughout.

Bayesian model selection using the median probability model

In the Bayesian approach to model selection, models and model specific parameters are treated as unknown quantities and uncertainty about them are expressed through prior distributions. Given the observed data, updating of the prior distribution to the posterior distribution occurs via Bayes' theorem. The posterior probability of a given model may be interpreted as the support it gets based on the observed data. The highest probability model (HPM) that receives the maximum support from the data is a possible choice for model selection. For large model spaces, Markov chain Monte Carlo (MCMC) algorithms are commonly used to estimate the posterior distribution over models. However, estimates of posterior probabilities of individual models based on MCMC output are not reliable because the number of MCMC samples is typically far smaller than the size of the model space. Thus, the HPM is difficult to estimate and for large model spaces it often has a very small posterior probability. An alternative to the HPM is the median probability model (MPM) of Barbieri and Berger, which has been shown to be the optimal model for prediction using a squared error loss function, under certain conditions. In this article we review some of the conditions for which the MPM is optimal, and provide real data examples to evaluate the performance of the MPM under small and large model spaces. We also discuss the behavior of the MPM under collinearity. WIREs Comput Stat 2015, 7:185–193. doi: 10.1002/wics.1352This article is categorized under: Statistical and Graphical Methods of Data Analysis > Bayesian Methods and Theory Statistical and Graphical Methods of Data Analysis > Markov Chain Monte Carlo (MCMC)