Markov Chain Monte Carlo Integration Research Articles

Black Box Variational Bayesian Model Averaging

For many decades now, Bayesian Model Averaging (BMA) has been a popular framework to systematically account for model uncertainty that arises in situations when multiple competing models are available to describe the same or similar physical process. The implementation of this framework, however, comes with a multitude of practical challenges including posterior approximation via Markov chain Monte Carlo and numerical integration. We present a Variational Bayesian Inference approach to BMA as a viable alternative to the standard solutions which avoids many of the aforementioned pitfalls. The proposed method is “black box” in the sense that it can be readily applied to many models with little to no model-specific derivation. We illustrate the utility of our variational approach on a suite of examples and discuss all the necessary implementation details. Fully documented Python code with all the examples is provided as well.

Dataset on the epidemiology and genetic diversification of dengue virus (DENV) serotypes and genotypes in Mexico.

Dengue virus (DENV) evolution has had a significant impact on disease pathogenesis, virulence, and epidemiology in Mexico. Novel genotypic variation in DENV serotypes and genotypes may influence the magnitude and severity of dengue epidemics, as evidenced by 2009 data from Veracruz State. The data presented herein is related to the publication entitled “Epidemiological Implications of the Genetic Diversification of Dengue Virus (DENV) Serotypes and Genotypes in Mexico” [1]. Raw data and trees provide epidemiological data on DENV prevalence and a comprehensive phylogeny of both representative sequences collected from an NCBI repository, and 28 additional isolates from acute-phase plasma samples diagnosed with dengue fever or severe dengue (Raw sequencing data is hosted in the public repository Mendeley Data (http://dx.doi.org/10.17632/bf2kdhhf6x.2). Phylogenetic trees for each DENV serotype (DENV-1, -2, -3 and -4) were constructed using these sequences by a maximum likelihood methodology as well as a Bayesian Markov chain Monte Carlo (MCMC) integration approach. Phylogenetic trees exhibited: (1) DENV-1, genotype V, (2) the DENV-2 Asian/American and Asian II genotypes, (3) DENV-3, genotype III, and (4) DENV-4, genotype I. This data can be beneficial for future analyses on DENV serotype and genotype structure and the introduction of novel DENV genotype sequences in the Americas, for the further elucidation of dengue etiology.

Statistical solution of inverse problems using a state reduction

PurposeThe application of statistical inversion theory provides a powerful approach for solving estimation problems including the ability for uncertainty quantification (UQ) by means of Markov chain Monte Carlo (MCMC) methods and Monte Carlo integration. This paper aims to analyze the application of a state reduction technique within different MCMC techniques to improve the computational efficiency and the tuning process of these algorithms.Design/methodology/approachA reduced state representation is constructed from a general prior distribution. For sampling the Metropolis Hastings (MH) Algorithm and the Gibbs sampler are used. Efficient proposal generation techniques and techniques for conditional sampling are proposed and evaluated for an exemplary inverse problem.FindingsFor the MH-algorithm, high acceptance rates can be obtained with a simple proposal kernel. For the Gibbs sampler, an efficient technique for conditional sampling was found. The state reduction scheme stabilizes the ill-posed inverse problem, allowing a solution without a dedicated prior distribution. The state reduction is suitable to represent general material distributions.Practical implicationsThe state reduction scheme and the MCMC techniques can be applied in different imaging problems. The stabilizing nature of the state reduction improves the solution of ill-posed problems. The tuning of the MCMC methods is simplified.Originality/valueThe paper presents a method to improve the solution process of inverse problems within the Bayesian framework. The stabilization of the inverse problem due to the state reduction improves the solution. The approach simplifies the tuning of MCMC methods.

Quantum parameter estimation via dispersive measurement in circuit QED

We investigate the quantum parameter estimation in circuit quantum electrodynamics via dispersive measurement. Based on the Metropolis–Hastings algorithm and the Markov chain Monte Carlo (MCMC) integration, a new algorithm is proposed to calculate the Fisher information by the stochastic master equation. The Fisher information is expressed in the form of log-likelihood functions and further approximated by the MCMC integration. Numerical results show that the evolution of the Fisher information can approach the quantum Fisher information in a short time interval. These results demonstrate the effectiveness of the proposed algorithm. Finally, based on the proposed algorithm, we consider the effects of the measurement operator and the measurement efficiency on the Fisher information.

No-U-turn sampling for fast Bayesian inference in ADMB and TMB: Introducing the adnuts and tmbstan R packages.

Statistical inference is a widely-used, powerful tool for learning about natural processes in diverse fields. The statistical software platforms AD Model Builder (ADMB) and Template Model Builder (TMB) are particularly popular in the ecological literature, where they are typically used to perform frequentist inference of complex models. However, both lack capabilities for flexible and efficient Markov chain Monte Carlo (MCMC) integration. Recently, the no-U-turn sampler (NUTS) MCMC algorithm has gained popularity for Bayesian inference through the software Stan because it is efficient for high dimensional, complex hierarchical models. Here, we introduce the R packages adnuts and tmbstan, which provide NUTS sampling in parallel and interactive diagnostics with ShinyStan. The ADMB source code was modified to provide NUTS, while TMB models are linked directly into Stan. We describe the packages, provide case studies demonstrating their use, and contrast performance against Stan. For TMB models, we show how to test the accuracy of the Laplace approximation using NUTS. For complex models, the performance of ADMB and TMB was typically within +/- 50% the speed of Stan. In one TMB case study we found inaccuracies in the Laplace approximation, potentially leading to biased inference. adnuts provides a new method for estimating hierarchical ADMB models which previously were infeasible. TMB users can fit the same model in both frequentist and Bayesian paradigms, including using NUTS to test the validity of the Laplace approximation of the marginal likelihood for arbitrary subsets of parameters. These software developments extend the available statistical methods of the ADMB and TMB user base with no additional effort by the user.

Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10.

The Bayesian Evolutionary Analysis by Sampling Trees (BEAST) software package has become a primary tool for Bayesian phylogenetic and phylodynamic inference from genetic sequence data. BEAST unifies molecular phylogenetic reconstruction with complex discrete and continuous trait evolution, divergence-time dating, and coalescent demographic models in an efficient statistical inference engine using Markov chain Monte Carlo integration. A convenient, cross-platform, graphical user interface allows the flexible construction of complex evolutionary analyses.

Sequential Hierarchical Regression Imputation

Skip patterns, bounds, and diverse measurement scales often exacerbate the problem of item nonresponse in the analysis of survey data. Sequential, or variable-by-variable imputation techniques have been quite successfully applied to overcome such problems. Most of these techniques have so far focused on relatively simple designs, and studies have demonstrated the consistency of these methods with techniques that draw from a joint posterior predictive distribution of missing data. Here we consider a sequential imputation technique based on a family of hierarchical regression models, extending the sequential approach to correlated data, (e.g., clustered data) and assess its performance. Each of the regression models is tailored to the variable being handled. Computational techniques used to approximate the posterior predictive distributions are based on Markov Chain Monte Carlo (MCMC) and numerical integration to overcome the problem of intractability. We present a simulation study assessing the compatibility of this approach with the joint data generation mechanism. In the scenarios studied, the sequential method leads to well-calibrated estimates and often performs better than methods that are currently available to practitioners.

Adaptive MCMC in Bayesian phylogenetics: an application to analyzing partitioned data in BEAST.

Advances in sequencing technology continue to deliver increasingly large molecular sequence datasets that are often heavily partitioned in order to accurately model the underlying evolutionary processes. In phylogenetic analyses, partitioning strategies involve estimating conditionally independent models of molecular evolution for different genes and different positions within those genes, requiring a large number of evolutionary parameters that have to be estimated, leading to an increased computational burden for such analyses. The past two decades have also seen the rise of multi-core processors, both in the central processing unit (CPU) and Graphics processing unit processor markets, enabling massively parallel computations that are not yet fully exploited by many software packages for multipartite analyses. We here propose a Markov chain Monte Carlo (MCMC) approach using an adaptive multivariate transition kernel to estimate in parallel a large number of parameters, split across partitioned data, by exploiting multi-core processing. Across several real-world examples, we demonstrate that our approach enables the estimation of these multipartite parameters more efficiently than standard approaches that typically use a mixture of univariate transition kernels. In one case, when estimating the relative rate parameter of the non-coding partition in a heterochronous dataset, MCMC integration efficiency improves by > 14-fold. Our implementation is part of the BEAST code base, a widely used open source software package to perform Bayesian phylogenetic inference. guy.baele@kuleuven.be. Supplementary data are available at Bioinformatics online.

The Hartle–Hawking wave function in 2D causal set quantum gravity

We define the Hartle–Hawking no-boundary wave function for causal set theory (CST) over the discrete analogs of spacelike hypersurfaces. Using Markov Chain Monte Carlo and numerical integration methods we analyze the wave function in non-perturbative 2D CST. We find that in the low-temperature regime it is dominated by causal sets which have no continuum counterparts but possess physically interesting geometric properties. Not only do they exhibit a rapid spatial expansion with respect to the discrete proper time, but a high degree of spatial homogeneity. The latter is due to the extensive overlap of the causal pasts of the elements in the final discrete hypersurface and corresponds to high graph connectivity. Our results thus suggest new possibilities for the role of quantum gravity in the observable Universe.

The symplectic geometry of closed equilateral random walks in 3-space

A closed equilateral random walk in 3-space is a selection of unit length vectors giving the steps of the walk conditioned on the assumption that the sum of the vectors is zero. The sample space of such walks with $n$ edges is the $(2n-3)$-dimensional Riemannian manifold of equilateral closed polygons in $\mathbb{R}^3$. We study closed random walks using the symplectic geometry of the $(2n-6)$-dimensional quotient of the manifold of polygons by the action of the rotation group $\operatorname {SO}(3)$. The basic objects of study are the moment maps on equilateral random polygon space given by the lengths of any $(n-3)$-tuple of nonintersecting diagonals. The Atiyah-Guillemin-Sternberg theorem shows that the image of such a moment map is a convex polytope in $(n-3)$-dimensional space, while the Duistermaat-Heckman theorem shows that the pushforward measure on this polytope is Lebesgue measure on $\mathbb{R}^{n-3}$. Together, these theorems allow us to define a measure-preserving set of "action-angle" coordinates on the space of closed equilateral polygons. The new coordinate system allows us to make explicit computations of exact expectations for total curvature and for some chord lengths of closed (and confined) equilateral random walks, to give statistical criteria for sampling algorithms on the space of polygons and to prove that the probability that a randomly chosen equilateral hexagon is unknotted is at least $\frac{1}{2}$. We then use our methods to construct a new Markov chain sampling algorithm for equilateral closed polygons, with a simple modification to sample (rooted) confined equilateral closed polygons. We prove rigorously that our algorithm converges geometrically to the standard measure on the space of closed random walks, give a theory of error estimators for Markov chain Monte Carlo integration using our method and analyze the performance of our method. Our methods also apply to open random walks in certain types of confinement, and in general to walks with arbitrary (fixed) edgelengths as well as equilateral walks.

Markov Chain Monte Carlo Effective Range Analysis of Low-Energy Electron Elastic Scattering from Xenon

Modified effective range theory formulated as a Bayesian statistical model through the combination with Markov Chain Monte Carlo integration and fitting tech- niques is used to check the compatibility of different e − −Xe scattering data such as the total cross-sections, the momentum transfer cross-sections, and the differential cross-sections that were determined experimentally in the region of Ramsauer-Townsend minimum. On the basis of this predictive approach, the most probable value of the scattering length, (−6.51 ± 0.05)a0, is proposed. The present analysis suggests that the non-relativistic spinless effective range theory is suitable for the description of angu- lar and energy dependencies of e − −Xe elastic scattering

Detecting Multivariate Gene Interactions in RNA-Seq Data Using Optimal Bayesian Classification.

Differential gene expression testing is an analysis commonly applied to RNA-Seq data. These statistical tests identify genes that are significantly different across phenotypes. We extend this testing paradigm to multivariate gene interactions from a classification perspective with the goal to detect novel gene interactions for the phenotypes of interest. This is achieved through our novel computational framework comprised of a hierarchical statistical model of the RNA-Seq processing pipeline and the corresponding optimal Bayesian classifier. Through Markov Chain Monte Carlo sampling and Monte Carlo integration, we compute quantities where no analytical formulation exists. The performance is then illustrated on an expression dataset from a dietary intervention study where we identify gene pairs that have low classification error yet were not identified as differentially expressed. Additionally, we have released the software package to perform OBC classification on RNA-Seq data under an open source license and is available at http://bit.ly/obc_package.

Bayesian Approach for the Seismic Fragility Estimation of a Containment Shell Based on the Formation of Through-Wall Cracks

AbstractThe primary containment (PC) in a nuclear power plant provides the last barrier against radiation leakage to the environment. Design specifications fail to assess the performance of a seismically-damaged containment against accidental overpressure. The present work focuses on the seismic fragility estimation of a PC considering formation of through-wall cracks, which may lead to radiation leakage. This is significantly challenging, considering the computation involved in the integration over a multidimensional probability space, detailed numerical modeling of the containment including concrete damage, and nonlinear dynamic analyses of the structure. A Bayesian approach using a Markov chain Monte Carlo (MCMC) integration scheme is adopted, along with more traditional reliability algorithms (point estimate, first-order second moment (FOSM), and response surface), to address the issue of computation while estimating fragility considering both record-to-record variability and structural parameter unce...

Species Delimitation using Genome-Wide SNP Data

The multispecies coalescent has provided important progress for evolutionary inferences, including increasing the statistical rigor and objectivity of comparisons among competing species delimitation models. However, Bayesian species delimitation methods typically require brute force integration over gene trees via Markov chain Monte Carlo (MCMC), which introduces a large computation burden and precludes their application to genomic-scale data. Here we combine a recently introduced dynamic programming algorithm for estimating species trees that bypasses MCMC integration over gene trees with sophisticated methods for estimating marginal likelihoods, needed for Bayesian model selection, to provide a rigorous and computationally tractable technique for genome-wide species delimitation. We provide a critical yet simple correction that brings the likelihoods of different species trees, and more importantly their corresponding marginal likelihoods, to the same common denominator, which enables direct and accurate comparisons of competing species delimitation models using Bayes factors. We test this approach, which we call Bayes factor delimitation (*with genomic data; BFD*), using common species delimitation scenarios with computer simulations. Varying the numbers of loci and the number of samples suggest that the approach can distinguish the true model even with few loci and limited samples per species. Misspecification of the prior for population size θ has little impact on support for the true model. We apply the approach to West African forest geckos (Hemidactylus fasciatus complex) using genome-wide SNP data. This new Bayesian method for species delimitation builds on a growing trend for objective species delimitation methods with explicit model assumptions that are easily tested. [Bayes factor; model testing; phylogeography; RADseq; simulation; speciation.].

Quantifying the natural history of breast cancer

Background:Natural history models of breast cancer progression provide an opportunity to evaluate and identify optimal screening scenarios. This paper describes a detailed Markov model characterising breast cancer tumour progression.Methods:Breast cancer is modelled by a 13-state continuous-time Markov model. The model differentiates between indolent and aggressive ductal carcinomas in situ tumours, and aggressive tumours of different sizes. We compared such aggressive cancers, that is, which are non-indolent, to those which are non-growing and regressing. Model input parameters and structure were informed by the 1978–1984 Ostergotland county breast screening randomised controlled trial. Overlaid on the natural history model is the effect of screening on diagnosis. Parameters were estimated using Bayesian methods. Markov chain Monte Carlo integration was used to sample the resulting posterior distribution.Results:The breast cancer incidence rate in the Ostergotland population was 21 (95% CI: 17–25) per 10 000 woman-years. Accounting for length-biased sampling, an estimated 91% (95% CI: 85–97%) of breast cancers were aggressive. Larger tumours, 21–50 mm, had an average sojourn of 6 years (95% CI: 3–16 years), whereas aggressive ductal carcinomas in situ took around half a month (95% CI: 0–1 month) to progress to the invasive ⩽10 mm state.Conclusion:These tumour progression rate estimates may facilitate future work analysing cost-effectiveness and quality-adjusted life years for various screening strategies.

Calibration and high fidelity measurement of a quantum photonic chip

Integrated quantum photonic circuits are becoming increasingly complex. Accurate calibration of device parameters and detailed characterization of the prepared quantum states are critically important for future progress. Here we report on an effective experimental calibration method based on Bayesian updating and Markov chain Monte Carlo integration. We use this calibration technique to characterize a two qubit chip and extract the reflectivities of its directional couplers. An average quantum state tomography fidelity of 93.79 ± 1.05% against the four Bell states is achieved. Furthermore, comparing the measured density matrices against a model using the non-ideal device parameters derived from the calibration we achieve an average fidelity of 97.57 ± 0.96%. This pinpoints non-ideality of chip parameters as a major factor in the decrease of Bell state fidelity. We also perform quantum state tomography for Bell states while continuously varying photon distinguishability and find excellent agreement with theory.

A Bayesian model for classifying all differentially expressed proteins simultaneously in 2D PAGE gels

BackgroundTwo-dimensional polyacrylamide gel electrophoresis (2D PAGE) is commonly used to identify differentially expressed proteins under two or more experimental or observational conditions. Wu et al (2009) developed a univariate probabilistic model which was used to identify differential expression between Case and Control groups, by applying a Likelihood Ratio Test (LRT) to each protein on a 2D PAGE. In contrast to commonly used statistical approaches, this model takes into account the two possible causes of missing values in 2D PAGE: either (1) the non-expression of a protein; or (2) a level of expression that falls below the limit of detection.ResultsWe develop a global Bayesian model which extends the previously described model. Unlike the univariate approach, the model reported here is able treat all differentially expressed proteins simultaneously. Whereas each protein is modelled by the univariate likelihood function previously described, several global distributions are used to model the underlying relationship between the parameters associated with individual proteins. These global distributions are able to combine information from each protein to give more accurate estimates of the true parameters. In our implementation of the procedure, all parameters are recovered by Markov chain Monte Carlo (MCMC) integration. The 95% highest posterior density (HPD) intervals for the marginal posterior distributions are used to determine whether differences in protein expression are due to differences in mean expression intensities, and/or differences in the probabilities of expression.ConclusionsSimulation analyses showed that the global model is able to accurately recover the underlying global distributions, and identify more differentially expressed proteins than the simple application of a LRT. Additionally, simulations also indicate that the probability of incorrectly identifying a protein as differentially expressed (i.e., the False Discovery Rate) is very low. The source code is available at https://github.com/stevenhwu/BIDE-2D.

Classification of newborn EEG maturity with Bayesian averaging over decision trees

EEG experts can assess a newborn’s brain maturity by visual analysis of age-related patterns in sleep EEG. It is highly desirable to make the results of assessment most accurate and reliable. However, the expert analysis is limited in capability to provide the estimate of uncertainty in assessments. Bayesian inference has been shown providing the most accurate estimates of uncertainty by using Markov Chain Monte Carlo (MCMC) integration over the posterior distribution. The use of MCMC enables to approximate the desired distribution by sampling the areas of interests in which the density of distribution is high. In practice, the posterior distribution can be multimodal, and so that the existing MCMC techniques cannot provide the proportional sampling from the areas of interest. The lack of prior information makes MCMC integration more difficult when a model parameter space is large and cannot be explored in detail within a reasonable time. In particular, the lack of information about EEG feature importance can affect the results of Bayesian assessment of EEG maturity. In this paper we explore how the posterior information about EEG feature importance can be used to reduce a negative influence of disproportional sampling on the results of Bayesian assessment. We found that the MCMC integration tends to oversample the areas in which a model parameter space includes one or more features, the importance of which counted in terms of their posterior use is low. Using this finding, we proposed to cure the results of MCMC integration and then described the results of testing the proposed method on a set of sleep EEG recordings.

Efficient direct sampling MCEM algorithm for latent variable models with binary responses

While latent variable models have been successfully applied in many fields and underpin various modeling techniques, their ability to incorporate categorical responses is hindered due to the lack of accurate and efficient estimation methods. Approximation procedures, such as penalized quasi-likelihood, are computationally efficient, but the resulting estimators can be seriously biased for binary responses. Gauss–Hermite quadrature and Markov Chain Monte Carlo (MCMC) integration based methods can yield more accurate estimation, but they are computationally much more intensive. Estimation methods that can achieve both computational efficiency and estimation accuracy are still under development. This paper proposes an efficient direct sampling based Monte Carlo EM algorithm (DSMCEM) for latent variable models with binary responses. Mixed effects and item factor analysis models with binary responses are used to illustrate this algorithm. Results from two simulation studies and a real data example suggest that, as compared with MCMC based EM, DSMCEM can significantly improve computational efficiency as well as produce equally accurate parameter estimates. Other aspects and extensions of the algorithm are discussed.

Bayesian Inference for Contact Networks Given Epidemic Data

Abstract. In this article, we estimate the parameters of a simple random network and a stochastic epidemic on that network using data consisting of recovery times of infected hosts. The SEIR epidemic model we fit has exponentially distributed transmission times with Gamma distributed exposed and infectious periods on a network where every edge exists with the same probability, independent of other edges. We employ a Bayesian framework and Markov chain Monte Carlo (MCMC) integration to make estimates of the joint posterior distribution of the model parameters. We discuss the accuracy of the parameter estimates under various prior assumptions and show that it is possible in many scientifically interesting cases to accurately recover the parameters. We demonstrate our approach by studying a measles outbreak in Hagelloch, Germany, in 1861 consisting of 188 affected individuals. We provide an R package to carry out these analyses, which is available publicly on the Comprehensive R Archive Network.