Markov Chain Monte Carlo Chains Research Articles

An automated convergence diagnostic for phylogenetic MCMC analyses.

Assessing convergence of Markov chain Monte Carlo (MCMC) based analyses is crucial but challenging, especially so in high dimensional and complex spaces such as the space of phylogenetic trees (treespace). In practice, it is assumed that the target distribution is the unique stationary distribution of the MCMC and convergence is achieved when samples appear to be stationary. Here we leverage recent advances in computational geometry of the treespace and introduce a method that combines classical statistical techniques and algorithms with geometric properties of the treespace to automatically evaluate and assess practical convergence of phylogenetic MCMC analyses. Our method monitors convergence across multiple MCMC chains and achieves high accuracy in detecting both practical convergence and convergence issues within treespace. Furthermore, our approach is developed to allow for real-time evaluation during the MCMC algorithm run, eliminating any of the chain post-processing steps that are currently required. Our tool therefore improves reliability and efficiency of MCMC based phylogenetic inference methods and makes analyses easier to reproduce and compare. We demonstrate the efficacy of our diagnostic via a well-calibrated simulation study and provide examples of its performance on real data sets. Although our method performs well in practice, a significant part of the underlying treespace probability theory is still missing, which creates an excellent opportunity for future mathematical research in this area. The open source package for the phylogenetic inference framework BEAST2, called ASM, that implements these methods, making them accessible through a user-friendly GUI, is available from https://github.com/rbouckaert/asm/. The open source Python package, called tetres, that provides an interface for these methods enabling their applications beyond BEAST2 can be accessed at https://github.com/bioDS/tetres/.

Tropical Logistic Regression Model on Space of Phylogenetic Trees

Classification of gene trees is an important task both in the analysis of multi-locus phylogenetic data, and assessment of the convergence of Markov Chain Monte Carlo (MCMC) analyses used in Bayesian phylogenetic tree reconstruction. The logistic regression model is one of the most popular classification models in statistical learning, thanks to its computational speed and interpretability. However, it is not appropriate to directly apply the standard logistic regression model to a set of phylogenetic trees, as the space of phylogenetic trees is non-Euclidean and thus contradicts the standard assumptions on covariates. It is well-known in tropical geometry and phylogenetics that the space of phylogenetic trees is a tropical linear space in terms of the max-plus algebra. Therefore, in this paper, we propose an analogue approach of the logistic regression model in the setting of tropical geometry. Our proposed method outperforms classical logistic regression in terms of Area under the ROC Curve in numerical examples, including with data generated by the multi-species coalescent model. Theoretical properties such as statistical consistency have been proved and generalization error rates have been derived. Finally, our classification algorithm is proposed as an MCMC convergence criterion for Mr Bayes. Unlike the convergence metric used by Mr Bayes which is only dependent on tree topologies, our method is sensitive to branch lengths and therefore provides a more robust metric for convergence. In a test case, it is illustrated that the tropical logistic regression can differentiate between two independently run MCMC chains, even when the standard metric cannot.

Single MCMC chain parallelisation on decision trees

AbstractDecision trees (DT) are highly famous in machine learning and usually acquire state-of-the-art performance. Despite that, well-known variants like CART, ID3, random forest, and boosted trees miss a probabilistic version that encodes prior assumptions about tree structures and shares statistical strength between node parameters. Existing work on Bayesian DT depends on Markov Chain Monte Carlo (MCMC), which can be computationally slow, especially on high dimensional data and expensive proposals. In this study, we propose a method to parallelise a single MCMC DT chain on an average laptop or personal computer that enables us to reduce its run-time through multi-core processing while the results are statistically identical to conventional sequential implementation. We also calculate the theoretical and practical reduction in run time, which can be obtained utilising our method on multi-processor architectures. Experiments showed that we could achieve 18 times faster running time provided that the serial and the parallel implementation are statistically identical.

Photometric redshift uncertainties in weak gravitational lensing shear analysis: models and marginalization

ABSTRACT Recovering credible cosmological parameter constraints in a weak lensing shear analysis requires an accurate model that can be used to marginalize over nuisance parameters describing potential sources of systematic uncertainty, such as the uncertainties on the sample redshift distribution n(z). Due to the challenge of running Markov chain Monte Carlo (MCMC) in the high-dimensional parameter spaces in which the n(z) uncertainties may be parametrized, it is common practice to simplify the n(z) parametrization or combine MCMC chains that each have a fixed n(z) resampled from the n(z) uncertainties. In this work, we propose a statistically principled Bayesian resampling approach for marginalizing over the n(z) uncertainty using multiple MCMC chains. We self-consistently compare the new method to existing ones from the literature in the context of a forecasted cosmic shear analysis for the HSC three-year shape catalogue, and find that these methods recover statistically consistent error bars for the cosmological parameter constraints for predicted HSC three-year analysis, implying that using the most computationally efficient of the approaches is appropriate. However, we find that for data sets with the constraining power of the full HSC survey data set (and, by implication, those upcoming surveys with even tighter constraints), the choice of method for marginalizing over n(z) uncertainty among the several methods from the literature may modify the 1σ uncertainties on Ωm–S8 constraints by ∼4 per cent, and a careful model selection is needed to ensure credible parameter intervals.

Reducing Geometric Uncertainty in Computational Hemodynamics by Deep Learning-Assisted Parallel-Chain MCMC.

Computational hemodynamic modeling has been widely used in cardiovascular research and healthcare. However, the reliability of model predictions is largely dependent on the uncertainties of modeling parameters and boundary conditions, which should be carefully quantified and further reduced with available measurements. In this work, we focus on propagating and reducing the uncertainty of vascular geometries within a Bayesian framework. A novel deep learning (DL)-assisted parallel Markov chain Monte Carlo (MCMC) method is presented to enable efficient Bayesian posterior sampling and geometric uncertainty reduction. A DL model is built to approximate the geometry-to-hemodynamic map, which is trained actively using online data collected from parallel MCMC chains and utilized for early rejection of unlikely proposals to facilitate convergence with less expensive full-order model evaluations. Numerical studies on two-dimensional aortic flows are conducted to demonstrate the effectiveness and merit of the proposed method.

Probabilistic site characterisation using monitored spudcan installation data

Predicting spudcan penetration is essential for reducing risk during the installation of an offshore mobile jack-up platform. However, predictions made prior to installation often deviate from those recorded in practice, mainly because inherent uncertainties in seabed soils are not adequately considered. In this paper, the characteristics of the seabed are quantified using monitored spudcan installation data by developing a Bayesian framework coupled with Metropolis algorithm-based Markov chain Monte Carlo (MCMC) simulation. Applied to sand overlying clay conditions multiple readings from the load–penetration curve are used to derive statistical interpretations of seabed parameters such as sand thickness, sand friction angle and the undrained shear strength profile of clay. Illustrative examples demonstrate the proposed method using data from geotechnical centrifuge tests and a monitored offshore jack-up. The accuracy of the method is also retrospectively examined against a database of 66 centrifuge tests. The results show that the proposed method can provide consistent soil parameters with reduced uncertainty regardless of the influence of the updating sequence, number of MCMC chains, scale factors and the prior distributions applied. This method offers the offshore industry a new methodology to statistically interpret seabed conditions by incorporating monitored data. Potentially it could be combined with spatial interpolation methods to infer ground conditions and more accurately predict future spudcan installations.

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.

Markov Chains for Horizons MARCH. I. Identifying Biases in Fitting Theoretical Models to Event Horizon Telescope Observations

Abstract We introduce a new Markov Chain Monte Carlo (MCMC) algorithm with parallel tempering for fitting theoretical models of horizon-scale images of black holes to the interferometric data from the Event Horizon Telescope (EHT). The algorithm implements forms of the noise distribution in the data that are accurate for all signal-to-noise ratios. In addition to being trivially parallelizable, the algorithm is optimized for high performance, achieving 1 million MCMC chain steps in under 20 s on a single processor. We use synthetic data for the 2017 EHT coverage of M87 that are generated based on analytic as well as General Relativistic Magnetohydrodynamic (GRMHD) model images to explore several potential sources of biases in fitting models to sparse interferometric data. We demonstrate that a very small number of data points that lie near salient features of the interferometric data exert disproportionate influence on the inferred model parameters. We also show that the preferred orientations of the EHT baselines introduce significant biases in the inference of the orientation of the model images. Finally, we discuss strategies that help identify the presence and severity of such biases in realistic applications.

The Analysis of the Aftershock Sequences of the Recent Mainshocks in Alaska

The forecasting of the evolution of natural hazards is an important and critical problem in natural sciences and engineering. Earthquake forecasting is one such example and is a difficult task due to the complexity of the occurrence of earthquakes. Since earthquake forecasting is typically based on the seismic history of a given region, the analysis of the past seismicity plays a critical role in modern statistical seismology. In this respect, the recent three significant mainshocks that occurred in Alaska (the 2002, Mw 7.9 Denali; the 2018, Mw 7.9 Kodiak; and the 2018, Mw 7.1 Anchorage earthquakes) presented an opportunity to analyze these sequences in detail. This included the modelling of the frequency-magnitude statistics of the corresponding aftershock sequences. In addition, the aftershock occurrence rates were modelled using the Omori–Utsu (OU) law and the Epidemic Type Aftershock Sequence (ETAS) model. For each sequence, the calculation of the probability to have the largest expected aftershock during a given forecasting time interval was performed using both the extreme value theory and the Bayesian predictive framework. For the Bayesian approach, the Markov Chain Monte Carlo (MCMC) sampling of the posterior distribution was performed to generate the chains of the model parameters. These MCMC chains were used to simulate the models forward in time to compute the predictive distributions. The calculation of the probabilities to have the largest expected aftershock to be above a certain magnitude after a mainshock using the Bayesian predictive framework fully takes into account the uncertainties of the model parameters. Moreover, in order to investigate the credibility of the obtained forecasts, several statistical tests were conducted to compare the performance of the earthquake rate models based on the OU formula and the ETAS model. The results indicate that the Bayesian approach combined with the ETAS model produced more robust results than the standard approach based on the extreme value distribution and the OU law.

Correction to 'Bayesian inference reveals positive but subtle effects of experimental fishery closures on marine predator demographics'.

Correction to 'Bayesian inference reveals positive but subtle effects of experimental fishery closures on marine predator demographics'.

Surgeon effects on cataract refractive outcomes are minimal compared with patient comorbidity and gender: an analysis of 490 987 cases

AimTo investigate effect of patient age, gender, comorbidities and surgeon on refractive outcomes following cataract surgery.MethodsStudy population: patients on UK national ophthalmic cataract database on cataract operations undertaken between 1...

Subhorizon linear Nash–Greene perturbations with constraints on H(z) and the deceleration parameter q(z)

The Hubble parameter H ( z ) , as a function of redshift, is modified by the presence of a new term originated from the extrinsic curvature in a four embedded space–time. Assuming an asymptotic expansion factor a ∼ 0 or, equivalently z → ∞ , it is obtained a nearly resemblance of the present model ( β -model) with the w CDM model with H ( z ) ∼ H w C D M ( z ) at background level. We test the models using a pack of recent datasets as the “Gold 2018” growth-rate data, the best-fit Planck2018/ Λ CDM parameters on the Cosmic Microwave Background (CMB), the Baryon Acoustic Oscillations (BAO) measurements, the Pantheon Supernovae type Ia and the Hubble parameter data with redshift ranging from 0 . 01 < z < 2 . 3 . Performing the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) to ascertain the statistical viability of the model from Jeffreys’ scale, we apply a joint likelihood analysis to the data with the Markov Chain Monte Carlo (MCMC) method. We find that the present model is in very good agreement with observations with a close statistical equivalence with the Λ CDM and w CDM cosmologies at 1- σ level. We also show that a mild alleviation of the σ 8 discrepancy between the growth amplitude factor and the matter content ( σ 8 - Ω m ) of the observations from CMB and Large Scale Structure (LSS) probes. A comparison of the aforementioned full pack of data in MCMC chains is made with the resulting MCMC from Pantheon SNIa+ H ( z ) in order to analyse the sensitivity of the models and how they respond to a cosmography analysis on the evolution of H ( z ) and the deceleration parameter q ( z ) . As a result, we find that only the β -model can be unaffected to the variations of the previous datasets due to its several minima in the likelihood (degeneracies) of the related parameters with an overall percentage relative difference up to 8%. • Subhorizon linear perturbations in the embedding context. • Statistical equivalence of the present model with the Λ CDM and w CDM cosmologies at 1- σ level. • Cosmography analysis on the evolution of H ( z ) and the deceleration parameter q ( z ) • Mild alleviation of the σ 8 discrepancy.

Accelerated Bayesian Inversion of Transient Electromagnetic Data Using MCMC Subposteriors

Transient electromagnetic method (TEM) is one of the major tools to image the subsurface resistivity. The gradient-based inversion of TEM data only provides a unique solution using a subjectively defined regularization penalty, leaving the uncertainty of the solution unaddressed. The Bayesian method can be used to estimate the model parameters, as well as quantify their uncertainty. However, it requires far higher computational costs than gradient-based inversion, which limits the Bayesian inversion of TEM data to 1-D assumptions. We propose an accelerated Bayesian method based on Markov chain Monte Carlo (MCMC) subposteriors to perform full 2-D inversion of TEM data. A robust scheme is designed to divide the model space of a TEM profile into subspaces so that independent MCMC chains can be used to update the parameters in each subspace in parallel. The division is based on the coverage of the source loops using a cumulative sensitivity matrix. Then, the subposteriors obtained at each subspace are merged to approximate the full posterior of the model space using a weighting strategy. A numerical test of a 2-D valley model is used to validate the proposed method. The Bayesian inversion successfully obtained posterior that converges to the true model. The median model makes a good inference of the model parameters, while the probability density function gives their uncertainty estimates. The statistical model of interest can be further extracted from the model ensemble. The proposed method provides an effective framework for Bayesian inversion of TEM data with a multidimensional forward operator.

Contaminant transport forecasting in the subsurface using a Bayesian framework

In monitoring subsurface aquifer contamination, we want to predict quantities—fractional flow curves of pollutant concentration—using subsurface fluid flow models with expertise and limited data. A Bayesian approach is considered here and the complexity associated with the simulation study presents an ongoing practical challenge. We use a Karhunen–Loève expansion for the permeability field in conjunction with GPU computing within a two–stage Markov Chain Monte Carlo (MCMC) method. Further reduction in computing costs is addressed by running several MCMC chains. We compare convergence criteria to quantify the uncertainty of predictions. Our contributions are two-fold: we first propose a fitting procedure for the Multivariate Potential Scale Reduction Factor (MPSRF) data that allows us to estimate the number of iterations for convergence. Then we present a careful analysis of ensembles of fractional flow curves suggesting that, for the problem at hand, the number of iterations required for convergence through the MPSRF analysis is excessive. Thus, for practical applications, our results provide an indication that an analysis of the posterior distributions of quantities of interest provides a reliable criterion to terminate MCMC simulations for quantifying uncertainty.

Scalable Bayesian Inference for Coupled Hidden Markov and Semi-Markov Models

Bayesian inference for coupled hidden Markov models frequently relies on data augmentation techniques for imputation of the hidden state processes. Considerable progress has been made on developing such techniques, mainly using Markov chain Monte Carlo (MCMC) methods. However, as the dimensionality and complexity of the hidden processes increase some of these methods become inefficient, either because they produce MCMC chains with high autocorrelation or because they become computationally intractable. Motivated by this fact we developed a novel MCMC algorithm, which is a modification of the forward filtering backward sampling algorithm, that achieves a good balance between computation and mixing properties, and thus can be used to analyze models with large numbers of hidden chains. Even though our approach is developed under the assumption of a Markovian model, we show how this assumption can be relaxed leading to minor modifications in the algorithm. Our approach is particularly well suited to epidemic models, where the hidden Markov chains represent the infection status of an individual through time. The performance of our method is assessed on simulated data on epidemic models for the spread of Escherichia coli O157:H7 in cattle. Supplementary materials for this article are available online.

A transport-based multifidelity preconditioner for Markov chain Monte Carlo

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

Parallel Metropolis–Hastings Coupler

Bayesian methods and their implementations by means of sophisticated Monte Carlo techniques, such as Markov chain Monte Carlo (MCMC) and particle filters, have become very popular in signal processing over the last years. In this letter, we present a novel interacting parallel MCMC scheme, called parallel Metropolis–Hastings coupler (PMHC), where the information provided by different parallel MCMC chains is properly combined by the use of another advanced MCMC method, called normal kernel coupler (NKC). The NKC employs a mixture of densities as proposal density, which is updated according to a population of states. The PMHC is particularly efficient in multimodal scenarios, since it obtains a faster exploration of the state space with respect to other benchmark techniques. Several numerical simulations are provided showing the efficiency and robustness of the proposed method.

Parallel multiple-chain DRAM MCMC for large-scale geosteering inversion and uncertainty quantification

Geosteering is the proactive control of a wellbore placement based on the downhole measurements, aiming at maximizing economic production from the well. The recent development of azimuthal resistivity logging-while-drilling (LWD) tools have a much larger depth of detection, thus sets higher demands for geosteering inversion as more unknown parameters need to be taken into account in the earth model to be inverted. For complicated non-linear problems, traditional deterministic inversion methods are more likely to be trapped near local optima. In general, Markov chain Monte Carlo (MCMC) inversion methods are more capable of finding the global optimal solution and providing additional uncertainty analyses by sampling from the target distribution. However, MCMC methods usually have the problem of slow convergence. Even though optimization methods like delayed rejection (DR) and adaptive Metropolis (AM) were proposed to speed up the convergence, using MCMC methods to solve geosteering inversion problems may still incur unacceptable time cost. To reduce the sampling cost of MCMC methods, we use parallel multiple-chain DRAM MCMC methods to solve geosteering inverse problems and estimate the corresponding uncertainty. A clustering method, density-based spatial clustering of applications with noise (DBSCAN), is applied to select the best solution among many results. The simulation results show that running many relatively short MCMC chains for one problem can obtain a similar result as running a long single chain. Besides, avoiding communications between multiple Markov chains during the sampling process yields almost linear scalability.

Performance test of QU-fitting in cosmic magnetism study

QU-fitting is a standard model-fitting method to reconstruct distribution of magnetic fields and polarized intensity along a line of sight (LOS) from an observed polarization spectrum. In this paper, we examine the performance of QU-fitting by simulating observations of two polarized sources located along the same LOS, varying the widths of the sources and the gap between them in Faraday depth space, systematically. Markov Chain Monte Carlo (MCMC) approach is used to obtain the best-fit parameters for a fitting model, and Akaike and Bayesian Information Criteria (AIC and BIC, respectively) are adopted to select the best model from four fitting models. We find that the combination of MCMC and AIC/BIC works fairly well in model selection and estimation of model parameters in the cases where two sources have relatively small widths and a larger gap in Faraday depth space. On the other hand, when two sources have large width in Faraday depth space, MCMC chain tends to be trapped in a local maximum so that AIC/BIC cannot select a correct model. We discuss the causes and the tendency of the failure of QU-fitting and suggest a way to improve it.

Comprehensive benchmarking of Markov chain Monte Carlo methods for dynamical systems

BackgroundIn quantitative biology, mathematical models are used to describe and analyze biological processes. The parameters of these models are usually unknown and need to be estimated from experimental data using statistical methods. In particular, Markov chain Monte Carlo (MCMC) methods have become increasingly popular as they allow for a rigorous analysis of parameter and prediction uncertainties without the need for assuming parameter identifiability or removing non-identifiable parameters. A broad spectrum of MCMC algorithms have been proposed, including single- and multi-chain approaches. However, selecting and tuning sampling algorithms suited for a given problem remains challenging and a comprehensive comparison of different methods is so far not available.ResultsWe present the results of a thorough benchmarking of state-of-the-art single- and multi-chain sampling methods, including Adaptive Metropolis, Delayed Rejection Adaptive Metropolis, Metropolis adjusted Langevin algorithm, Parallel Tempering and Parallel Hierarchical Sampling. Different initialization and adaptation schemes are considered. To ensure a comprehensive and fair comparison, we consider problems with a range of features such as bifurcations, periodical orbits, multistability of steady-state solutions and chaotic regimes. These problem properties give rise to various posterior distributions including uni- and multi-modal distributions and non-normally distributed mode tails. For an objective comparison, we developed a pipeline for the semi-automatic comparison of sampling results.ConclusionThe comparison of MCMC algorithms, initialization and adaptation schemes revealed that overall multi-chain algorithms perform better than single-chain algorithms. In some cases this performance can be further increased by using a preceding multi-start local optimization scheme. These results can inform the selection of sampling methods and the benchmark collection can serve for the evaluation of new algorithms. Furthermore, our results confirm the need to address exploration quality of MCMC chains before applying the commonly used quality measure of effective sample size to prevent false analysis conclusions.