Markov Chain Monte Carlo Scheme Research Articles

Computational efficiency study of a micro-macro Markov chain Monte Carlo method for molecular dynamics

We study the numerical properties of a recently introduced micro-macro Markov chain Monte Carlo (mM-MCMC) scheme that accelerates sampling of Gibbs distributions when there is a time-scale separation between the complete (microscopic) molecular dynamics and the slow dynamics of a low dimensional reaction coordinate. The micro-macro Markov chain works in three steps: 1) compute the reaction coordinate value associated to the current molecular state; 2) generate a new macroscopic proposal using a proposal distribution; 3) reconstruct a molecular configuration that is consistent with the newly sampled macroscopic value. Here, we systematically study the impact of method parameters on the efficiency of the micro-macro Markov chain Monte Carlo method, namely the chosen macroscopic dynamics and the microscopic reconstruction distribution. We specifically investigate the impact of the macroscopic and reconstruction proposal distributions and their parameters on this efficiency on two prototypical molecular systems, a three-atom molecule and butane. We find, through detailed computational experiments, that the efficiency is largest when these proposal distributions are close to their time-invariant counterparts and that it is robust with respect to parameter changes.

Field-level Lyman-α forest modeling in redshift space via augmented nonlocal Fluctuating Gunn-Peterson Approximation

Context. Devising fast and accurate methods of predicting the Lyman-α forest at the field level, avoiding the computational burden of running large-volume cosmological hydrodynamic simulations, is of fundamental importance to quickly generate the massive set of simulations needed by the state-of-the-art galaxy and Lyα forest spectroscopic surveys. Aims. We present an improved analytical model to predict the Lyα forest at the field level in redshift space from the dark matter field, expanding upon the widely used Fluctuating Gunn-Peterson Approximation (FGPA). Instead of assuming a unique universal relation over the whole considered cosmic volume, we introduce a dependence on the cosmic web environment (knots, filaments, sheets, and voids) in the model, thereby effectively accounting for nonlocal bias. Furthermore, we include a detailed treatment of velocity bias in the redshift space distortion modeling, allowing the velocity bias to be cosmic-web-dependent. Methods. We first mapped the dark matter field from real to redshift space through a particle-based relation including velocity bias, depending on the cosmic web classification of the dark matter field in real space. We then formalized an appropriate functional form for our model, building upon the traditional FGPA and including a cutoff and a boosting factor mimicking a threshold and inverse-threshold bias effect, respectively, with model parameters depending on the cosmic web classification in redshift space. Eventually, we fit the coefficients of the model via an efficient Markov chain Monte Carlo scheme. Results. We find evidence for a significant difference between the same model parameters in different environments, suggesting that for the investigated setup the simple standard FGPA is not able to adequately predict the Lyα forest in the different cosmic web regimes. We reproduce the summary statistics of the reference cosmological hydrodynamic simulation that we use for comparison, yielding an accurate mean transmitted flux, probability distribution function, 3D power spectrum, and bispectrum. In particular, we achieve maximum deviation and average deviation accuracy in the Lyα forest 3D power spectrum of ∼3% and ∼0.1% up to k ∼ 0.4 h Mpc−1, and ∼5% and ∼1.8% up to k ∼ 1.4 h Mpc−1. Conclusions. Our new model outperforms previous analytical efforts to predict the Lyα forest at the field level in all the probed summary statistics, and has the potential to become instrumental in the generation of fast accurate mocks for covariance matrices estimation in the context of current and forthcoming Lyα forest surveys.

Bayesian tomography using polynomial chaos expansion and deep generative networks

SUMMARY Implementations of Markov chain Monte Carlo (MCMC) methods need to confront two fundamental challenges: accurate representation of prior information and efficient evaluation of likelihood functions. The definition and sampling of the prior distribution can often be facilitated by standard dimensionality-reduction techniques such as Principal Component Analysis (PCA). Additionally, PCA-based decompositions can enable the implementation of accurate surrogate models, for instance, based on polynomial chaos expansion (PCE). However, intricate geological priors with sharp contrasts may demand advanced dimensionality-reduction techniques, such as deep generative models (DGMs). Although suitable for prior sampling, these DGMs pose challenges for surrogate modelling. In this contribution, we present a MCMC strategy that combines the high reconstruction performance of a DGM in the form of a variational autoencoder with the accuracy of PCA–PCE surrogate modelling. Additionally, we introduce a physics-informed PCA decomposition to improve accuracy and reduce the computational burden associated with surrogate modelling. Our methodology is exemplified in the context of Bayesian ground-penetrating radar traveltime tomography using channelized subsurface structures, providing accurate reconstructions and significant speed-ups, particularly when the computation of the full-physics forward model is costly.

A Bayesian Semi-parametric Modelling Approach for Area Level Small Area Studies

We present a new semiparametric extension of the Fay-Herriot model, termed the agnostic Fay-Herriot model (AGFH), in which the sampling-level model is expressed in terms of an unknown general function [Formula: see text]. Thus, the AGFH model can express any distribution in the sampling model since the choice of [Formula: see text] is extremely broad. We propose a Bayesian modelling scheme for AGFH where the unknown function [Formula: see text] is assigned a Gaussian Process prior. Using a Metropolis within Gibbs sampling Markov Chain Monte Carlo scheme, we study the performance of the AGFH model, along with that of a hierarchical Bayesian extension of the Fay-Herriot model. Our analysis shows that the AGFH is an excellent modelling alternative when the sampling distribution is non-Normal, especially in the case where the sampling distribution is bounded. It is also the best choice when the sampling variance is high. However, the hierarchical Bayesian framework and the traditional empirical Bayesian framework can be good modelling alternatives when the signal-to-noise ratio is high, and there are computational constraints. AMS subject classification: 62D05; 62F15

Adaptive Bayesian learning for making risk-aware decisions: A case of trauma survival prediction

Decision tree (DT) models provide a transparent approach to prediction of patient’s outcomes within a probabilistic framework. Averaging over DT models under certain conditions can deliver reliable estimates of predictive posterior probability distributions, which is of critical importance in the case of predicting an individual patient’s outcome. Reliable estimations of the distribution can be achieved within the Bayesian framework using Markov chain Monte Carlo (MCMC) and its Reversible Jump extension enabling DT models to grow to a reasonable size. Existing MCMC strategies however have limited ability to control DT structures and tend to sample overgrown DT models, making unreasonably small partitions, thus deteriorating the uncertainty calibration. This happens because the MCMC explores a DT model parameter space within a limited knowledge of the distribution of data partitions. We propose a new adaptive strategy which overcomes this limitation, and show that in the case of predicting trauma outcomes the number of data partitions can be significantly reduced, so that the unnecessary uncertainty of estimating the predictive posterior density is avoided. The proposed and existing strategies are compared in terms of entropy which, being calculated for predicted posterior distributions, represents the uncertainty in decisions. In this framework, the proposed method has outperformed the existing sampling strategies, so that the unnecessary uncertainty in decisions is efficiently avoided.

A Bayesian hierarchical sparse factor model for estimating simultaneous covariance matrices for gestational outcomes in consecutive pregnancies.

Covariance estimation for multiple groups is a key feature for drawing inference from a heterogeneous population. One should seek to share information about common features in the dependence structures across the various groups. In this paper, we introduce a novel approach for estimating the covariance matrices for multiple groups using a hierarchical latent factor model that shrinks the factor loadings across groups toward a global value. Using a sparse spike and slab model on these loading coefficients allows for a sparse formulation of our model. Parameter estimation is accomplished through a Markov chain Monte Carlo scheme, and a model selection approach is used to select the number of factors to use. We validate our model through extensive simulation studies. Finally, we apply our methodology to the NICHD Consecutive Pregnancies Study to estimate the correlations between birth weights and gestational ages of three consecutive birth within four different subgroups (underweight, normal, overweight, and obese) of women.

Kernel-based learning of orthogonal functions

Estimating a set of orthogonal functions from a finite set of noisy data plays a crucial role in several areas such as imaging, dictionary learning and compressed sensing. The problem turns out especially hard due to its intrinsic non-convexity. In this paper, we solve it by recasting it in the framework of multi-task learning in Hilbert spaces, where orthogonality plays a role as inductive bias. Two perspectives are analyzed. The first one is mainly theoretic. It considers a formulation of the problem where non-orthogonal function estimates are seen as noisy data belonging to an infinite-dimensional space from which orthogonal functions have to be reconstructed. We then provide results concerning the existence and the convergence of the optimizers. The second one is more oriented towards applications. It consists in a learning scheme where orthogonal functions are directly inferred from a finite amount of noisy data. It relies on regularization in reproducing kernel Hilbert spaces and on the introduction of special penalty terms promoting orthogonality among tasks. The problem is then cast in a Bayesian framework, overcoming non-convexity through an efficient Markov chain Monte Carlo scheme. If orthogonality is not certain, our scheme can also understand from data if such form of task interaction really holds.

Full Bayesian identification of linear dynamic systems using stable kernels

System identification learns mathematical models of dynamic systems starting from input-output data. Despite its long history, such research area is still extremely active. New challenges are posed by identification of complex physical processes given by the interconnection of dynamic systems. Examples arise in biology and industry, e.g., in the study of brain dynamics or sensor networks. In the last years, regularized kernel-based identification, with inspiration from machine learning, has emerged as an interesting alternative to the classical approach commonly adopted in the literature. In the linear setting, it uses the class of stable kernels to include fundamental features of physical dynamical systems, e.g., smooth exponential decay of impulse responses. Such class includes also unknown continuous parameters, called hyperparameters, which play a similar role as the model discrete order in controlling complexity. In this paper, we develop a linear system identification procedure by casting stable kernels in a full Bayesian framework. Our models incorporate hyperparameters uncertainty and consist of a mixture of dynamic systems over a continuum spectrum of dimensions. They are obtained by overcoming drawbacks related to classical Markov chain Monte Carlo schemes that, when applied to stable kernels, are proved to become nearly reducible (i.e., unable to reconstruct posteriors of interest in reasonable time). Numerical experiments show that full Bayes frequently outperforms the state-of-the-art results on typical benchmark problems. Two real applications related to brain dynamics (neural activity) and sensor networks are also included.

A Bayesian model for multivariate discrete data using spatial and expert information with application to inferring building attributes

When modeling sparsely observed multivariate data, strong prior information elicited from experts can be used to bolster predictive accuracy and counteract sampling bias. Similarly, modeling autocorrelation in space can help make use of co-occurrence patterns present in many types of spatial data. To make use of both expert prior information and spatial structure, we propose a novel graphical model for a spatial Bayesian network developed specifically to address challenges in inferring the attributes of buildings from geographically sparse observational data. This model is implemented as the sum of a spatial multivariate Gaussian random field and a tabular conditional probability function in real-valued space prior to projection onto the probability simplex. This modeling form is especially suitable for the usage of prior information in the form of sets of atomic rules obtained from experts. To perform inference with missing data, we implement a Markov chain Monte Carlo scheme composed of alternating steps of Gibbs sampling of missing entries and Hamiltonian Monte Carlo for model parameters. A case study in building attribution is presented to highlight the advantages and limitations of this approach.

Bayesian clustering of multiple zero-inflated outcomes.

Several applications involving counts present a large proportion of zeros (excess-of-zeros data). A popular model for such data is the hurdle model, which explicitly models the probability of a zero count, while assuming a sampling distribution on the positive integers. We consider data from multiple count processes. In this context, it is of interest to study the patterns of counts and cluster the subjects accordingly. We introduce a novel Bayesian approach to cluster multiple, possibly related, zero-inflated processes. We propose a joint model for zero-inflated counts, specifying a hurdle model for each process with a shifted Negative Binomial sampling distribution. Conditionally on the model parameters, the different processes are assumed independent, leading to a substantial reduction in the number of parameters as compared with traditional multivariate approaches. The subject-specific probabilities of zero-inflation and the parameters of the sampling distribution are flexibly modelled via an enriched finite mixture with random number of components. This induces a two-level clustering of the subjects based on the zero/non-zero patterns (outer clustering) and on the sampling distribution (inner clustering). Posterior inference is performed through tailored Markov chain Monte Carlo schemes. We demonstrate the proposed approach on an application involving the use of the messaging service WhatsApp. This article is part of the theme issue 'Bayesian inference: challenges, perspectives, and prospects'.

Differentiable samplers for deep latent variable models.

Latent variable models are a popular class of models in statistics. Combined with neural networks to improve their expressivity, the resulting deep latent variable models have also found numerous applications in machine learning. A drawback of these models is that their likelihood function is intractable so approximations have to be carried out to perform inference. A standard approach consists of maximizing instead an evidence lower bound (ELBO) obtained based on a variational approximation of the posterior distribution of the latent variables. The standard ELBO can, however, be a very loose bound if the variational family is not rich enough. A generic strategy to tighten such bounds is to rely on an unbiased low-variance Monte Carlo estimate of the evidence. We review here some recent importance sampling, Markov chain Monte Carlo and sequential Monte Carlo strategies that have been proposed to achieve this. This article is part of the theme issue 'Bayesian inference: challenges, perspectives, and prospects'.

On free energy barriers in Gaussian priors and failure of cold start MCMC for high-dimensional unimodal distributions.

We exhibit examples of high-dimensional unimodal posterior distributions arising in nonlinear regression models with Gaussian process priors for which Markov chain Monte Carlo (MCMC) methods can take an exponential run-time to enter the regions where the bulk of the posterior measure concentrates. Our results apply to worst-case initialized ('cold start') algorithms that are local in the sense that their step sizes cannot be too large on average. The counter-examples hold for general MCMC schemes based on gradient or random walk steps, and the theory is illustrated for Metropolis-Hastings adjusted methods such as preconditioned Crank-Nicolson and Metropolis-adjusted Langevin algorithm. This article is part of the theme issue 'Bayesian inference: challenges, perspectives, and prospects'.

Degradation modelling and lifetime assessment for boiler waterwall with incomplete inspection data

Thermal fatigue cracking is a common problem for boiler waterwall tubes in power plants, resulting the risk of unscheduled plant shutdown. Developing an effective degradation model to predict the future condition of waterwall tubes is vital to assess boiler remaining life. Even though many numerical studies are available in the literature for understanding the physical mechanism (thermal fatigue cracking), the development of predictive models for future risks that can inform decision making are limited. One of the key impediments for developing predictive models is the lack of data. In the case of boiler critical subsystems, very few failures (if any) are recorded in its history. In addition, the number of inspections is limited due to access difficulties. Thus, in many practical cases, the condition data available for degradation modelling is sparse with both infrequent and (spatially) incomplete inspections. This paper presents the development of a stochastic degradation model prediction of thermal fatigue cracking severity in a boiler waterwall when the data is limited. The time evolution of the condition indicator for each waterwall tube is modelled using Markov Models and Bayesian identification algorithms are used to tractably address the large amount of missing data. The methodology presented in this paper employs a Markov Chain Monte Carlo (MCMC) scheme to account for the parameter uncertainty in the presence of missing data. Moreover, this paper also incorporates a novel grouping strategy in which the neighbouring components are grouped to tackle the problem of (spatial) data sparsity. The degradation modelling and prediction tools developed in this paper will support renewal decision making by using available onsite asset data and knowledge. A real-world case study of a boiler waterwall operating in an Australian power industry is also presented. It is found that the waterwall degradation appears to be slowing, but an additional inspection is recommended to confirm the trend.

TOI-1055 b: Neptunian planet characterised with HARPS, TESS, and CHEOPS

Context. TOI-1055 is a Sun-like star known to host a transiting Neptune-sized planet on a 17.5-day orbit (TOI-1055 b). Radial velocity (RV) analyses carried out by two independent groups using nearly the same set of HARPS spectra have provided measurements of planetary masses that differ by ∼2σ. Aims. Our aim in this work is to solve the inconsistency in the published planetary masses by significantly extending the set of HARPS RV measurements and employing a new analysis tool that is able to account and correct for stellar activity. Our further aim was to improve the precision on measurements of the planetary radius by observing two transits of the planet with the CHEOPS space telescope. Methods. We fit a skew normal function to each cross correlation function extracted from the HARPS spectra to obtain RV measurements and hyperparameters to be used for the detrending. We evaluated the correlation changes of the hyperparameters along the RV time series using the breakpoint technique. We performed a joint photometric and RV analysis using a Markov chain Monte Carlo scheme to simultaneously detrend the light curves and the RV time series. Results. We firmly detected the Keplerian signal of TOI-1055 b, deriving a planetary mass of Mb = 20.4−2.5+2.6 M⊕ (∼12%). This value is in agreement with one of the two estimates in the literature, but it is significantly more precise. Thanks to the TESS transit light curves combined with exquisite CHEOPS photometry, we also derived a planetary radius of Rb = 3.490−0.064+0.070 R⊕ (∼1.9%). Our mass and radius measurements imply a mean density of ρb = 2.65−0.35+0.37 g cm−3 (∼14%). We further inferred the planetary structure and found that TOI-1055 b is very likely to host a substantial gas envelope with a mass of 0.41−0.20+0.34 M⊕ and a thickness of 1.05−0.29+0.30 R⊕. Conclusions. Our RV extraction combined with the breakpoint technique has played a key role in the optimal removal of stellar activity from the HARPS time series, enabling us to solve the tension in the planetary mass values published so far for TOI-1055 b.

Bayesian Exploration of Phenomenological EoS of Neutron/Hybrid Stars with Recent Observations

The description of the stellar interior of compact stars remains as a big challenge for the nuclear astrophysics community. The consolidated knowledge is restricted to density regions around the saturation of hadronic matter ρ0=2.8×1014gcm−3, regimes where our nuclear models are successfully applied. As one moves towards higher densities and extreme conditions up to the quark/gluons deconfinement, little can be said about the microphysics of the equation of state (EoS). Here, we employ a Markov Chain Monte Carlo (MCMC) strategy to access the variability at high density regions of polytropic piecewise models for neutron star (NS) EoS or possible hybrid stars, i.e., a NS with a small quark-matter core. With a fixed description of the hadronic matter for low density, below the nuclear saturation density, we explore a variety of models for the high density regimes leading to stellar masses near to 2.5M⊙, in accordance with the observations of massive pulsars. The models are constrained, including the observation of the merger of neutrons stars from VIRGO-LIGO and with the pulsar observed by NICER. In addition, we also discuss the possibility of the use of a Bayesian power regression model with heteroscedastic error. The set of EoS from the Laser Interferometer Gravitational-Wave Observatory (LIGO) was used as input and treated as the data set for the testing case.

Bayesian Structure Learning and Sampling of Bayesian Networks with the R Package BiDAG

The R package BiDAG implements Markov chain Monte Carlo (MCMC) methods for structure learning and sampling of Bayesian networks. The package includes tools to search for a maximum a posteriori (MAP) graph and to sample graphs from the posterior distribution given the data. A new hybrid approach to structure learning enables inference in large graphs. In the first step, we define a reduced search space by means of the PC algorithm or based on prior knowledge. In the second step, an iterative order MCMC scheme proceeds to optimize the restricted search space and estimate the MAP graph. Sampling from the posterior distribution is implemented using either order or partition MCMC. The models and algorithms can handle both discrete and continuous data. The BiDAG package also provides an implementation of MCMC schemes for structure learning and sampling of dynamic Bayesian networks.

Block-Wisely Supervised Network Pruning with Knowledge Distillation and Markov Chain Monte Carlo

Structural network pruning is an effective way to reduce network size for deploying deep networks to resource-constrained devices. Existing methods mainly employ knowledge distillation from the last layer of network to guide pruning of the whole network, and informative features from intermediate layers are not yet fully exploited to improve pruning efficiency and accuracy. In this paper, we propose a block-wisely supervised network pruning (BNP) approach to find the optimal subnet from a baseline network based on knowledge distillation and Markov Chain Monte Carlo. To achieve this, the baseline network is divided into small blocks, and block shrinkage can be independently applied to each block under a same manner. Specifically, block-wise representations of the baseline network are exploited to supervise subnet search by encouraging each block of student network to imitate the behavior of the corresponding baseline block. A score metric measuring block accuracy and efficiency is assigned to each block, and block search is conducted under a Markov Chain Monte Carlo scheme to sample blocks from the posterior. Knowledge distillation enables effective feature representations of the student network, and Markov Chain Monte Carlo provides a sampling scheme to find the optimal solution. Extensive evaluations on multiple network architectures and datasets show BNP outperforms the state of the art. For instance, with 0.16% accuracy improvement on the CIFAR-10 dataset, it yields a more compact subnet of ResNet-110 than other methods by reducing 61.24% FLOPs.

Sparse estimation in linear dynamic networks using the stable spline horseshoe prior

Identification of the so-called dynamic networks is one of the most challenging problems appeared recently in control literature. Such systems consist of large-scale interconnected systems, also called modules. To recover full networks dynamics the two crucial steps are topology detection, where one has to infer from data which connections are active, and modules estimation. Since a small percentage of connections are effective in many real systems, the problem finds also fundamental connections with group-sparse estimation. In particular, in the linear setting modules correspond to unknown impulse responses expected to have null norm but in a small fraction of samples. This paper introduces a new Bayesian approach for linear dynamic networks identification where impulse responses are described through the combination of two particular prior distributions. The first one is a block version of the horseshoe prior, a model possessing important global–local shrinkage features. The second one is the stable spline prior, that encodes information on smooth-exponential decay of the modules. The resulting model is called stable spline horseshoe (SSH) prior. It implements aggressive shrinkage of small impulse responses while larger impulse responses are conveniently subject to stable spline regularization. Inference is performed by a Markov Chain Monte Carlo scheme, tailored to the dynamic context and able to efficiently return the posterior of the modules in sampled form. Numerical studies show that the new approach can accurately reconstruct “line by line” networks dynamics without assuming any knowledge on the topology also when thousands of unknown impulse response coefficients must be inferred from data sets of relatively small size.

EpiLPS: A fast and flexible Bayesian tool for estimation of the time-varying reproduction number.

In infectious disease epidemiology, the instantaneous reproduction number [Formula: see text] is a time-varying parameter defined as the average number of secondary infections generated by an infected individual at time t. It is therefore a crucial epidemiological statistic that assists public health decision makers in the management of an epidemic. We present a new Bayesian tool (EpiLPS) for robust estimation of the time-varying reproduction number. The proposed methodology smooths the epidemic curve and allows to obtain (approximate) point estimates and credible intervals of [Formula: see text] by employing the renewal equation, using Bayesian P-splines coupled with Laplace approximations of the conditional posterior of the spline vector. Two alternative approaches for inference are presented: (1) an approach based on a maximum a posteriori argument for the model hyperparameters, delivering estimates of [Formula: see text] in only a few seconds; and (2) an approach based on a Markov chain Monte Carlo (MCMC) scheme with underlying Langevin dynamics for efficient sampling of the posterior target distribution. Case counts per unit of time are assumed to follow a negative binomial distribution to account for potential overdispersion in the data that would not be captured by a classic Poisson model. Furthermore, after smoothing the epidemic curve, a "plug-in'' estimate of the reproduction number can be obtained from the renewal equation yielding a closed form expression of [Formula: see text] as a function of the spline parameters. The approach is extremely fast and free of arbitrary smoothing assumptions. EpiLPS is applied on data of SARS-CoV-1 in Hong-Kong (2003), influenza A H1N1 (2009) in the USA and on the SARS-CoV-2 pandemic (2020-2021) for Belgium, Portugal, Denmark and France.

Monte Carlo sampling of complex actions in extended state spaces

Path integrals with complex actions are encountered for many physical systems ranging from spin- or mass-imbalanced atomic gases and graphene to quantum chromodynamics at finite density to the nonequilibrium evolution of quantum systems. Many computational approaches have been developed for tackling the sign problem emerging for complex actions. Among these, complex Langevin dynamics has the appeal of general applicability. One of its key challenges is the potential convergence of the dynamics to unphysical fixed points. The statistical sampling process at such a fixed point is not based on the physical action and hence leads to wrong predictions. Moreover, its unphysical nature is hard to detect due to the implicit nature of the process. In the present work we set up a general approach based on a Markov chain Monte Carlo scheme in an extended state space. In this approach we derive an explicit real sampling process for generalized complex Langevin dynamics. Subject to a set of constraints, this sampling process is the physical one. These constraints originate from the detailed-balance equationssatisfied by the Monte Carlo scheme. This allows us to rederive complex Langevin dynamics from a new perspective and establishes a framework for the explicit construction of new sampling schemes for complex actions.