Markov Chain Monte Carlo Convergence Research Articles

Reparameterization of extreme value framework for improved Bayesian workflow

Using Bayesian methods for extreme value analysis offers an alternative to frequentist ones, with several advantages such as easily dealing with parametric uncertainty or studying irregular models. However, computations can be challenging and the efficiency of algorithms can be altered by poor parametrization choices. The focus is on the Poisson process characterization of univariate extremes and outline two key benefits of an orthogonal parameterization. First, Markov chain Monte Carlo convergence is improved when applied on orthogonal parameters. This analysis relies on convergence diagnostics computed on several simulations. Second, orthogonalization also helps deriving Jeffreys and penalized complexity priors, and establishing posterior propriety thereof. The proposed framework is applied to return level estimation of Garonne flow data (France).

Improving MCMC Convergence for Joint Astrometry and Radial Velocity Orbit-fits Through Reparameterization

The exoplanet orbit-fitting software package orbitize! was initially designed to fit the orbits of directly imaged planets with relative astrometric measurements using a Markov Chain Monte Carlo (MCMC) algorithm. Since the publication of orbitize! v1.0, the ability to jointly fit radial velocities and astrometry has been incorporated. We first implemented a Basis class into orbitize! that enables users to add and fit in various orbit parameterizations. We then introduced a radial velocity-focused parameterization of the Keplerian orbital elements for the joint radial velocity and astrometry fits. We compared MCMC convergence speeds of the new radial velocity-focused basis to the original orbitize! standard basis for the system HD 190771, which has full orbital coverage in radial velocity data. We found a 16% faster convergence in time with the radial velocity-focused basis. We encourage users to consider using this basis when doing joint radial velocity and astrometry fits.

The G-Wishart Weighted Proposal Algorithm: Efficient Posterior Computation for Gaussian Graphical Models

Gaussian graphical models can capture complex dependency structures among variables. For such models, Bayesian inference is attractive as it provides principled ways to incorporate prior information and to quantify uncertainty through the posterior distribution. However, posterior computation under the conjugate G-Wishart prior distribution on the precision matrix is expensive for general nondecomposable graphs. We therefore propose a new Markov chain Monte Carlo (MCMC) method named the G-Wishart weighted proposal algorithm (WWA). WWA’s distinctive features include delayed acceptance MCMC, Gibbs updates for the precision matrix and an informed proposal distribution on the graph space that enables embarrassingly parallel computations. Compared to existing approaches, WWA reduces the frequency of the relatively expensive sampling from the G-Wishart distribution. This results in faster MCMC convergence, improved MCMC mixing and reduced computing time. Numerical studies on simulated and real data show that WWA provides a more efficient tool for posterior inference than competing state-of-the-art MCMC algorithms. Supplemental materials for the article are available online.

Using Bayesian inference to estimate plausible muscle forces in musculoskeletal models

BackgroundMusculoskeletal modeling is currently a preferred method for estimating the muscle forces that underlie observed movements. However, these estimates are sensitive to a variety of assumptions and uncertainties, which creates difficulty when trying to interpret the muscle forces from musculoskeletal simulations. Here, we describe an approach that uses Bayesian inference to identify plausible ranges of muscle forces for a simple motion while representing uncertainty in the measurement of the motion and the objective function used to solve the muscle redundancy problem.MethodsWe generated a reference elbow flexion–extension motion and computed a set of reference forces that would produce the motion while minimizing muscle excitations cubed via OpenSim Moco. We then used a Markov Chain Monte Carlo (MCMC) algorithm to sample from a posterior probability distribution of muscle excitations that would result in the reference elbow motion. We constructed a prior over the excitation parameters which down-weighted regions of the parameter space with greater muscle excitations. We used muscle excitations to find the corresponding kinematics using OpenSim, where the error in position and velocity trajectories (likelihood function) was combined with the sum of the cubed muscle excitations integrated over time (prior function) to compute the posterior probability density.ResultsWe evaluated the muscle forces that resulted from the set of excitations that were visited in the MCMC chain (seven parallel chains, 500,000 iterations per chain). The estimated muscle forces compared favorably with the reference forces generated with OpenSim Moco, while the elbow angle and velocity from MCMC matched closely with the reference (average RMSE for elbow angle = 2°; and angular velocity = 32°/s). However, our rank plot analyses and potential scale reduction statistics, which we used to evaluate convergence of the algorithm, indicated that the chains did not fully mix.ConclusionsWhile the results from this process are a promising step towards characterizing uncertainty in muscle force estimation, the computational time required to search the solution space with, and the lack of MCMC convergence indicates that further developments in MCMC algorithms are necessary for this process to become feasible for larger-scale models.

Using deep generative neural networks to account for model errors in Markov chain Monte Carlo inversion

SUMMARY Most geophysical inverse problems are non-linear and rely upon numerical forward solvers involving discretization and simplified representations of the underlying physics. As a result, forward modelling errors are inevitable. In practice, such model errors tend to be either completely ignored, which leads to biased and over-confident inversion results, or only partly taken into account using restrictive Gaussian assumptions. Here, we rely on deep generative neural networks to learn problem-specific low-dimensional probabilistic representations of the discrepancy between high-fidelity and low-fidelity forward solvers. These representations are then used to probabilistically invert for the model error jointly with the target geophysical property field, using the computationally cheap, low-fidelity forward solver. To this end, we combine a Markov chain Monte Carlo (MCMC) inversion algorithm with a trained convolutional neural network of the spatial generative adversarial network (SGAN) type, whereby at each MCMC step, the simulated low-fidelity forward response is corrected using a proposed model-error realization. Considering the crosshole ground-penetrating radar traveltime tomography inverse problem, we train SGAN networks on traveltime discrepancy images between: (1) curved-ray (high fidelity) and straight-ray (low fidelity) forward solvers; and (2) finite-difference-time-domain (high fidelity) and straight-ray (low fidelity) forward solvers. We demonstrate that the SGAN is able to learn the spatial statistics of the model error and that suitable representations of both the subsurface model and model error can be recovered by MCMC. In comparison with inversion results obtained when model errors are either ignored or approximated by a Gaussian distribution, we find that our method has lower posterior parameter bias and better explains the observed traveltime data. Our method is most advantageous when high-fidelity forward solvers involve heavy computational costs and the Gaussian assumption of model errors is inappropriate. Unstable MCMC convergence due to non-linearities introduced by our method remain a challenge to be addressed in future work.

R∗: A Robust MCMC Convergence Diagnostic with Uncertainty Using Decision Tree Classifiers

Markov chain Monte Carlo (MCMC) has transformed Bayesian model inference over the past three decades: mainly because of this, Bayesian inference is now a workhorse of applied scientists. Under general conditions, MCMC sampling converges asymptotically to the posterior distribution, but this provides no guarantees about its performance in finite time. The predominant method for monitoring convergence is to run multiple chains and monitor individual chains’ characteristics and compare these to the population as a whole: if within-chain and between-chain summaries are comparable, then this is taken to indicate that the chains have converged to a common stationary distribution. Here, we introduce a new method for diagnosing convergence based on how well a machine learning classifier model can successfully discriminate the individual chains. We call this convergence measure R ∗ . In contrast to the predominant R ˆ , R ∗ is a single statistic across all parameters that indicates lack of mixing, although individual variables’ importance for this metric can also be determined. Additionally, R ∗ is not based on any single characteristic of the sampling distribution; instead it uses all the information in the chain, including that given by the joint sampling distribution, which is currently largely overlooked by existing approaches. We recommend calculating R ∗ using two different machine learning classifiers — gradient-boosted regression trees and random forests — which each work well in models of different dimensions. Because each of these methods outputs a classification probability, as a byproduct, we obtain uncertainty in R ∗ . The method is straightforward to implement and could be a complementary additional check on MCMC convergence for applied analyses.

A convergence diagnostic for Bayesian clustering

AbstractIn many applications of Bayesian clustering, posterior sampling on the discrete state space of cluster allocations is achieved via Markov chain Monte Carlo (MCMC) techniques. As it is typically challenging to design transition kernels to explore this state space efficiently, MCMC convergence diagnostics for clustering applications are especially important. Here we propose a diagnostic tool for discrete‐space MCMC, focusing on Bayesian clustering applications where the model parameters have been integrated out. We construct a Hotelling‐type statistic on the highest probability states, and use regenerative sampling theory to derive its equilibrium distribution. By leveraging information from the unnormalized posterior, our diagnostic offers added protection against seemingly convergent chains in which the relative frequency of visited states is incorrect. The methodology is illustrated with a Bayesian clustering analysis of genetic mutants of the flowering plant Arabidopsis thaliana.This article is categorized under: Statistical Learning and Exploratory Methods of the Data Sciences > Clustering and Classification Statistical Learning and Exploratory Methods of the Data Sciences > Knowledge Discovery Statistical and Graphical Methods of Data Analysis > Markov Chain Monte Carlo

Convergence Diagnostics for Markov Chain Monte Carlo

Markov chain Monte Carlo (MCMC) is one of the most useful approaches to scientific computing because of its flexible construction, ease of use, and generality. Indeed, MCMC is indispensable for performing Bayesian analysis. Two critical questions that MCMC practitioners need to address are where to start and when to stop the simulation. Although a great amount of research has gone into establishing convergence criteria and stopping rules with sound theoretical foundation, in practice, MCMC users often decide convergence by applying empirical diagnostic tools. This review article discusses the most widely used MCMC convergence diagnostic tools. Some recently proposed stopping rules with firm theoretical footing are also presented. The convergence diagnostics and stopping rules are illustrated using three detailed examples.

Non-point source evaluation of groundwater nitrate contamination from agriculture under geologic uncertainty

The long-term effect of non-point source pollution on groundwater from agricultural practices is a major concern globally. Non-point source pollutants such as nitrate that occurs through fertilizers and animal waste eventually make their way into the aquifer by infiltrating soil. The goal of this study is to develop an approach for characterization of nitrate concentrations at potential source locations under conditions of geologic uncertainty. A Bayesian framework using the Markov Chain Monte Carlo (MCMC) approach is developed to estimate posterior probability distributions of non-point sources by incorporating nitrate concentration data as well as geologic uncertainties. The proposed approach is tested using hypothetical contamination scenarios and then validated using an application case study in North Carolina. Uncertainty existing in geologic formation (i.e., heterogeneous hydraulic conductivity field) is treated as prior and used in evaluating the likelihood function that measures the match between observed and simulated concentrations. The likelihood function computation involves a numerical model that simulates nitrate transport in groundwater from non-point agricultural sources and predicts nitrate concentrations at observation wells. Effectiveness of the MCMC approach is evaluated through a convergence analysis. Comparison among different sampling algorithms is carried out with respect to MCMC convergence diagnostics and making inference. The Bayesian inference analysis methodology developed in this research will help decision makers and water managers to identify potential areas for source containment and decide if further sampling is required.

An eco-epidemiological model for Newcastle disease in central zone Tanzania

Newcastle disease is contagious bird disease which affects main domestic and wild avian species. A deterministic compartmental model for Newcastle disease (ND) is developed and analysed using ordinary differential equation theory. The uncertainties of model parameters were therefore examined using Markov Chain Monte Carlo (MCMC) simulations for the data of chicken death cases due to Newcastle disease from five districts in two regions in Tanzania. The parameter distribution was tested using MCMC convergence diagnostics. The graphical diagnostic test for MCMC used include Trace plots or time series plot, two-dimensional parameter plots and autocorrelation function plots. Hence, model parameters were successfully estimated for numerical simulations and the results of simulations were presented.

Bayesian realized-GARCH models for financial tail risk forecasting incorporating the two-sided Weibull distribution

The realized-GARCH framework is extended to incorporate the two-sided Weibull distribution, for the purpose of volatility and tail risk forecasting in a financial time series. Further, the realized range, as a competitor for realized variance or daily returns, is employed as the realized measure in the realized-GARCH framework. Sub-sampling and scaling methods are applied to both the realized range and realized variance, to help deal with inherent micro-structure noise and inefficiency. A Bayesian Markov Chain Monte Carlo (MCMC) method is adapted and employed for estimation and forecasting, while various MCMC efficiency and convergence measures are employed to assess the validity of the method. In addition, the properties of the MCMC estimator are assessed and compared with maximum likelihood, via a simulation study. Compared to a range of well-known parametric GARCH and realized-GARCH models, tail risk forecasting results across seven market indices, as well as two individual assets, clearly favour the proposed realized-GARCH model incorporating the two-sided Weibull distribution; especially those employing the sub-sampled realized variance and sub-sampled realized range.

On The Markov Chain Monte Carlo Convergence Diagnostic of Bayesian Finite Mixture Model for Income Distribution

The assessment and comparison of income inequality and poverty can be supported by estimating the probability distribution of income. Income distributions which are typically heavy-tailed and positively skewed have been estimated both parametric and nonparametric approach. In parametric approach, finite mixtures distributions have been usefully implemented in the modelling of income distributions which has the multimodal characteristic. The Markov Chain Monte Carlo (MCMC) approach is one of the estimation methods which has a good performance in estimating the parameter of Bayesian finite mixture model. The convergence of the MCMC sampler to the posterior distribution is typically assessed using standard diagnostics methods, i.e., Gelman-Rubin method, Geweke method, Raftery-Lewis method and Heidelberger-Welch method. Those methods can give different results to conclude MCMC convergence condition. In this paper, a real sample income data from the Indonesian Family Life Survey (IFLS) 2015 and BidikMisi 2015 are employed to demonstrate the performance of diagnostics tools that assess convergence of the MCMC algorithm in estimating the parameter of Bayesian finite mixture models.

RadVel: The Radial Velocity Modeling Toolkit

RadVel is an open-source Python package for modeling Keplerian orbits in radial velocity (RV) timeseries. RadVel provides a convenient framework to fit RVs using maximum a posteriori optimization and to compute robust confidence intervals by sampling the posterior probability density via Markov Chain Monte Carlo (MCMC). RadVel allows users to float or fix parameters, impose priors, and perform Bayesian model comparison. We have implemented real-time MCMC convergence tests to ensure adequate sampling of the posterior. RadVel can output a number of publication-quality plots and tables. Users may interface with RadVel through a convenient command-line interface or directly from Python. The code is object-oriented and thus naturally extensible. We encourage contributions from the community. Documentation is available at http://radvel.readthedocs.io.

Simultaneous estimation of population receptive field and hemodynamic parameters from single point BOLD responses using Metropolis-Hastings sampling

We introduce a new approach to Bayesian pRF model estimation using Markov Chain Monte Carlo (MCMC) sampling for simultaneous estimation of pRF and hemodynamic parameters. To obtain high performance on commonly accessible hardware we present a novel heuristic consisting of interpolation between precomputed responses for predetermined stimuli and a large cross-section of receptive field parameters. We investigate the validity of the proposed approach with respect to MCMC convergence, tuning and biases. We compare different combinations of pRF - Compressive Spatial Summation (CSS), Dumoulin-Wandell (DW) and hemodynamic (5-parameter and 3-parameter Balloon-Windkessel) models within our framework with and without the usage of the new heuristic. We evaluate estimation consistency and log probability across models. We perform as well a comparison of one model with and without lookup table within the RStan framework using its No-U-Turn Sampler. We present accelerated computation of whole-ROI parameters for one subject. Finally, we discuss risks and limitations associated with the usage of the new heuristic as well as the means of resolving them. We found that the new algorithm is a valid sampling approach to joint pRF/hemodynamic parameter estimation and that it exhibits very high performance.

Automatic Optimal Batch Size Selection for Recursive Estimators of Time-Average Covariance Matrix

ABSTRACTThe time-average covariance matrix (TACM) , where Γk is the auto-covariance function, is an important quantity for the inference of the mean of an -valued stationary process (d ⩾ 1). This article proposes two recursive estimators for Σ with optimal asymptotic mean square error (AMSE) under different strengths of serial dependence. The optimal estimator involves a batch size selection, which requires knowledge of a smoothness parameter , for some β. This article also develops recursive estimators for ϒβ. Combining these two estimators, we obtain a fully automatic procedure for optimal online estimation for Σ. Consistency and convergence rates of the proposed estimators are derived. Applications to confidence region construction and Markov chain Monte Carlo convergence diagnosis are discussed. Supplementary materials for this article are available online.

RWTY (R We There Yet): An R Package for Examining Convergence of Bayesian Phylogenetic Analyses.

Bayesian inference using Markov chain Monte Carlo (MCMC) has become one of the primary methods used to infer phylogenies from sequence data. Assessing convergence is a crucial component of these analyses, as it establishes the reliability of the posterior distribution estimates of the tree topology and model parameters sampled from the MCMC. Numerous tests and visualizations have been developed for this purpose, but many of the most popular methods are implemented in ways that make them inconvenient to use for large data sets. RWTY is an R package that implements established and new methods for diagnosing phylogenetic MCMC convergence in a single convenient interface.

Diagnosing Undersampling Biases in Monte Carlo Eigenvalue and Flux Tally Estimates

This study focuses on understanding the phenomena in Monte Carlo simulations known as undersampling, in which Monte Carlo tally estimates may not encounter a sufficient number of particles during each generation to obtain unbiased tally estimates. Steady-state Monte Carlo simulations were performed using the KENO Monte Carlo tools within the SCALE code system for models of several burnup credit applications with varying degrees of spatial and isotopic complexities, and the incidence and impact of undersampling on eigenvalue and flux estimates were examined. Using an inadequate number of particle histories in each generation was found to produce a maximum bias of ∼100 pcm in eigenvalue estimates and biases that exceeded 10% in fuel pin flux tally estimates. Having quantified the potential magnitude of undersampling biases in eigenvalue and flux tally estimates in these systems, this study then investigated whether Markov Chain Monte Carlo convergence metrics could be integrated into Monte Carlo simulations to predict the onset and magnitude of undersampling biases. Five potential metrics for identifying undersampling biases were implemented in the SCALE code system and evaluated for their ability to predict undersampling biases by comparing the test metric scores with the observed undersampling biases. Of the five convergence metrics that were investigated, three (the Heidelberger-Welch relative half-width, the Gelman-Rubin diagnostic, and tally entropy) showed the potential to accurately predict the behavior of undersampling biases in the responses examined.

Bayesian Inference on Malignant Breast Cancer in Nigeria: A Diagnosis of MCMC Convergence

Background:There has been no previous study to classify malignant breast tumor in details based on Markov Chain Monte Carlo (MCMC) convergence in Western, Nigeria. This study therefore aims to profile patients living with benign and malignant breast tumor in two different hospitals among women of Western Nigeria, with a focus on prognostic factors and MCMC convergence.Materials and Methods:A hospital-based record was used to identify prognostic factors for malignant breast cancer among women of Western Nigeria. This paper describes Bayesian inference and demonstrates its usage to estimation of parameters of the logistic regression via Markov Chain Monte Carlo (MCMC) algorithm. The result of the Bayesian approach is compared with the classical statistics.Results:The mean age of the respondents was 42.2 ±16.6 years with 52% of the women aged between 35-49 years. The results of both techniques suggest that age and women with at least high school education have a significantly higher risk of being diagnosed with malignant breast tumors than benign breast tumors. The results also indicate a reduction of standard errors is associated with the coefficients obtained from the Bayesian approach. In addition, simulation result reveal that women with at least high school are 1.3 times more at risk of having malignant breast lesion in western Nigeria compared to benign breast lesion.Conclusion:We concluded that more efforts are required towards creating awareness and advocacy campaigns on how the prevalence of malignant breast lesions can be reduced, especially among women. The application of Bayesian produces precise estimates for modeling malignant breast cancer.

A comparison of centring parameterisations of Gaussian process-based models for Bayesian computation using MCMC

Markov chain Monte Carlo (MCMC) algorithms for Bayesian computation for Gaussian process-based models under default parameterisations are slow to converge due to the presence of spatial- and other-induced dependence structures. The main focus of this paper is to study the effect of the assumed spatial correlation structure on the convergence properties of the Gibbs sampler under the default non-centred parameterisation and a rival centred parameterisation (CP), for the mean structure of a general multi-process Gaussian spatial model. Our investigation finds answers to many pertinent, but as yet unanswered, questions on the choice between the two. Assuming the covariance parameters to be known, we compare the exact rates of convergence of the two by varying the strength of the spatial correlation, the level of covariance tapering, the scale of the spatially varying covariates, the number of data points, the number and the structure of block updating of the spatial effects and the amount of smoothness assumed in a Matérn covariance function. We also study the effects of introducing differing levels of geometric anisotropy in the spatial model. The case of unknown variance parameters is investigated using well-known MCMC convergence diagnostics. A simulation study and a real-data example on modelling air pollution levels in London are used for illustrations. A generic pattern emerges that the CP is preferable in the presence of more spatial correlation or more information obtained through, for example, additional data points or by increased covariate variability.

A menu-driven software package of Bayesian nonparametric (and parametric) mixed models for regression analysis and density estimation.

Most of applied statistics involves regression analysis of data. In practice, it is important to specify a regression model that has minimal assumptions which are not violated by data, to ensure that statistical inferences from the model are informative and not misleading. This paper presents a stand-alone and menu-driven software package, Bayesian Regression: Nonparametric and Parametric Models, constructed from MATLAB Compiler. Currently, this package gives the user a choice from 83 Bayesian models for data analysis. They include 47 Bayesian nonparametric (BNP) infinite-mixture regression models; 5 BNP infinite-mixture models for density estimation; and 31 normal random effects models (HLMs), including normal linear models. Each of the 78 regression models handles either a continuous, binary, or ordinal dependent variable, and can handle multi-level (grouped) data. All 83 Bayesian models can handle the analysis of weighted observations (e.g., for meta-analysis), and the analysis of left-censored, right-censored, and/or interval-censored data. Each BNP infinite-mixture model has a mixture distribution assigned one of various BNP prior distributions, including priors defined by either the Dirichlet process, Pitman-Yor process (including the normalized stable process), beta (two-parameter) process, normalized inverse-Gaussian process, geometric weights prior, dependent Dirichlet process, or the dependent infinite-probits prior. The software user can mouse-click to select a Bayesian model and perform data analysis via Markov chain Monte Carlo (MCMC) sampling. After the sampling completes, the software automatically opens text output that reports MCMC-based estimates of the model's posterior distribution and model predictive fit to the data. Additional text and/or graphical output can be generated by mouse-clicking other menu options. This includes output of MCMC convergence analyses, and estimates of the model's posterior predictive distribution, for selected functionals and values of covariates. The software is illustrated through the BNP regression analysis of real data.