MCMC Sampling Research Articles

Variance Matrix Priors for Dirichlet Process Mixture Models With Gaussian Kernels

SummaryBayesian mixture modelling is widely used for density estimation and clustering. The Dirichlet process mixture model (DPMM) is the most popular Bayesian non‐parametric mixture modelling approach. In this manuscript, we study the choice of prior for the variance or precision matrix when Gaussian kernels are adopted. Typically, in the relevant literature, the assessment of mixture models is done by considering observations in a space of only a handful of dimensions. Instead, we are concerned with more realistic problems of higher dimensionality, in a space of up to 20 dimensions. We observe that the choice of prior is increasingly important as the dimensionality of the problem increases. After identifying certain undesirable properties of standard priors in problems of higher dimensionality, we review and implement possible alternative priors. The most promising priors are identified, as well as other factors that affect the convergence of MCMC samplers. Our results show that the choice of prior is critical for deriving reliable posterior inferences. This manuscript offers a thorough overview and comparative investigation into possible priors, with detailed guidelines for their implementation. Although our work focuses on the use of the DPMM in clustering, it is also applicable to density estimation.

Computational Methods for Fast Bayesian Model Assessment via Calibrated Posterior p-values

Posterior predictive p-values (ppps) have become popular tools for Bayesian model assessment, being general-purpose and easy to use. However, interpretation can be difficult because their distribution is not uniform under the hypothesis that the model did generate the data. Calibrated ppps (cppps) can be obtained via a bootstrap-like procedure, yet remain unavailable in practice due to high computational cost. This article introduces methods to enable efficient approximation of cppps and their uncertainty for fast model assessment. We first investigate the computational tradeoff between the number of calibration replicates and the number of MCMC samples per replicate. Provided that the MCMC chain from the real data has converged, using short MCMC chains per calibration replicate can save significant computation time compared to naive implementations, without significant loss in accuracy. We propose different variance estimators for the cppp approximation, which can be used to confirm quickly the lack of evidence against model misspecification. As variance estimation uses effective sample sizes of many short MCMC chains, we show these can be approximated well from the real-data MCMC chain. The procedure for cppp is implemented in NIMBLE, a flexible framework for hierarchical modeling that supports many models and discrepancy measures. Supplementary materials for this article are available online.

SB-ETAS: using simulation based inference for scalable, likelihood-free inference for the ETAS model of earthquake occurrences

The rapid growth of earthquake catalogs, driven by machine learning-based phase picking and denser seismic networks, calls for the application of a broader range of models to determine whether the new data enhances earthquake forecasting capabilities. Additionally, this growth demands that existing forecasting models efficiently scale to handle the increased data volume. Approximate inference methods such as inlabru, which is based on the Integrated nested Laplace approximation, offer improved computational efficiencies and the ability to perform inference on more complex point-process models compared to traditional MCMC approaches. We present SB-ETAS: a simulation based inference procedure for the epidemic-type aftershock sequence (ETAS) model. This approximate Bayesian method uses sequential neural posterior estimation (SNPE) to learn posterior distributions from simulations, rather than typical MCMC sampling using the likelihood. On synthetic earthquake catalogs, SB-ETAS provides better coverage of ETAS posterior distributions compared with inlabru. Furthermore, we demonstrate that using a simulation based procedure for inference improves the scalability from O(n2)\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\mathcal {O}(n^2)$$\\end{document} to O(nlogn)\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\mathcal {O}(n\\log n)$$\\end{document}. This makes it feasible to fit to very large earthquake catalogs, such as one for Southern California dating back to 1981. SB-ETAS can find Bayesian estimates of ETAS parameters for this catalog in less than 10 h on a standard laptop, a task that would have taken over 2 weeks using MCMC. Beyond the standard ETAS model, this simulation based framework allows earthquake modellers to define and infer parameters for much more complex models by removing the need to define a likelihood function.

Scalable Bayesian uncertainty quantification with data-driven priors for radio interferometric imaging

Abstract Next-generation radio interferometers like the Square Kilometer Array have the potential to unlock scientific discoveries thanks to their unprecedented angular resolution and sensitivity. One key to unlocking their potential resides in handling the deluge and complexity of incoming data. This challenge requires building radio interferometric (RI) imaging methods that can cope with the massive data sizes and provide high-quality image reconstructions with uncertainty quantification (UQ). This work proposes a method coined QuantifAI to address UQ in RI imaging with data-driven (learned) priors for high-dimensional settings. Our model, rooted in the Bayesian framework, uses a physically motivated model for the likelihood. The model exploits a data-driven convex prior potential, which can encode complex information learned implicitly from simulations and guarantee the log-concavity of the posterior. We leverage probability concentration phenomena of high-dimensional log-concave posteriors to obtain information about the posterior, avoiding MCMC sampling techniques. We rely on convex optimisation methods to compute the MAP estimation, which is known to be faster and better scale with dimension than MCMC strategies. QuantifAI allows us to compute local credible intervals and perform hypothesis testing of structure on the reconstructed image. We propose a novel fast method to compute pixel-wise uncertainties at different scales, which uses 3 and 6 orders of magnitude less likelihood evaluations than other UQ methods like length of the credible intervals and Monte Carlo posterior sampling, respectively. We demonstrate our method by reconstructing RI images in a simulated setting and carrying out fast and scalable UQ, which we validate with MCMC sampling. Our method shows an improved image quality and more meaningful uncertainties than the benchmark method based on a sparsity-promoting prior. QuantifAI’s source code is available from https://github.com/astro-informatics/QuantifAI.

Detecting univariate, bivariate, and overall effects of drug mixtures using Bayesian kernel machine regression

ABSTRACT Background: Innovative analytic approaches to drug studies are needed to understand better the co-use of opioids with non-opioids among people using illicit drugs. One approach is the Bayesian kernel machine regression (BKMR), widely applied in environmental epidemiology to study exposure mixtures but has received far less attention in substance use research. Objective: To describe the utility of the BKMR approach to study the effects of drug substance mixtures on health outcomes. Methods: We simulated data for 200 individuals. Using the Vale and Maurelli method, we simulated multivariate non-normal drug exposure data: xylazine (mean = 300 ng/mL, SD = 100 ng/mL), fentanyl (mean = 200 ng/mL, SD = 71 ng/mL), benzodiazepine (mean = 300 ng/mL, SD = 55 ng/mL), and nitazene (mean = 200 ng/mL, SD = 141 ng/mL) concentrations. We performed 10,000 MCMC sampling iterations with three Markov chains. Model diagnostics included trace plots, r-hat values, and effective sample sizes. We also provided visual relationships of the univariate and bivariate exposure-response and the overall mixture effect. Results: Higher levels of fentanyl and nitazene concentrations were associated with higher levels of the simulated health outcome, controlling for age. Trace plots, r-hat values, and effective sample size statistics demonstrated BKMR stability across multiple Markov chains. Conclusions: Our understanding of drug mixtures tends to be limited to studies of single-drug models. BKMR offers an innovative way to discern which substances pose a greater health risk than other substances and can be applied to assess univariate, bivariate, and cumulative drug effects on health outcomes.

Fast maximum likelihood estimation for general hierarchical models

Hierarchical statistical models are important in applied sciences because they capture complex relationships in data, especially when variables are related by space, time, sampling unit, or other shared features. Existing methods for maximum likelihood estimation that rely on Monte Carlo integration over latent variables, such as Monte Carlo Expectation Maximization (MCEM), suffer from drawbacks in efficiency and/or generality. We harness a connection between sampling-stepping iterations for such methods and stochastic gradient descent methods for non-hierarchical models: many noisier steps can do better than few cleaner steps. We call the resulting methods Hierarchical Model Stochastic Gradient Descent (HMSGD) and show that combining efficient, adaptive step-size algorithms with HMSGD yields efficiency gains. We introduce a one-dimensional sampling-based greedy line search for step-size determination. We implement these methods and conduct numerical experiments for a Gamma-Poisson mixture model, a generalized linear mixed models (GLMMs) with single and crossed random effects, and a multi-species ecological occupancy model with over 3000 latent variables. Our experiments show that the accelerated HMSGD methods provide faster convergence than commonly used methods and are robust to reasonable choices of MCMC sample size.

Ultra-Efficient MCMC for Bayesian Longitudinal Functional Data Analysis

Functional mixed models are widely useful for regression analysis with dependent functional data, including longitudinal functional data with scalar predictors. However, existing algorithms for Bayesian inference with these models only provide either scalable computing or accurate approximations to the posterior distribution, but not both. We introduce a new MCMC sampling strategy for highly efficient and fully Bayesian regression with longitudinal functional data. Using a novel blocking structure paired with an orthogonalized basis reparameterization, our algorithm jointly samples the fixed effects regression functions together with all subject- and replicate-specific random effects functions. Crucially, the joint sampler optimizes sampling efficiency for these key parameters while preserving computational scalability. Perhaps surprisingly, our new MCMC sampling algorithm even surpasses state-of-the-art algorithms for frequentist estimation and variational Bayes approximations for functional mixed models—while also providing accurate posterior uncertainty quantification—and is orders of magnitude faster than existing Gibbs samplers. Simulation studies show improved point estimation and interval coverage in nearly all simulation settings over competing approaches. We apply our method to a large physical activity dataset to study how various demographic and health factors associate with intraday activity. Supplementary materials for this article are available online.

Bayesian Statistical Inference for Factor Analysis Models with Clustered Data

Clustered data are a complex and frequently used type of data. Traditional factor analysis methods are effective for non-clustered data, but they do not adequately capture correlations between multiple observed individuals or variables in clustered data. This paper proposes a Bayesian approach utilizing MCMC and Gibbs sampling algorithms to accurately estimate parameters of interest within the clustered factor analysis model. The mean traversal graph of parameters ensures that the Markov chain converges, and the Bayesian case-deletion model is used to analyze the model’s impact and identify outliers in clustered data using Cook’s posterior mean distance. The applicability and validity of the principal-component-method-based factor analysis model for clustered data are demonstrated by comparing the parameter estimation of this method with the principal component method, the clustered data with and without internal relationships are compared by example analysis, and the anomalous groups are identified by the Cook’s posterior mean distance.

Spatial Prediction of Diameter Distributions for the Alpine Protection Forests in Ebensee, Austria, Using ALS/PLS and Spatial Distributional Regression Models

A novel Bayesian spatial distributional regression model is presented to predict forest structural diversity in terms of the distributions of the stem diameter at breast height (DBH) in the protection forests in Ebensee, Austria. The distributional regression approach overcomes the limitations and uncertainties of traditional regression modeling, in which the conditional mean of the response is regressed against explanatory variables. The distributional regression addresses the complete conditional response distribution, instead. In total 36,338 sample trees were measured via a handheld mobile personal laser scanning system (PLS) on 273 sample plots each having a 20 m radius. Recent airborne laser scanning (ALS) data were used to derive regression covariates from the normalized digital vegetation height model (DVHM) and the digital terrain model (DTM). Candidate models were constructed that differed in their linear predictors of the two gamma distribution parameters. In the distributional regression approach, covariates can enter the model in a flexible form, such as via nonlinear smooth curves, cyclic smooths, or spatial effects. Supported by Bayesian diagnostics DIC and WAIC, nonlinear smoothing splines outperformed linear parametric slope coefficients, and the best implementation of spatial structured effects was achieved by a Gaussian process smooth. Model fitting and posterior parameter inference was achieved by using full Bayesian methodology and MCMC sampling algorithms implemented in the R-package BAMLSS. With BAMLSS, spatial interval predictions of the DBH distribution at any new geo-locations were enabled via straightforward access to the posterior predictive distributions of the model terms and by offering simple plug-in solutions for new covariate values. A cross-validation analysis validated the robustness of the proposed method’s parameter estimation and out-of-sample prediction. Spatial predictions of stem count proportions per DBH classes revealed that regeneration of smaller trees was lacking in certain areas of the protection forest landscape. Therefore, the intensity of final felling needs to be increased to reduce shading from the dense, overmature shelter trees and to promote sunlight for the young regeneration trees.

MCMC sampling of directed flag complexes with fixed undirected graphs

Constructing null models to test the significance of extracted information is a crucial step in data analysis. In this work, we provide a uniformly sampleable null model of directed graphs with the same (or similar) number of simplices in the flag complex, with the restriction of retaining the underlying undirected graph. We describe an MCMC-based algorithm to sample from this null model and statistically investigate the mixing behaviour. This is paired with a high-performance, Rust-based, publicly available implementation. The motivation comes from topological data analysis of connectomes in neuroscience. In particular, we answer the fundamental question: are the high Betti numbers observed in the investigated graphs evidence of an interesting topology, or are they merely a byproduct of the high numbers of simplices? Indeed, by applying our new tool on the connectome of C. elegans and parts of the statistical reconstructions of the Blue Brain Project, we find that the Betti numbers observed are considerable statistical outliers with respect to this new null model. We thus, for the first time, statistically confirm that topological data analysis in microscale connectome research is extracting statistically meaningful information.

Time-domain structural model updating following the Bayesian approach in the absence of system input information

The modal-domain structural model updating and damage detection methods own the advantage that they perform the whole computation without the information of system input, when compared to the time-domain. The system input required in the time-domain methods is hardly measured in realistic field tests of civil engineering structures. However, the time-domain methods own the advantage that they carry out model updating based on raw measured dynamic responses but not the processed (and compressed) modal parameters. The modal-domain methods only uses the compressed modal parameters, and some useful system information may be lost in the compressing process. Thus, higher errors and uncertainties are induced in the model updating results from the modal-domain methods. To take advantages and avoid disadvantages from the traditional time-domain and modal-domain model updating methods, a novel time-domain MCMC-based Bayesian model updating method without the information of the system input is developed in this paper by incorporating a newly developed dynamic response reconstruction technique and extending the formulation of the posterior probability density function (PDF) and the corresponding bridge PDF in the MCMC sampling process. A new modal decomposition algorithm is proposed in this paper for accurate modal response decomposition even in case of dense or sparse modes, and filters are not required for higher computation accuracy. Experimental case studies on a two-story steel frame were conducted to illustrate the proposed methodology in both structural model updating and damage detection. According to the analysis results, it is concluded that the proposed methodology outperforms the modal-domain methods especially in the cases with small number of measured DOFs.

Posterior Representations for Bayesian Context Trees: Sampling, Estimation and Convergence

We revisit the Bayesian Context Trees (BCT) modelling framework for discrete time series, which was recently found to be very effective in numerous tasks including model selection, estimation and prediction. A novel representation of the induced posterior distribution on model space is derived in terms of a simple branching process, and several consequences of this are explored in theory and in practice. First, it is shown that the branching process representation leads to a simple variable-dimensional Monte Carlo sampler for the joint posterior distribution on models and parameters, which can efficiently produce independent samples. This sampler is found to be more efficient than earlier MCMC samplers for the same tasks. Then, the branching process representation is used to establish the asymptotic consistency of the BCT posterior, including the derivation of an almost-sure convergence rate. Finally, an extensive study is carried out on the performance of the induced Bayesian entropy estimator. Its utility is illustrated through both simulation experiments and real-world applications, where it is found to outperform several state-of-the-art methods.

EPSOM-Hyb: A General Purpose Estimator of Log-Marginal Likelihoods with Applications in Probabilistic Graphical Models

We consider the estimation of the marginal likelihood in Bayesian statistics, with primary emphasis on Gaussian graphical models, where the intractability of the marginal likelihood in high dimensions is a frequently researched problem. We propose a general algorithm that can be widely applied to a variety of problem settings and excels particularly when dealing with near log-concave posteriors. Our method builds upon a previously posited algorithm that uses MCMC samples to partition the parameter space and forms piecewise constant approximations over these partition sets as a means of estimating the normalizing constant. In this paper, we refine the aforementioned local approximations by taking advantage of the shape of the target distribution and leveraging an expectation propagation algorithm to approximate Gaussian integrals over rectangular polytopes. Our numerical experiments show the versatility and accuracy of the proposed estimator, even as the parameter space increases in dimension and becomes more complicated.

Being Bayesian about learning Bayesian networks from ordinal data

In this paper we propose a Bayesian approach for inferring Bayesian network (BN) structures from ordinal data. Our approach can be seen as the Bayesian counterpart of a recently proposed frequentist approach, referred to as the ‘ordinal structure expectation maximization’ (OSEM) method. Like for the OSEM method, the key idea is to assume that each ordinal variable originates from a Gaussian variable that can only be observed in discretized form, and that the dependencies in the latent Gaussian space can be modeled by BNs; i.e. by directed acyclic graphs (DAGs). Our Bayesian method combines the ‘structure MCMC sampler’ for DAG posterior sampling, a slightly modified version of the ‘Bayesian metric for Gaussian networks having score equivalence’ (BGe score), the concept of the ‘extended rank likelihood’, and a recently proposed algorithm for posterior sampling the parameters of Gaussian BNs. In simulation studies we compare the new Bayesian approach and the OSEM method in terms of the network reconstruction accuracy. The empirical results show that the new Bayesian approach leads to significantly improved network reconstruction accuracies.

Scores for Learning Discrete Causal Graphs with Unobserved Confounders

Structural learning is arguably one of the most challenging and pervasive tasks found throughout the data sciences. There exists a growing literature that studies structural learning in non-parametric settings where conditional independence constraints are taken to define the equivalence class. In the presence of unobserved confounders, it is understood that non-conditional independence constraints are imposed over the observational distribution, including certain equalities and inequalities between functionals of the joint distribution. In this paper, we develop structural learning methods that leverage additional constraints beyond conditional independences. Specifically, we first introduce a score for arbitrary graphs combining Watanabe's asymptotic expansion of the marginal likelihood and new bounds over the cardinality of the exogenous variables. Second, we show that the new score has desirable properties in terms of expressiveness and computability. In terms of expressiveness, we prove that the score captures distinct constraints imprinted in the data, including Verma's and inequalities'. In terms of computability, we show properties of score equivalence and decomposability, which allows, in principle, to break the problem of structural learning in smaller and more manageable pieces. Third, we implement this score using an MCMC sampling algorithm and test its properties in several simulation scenarios.

Convergence of stratified MCMC sampling of non-reversible dynamics

Convergence of stratified MCMC sampling of non-reversible dynamics

A Dirichlet process mixture regression model for the analysis of competing risk events

We develop a regression model for the analysis of competing risk events. The joint distribution of the time to these events is flexibly characterized by a random effect which follows a discrete probability distribution drawn from a Dirichlet Process, explaining their variability. This entails an additional layer of flexibility of this joint model, whose inference is robust with respect to the misspecification of the distribution of the random effects. The model is analysed in a fully Bayesian setting, yielding a flexible Dirichlet Process Mixture model for the joint distribution of the time to events. An efficient MCMC sampler is developed for inference. The modelling approach is applied to the empirical analysis of the surrending risk in a US life insurance portfolio previously analysed by Milhaud and Dutang (2018). The approach yields an improved predictive performance of the surrending rates.

MCMC algorithm based on Markov random field in image segmentation.

In the realm of digital image applications, image processing technology occupies a pivotal position, with image segmentation serving as a foundational component. As the digital image application domain expands across industries, the conventional segmentation techniques increasingly challenge to cater to modern demands. To address this gap, this paper introduces an MCMC-based image segmentation algorithm based on the Markov Random Field (MRF) model, marking a significant stride in the field. The novelty of this research lies in its method that capitalizes on domain information in pixel space, amplifying the local segmentation precision of image segmentation algorithms. Further innovation is manifested in the development of an adaptive segmentation image denoising algorithm based on MCMC sampling. This algorithm not only elevates image segmentation outcomes, but also proficiently denoises the image. In the experimental results, MRF-MCMC achieves better segmentation performance, with an average segmentation accuracy of 94.26% in Lena images, significantly superior to other common image segmentation algorithms. In addition, the study proposes that the denoising model outperforms other algorithms in peak signal-to-noise ratio and structural similarity in environments with noise standard deviations of 15, 25, and 50. In essence, these experimental findings affirm the efficacy of this study, opening avenues for refining digital image segmentation methodologies.

Screening the Discrepancy Function of a Computer Model

Traditionally, screening refers to the problem of detecting influential (active) inputs in the computer model. We develop methodology that applies to screening, but the main focus is on detecting active inputs not in the computer model itself but rather on the discrepancy function that is introduced to account for model inadequacy when linking the computer model with field observations. We contend this is an important problem as it informs the modeler which are the inputs that are potentially being mishandled in the model, but also along which directions it may be less recommendable to use the model for prediction. The methodology is Bayesian and is inspired by the continuous spike-and-slab prior popularized by the literature on Bayesian variable selection. In our approach, and in contrast with previous proposals, a single MCMC sample from the full model allows us to compute the posterior probabilities of all the competing models, resulting in a methodology that is computationally very fast. The approach hinges on the ability to obtain posterior inclusion probabilities of the inputs, which are easy to interpret quantities, as the basis for selecting active inputs. For that reason, we name the methodology PIPS—posterior inclusion probability screening.

Uncertainty quantification of reservoir performance using machine learning algorithms and structured expert judgment

The increasing demand for fossil energy necessitates forecasting of reservoir performance and informed decision-making under various production scenarios. Although reservoir models are used to make such forecasts, neglecting geological uncertainties and history matching can limit the understanding of reservoir production behavior. While performing reservoir simulation on multiple models for uncertainty quantification is a direct approach, it is associated with significant time and computational costs and low convergence rates. This study proposes a method to address these limitations using Bayes' theorem integrated with Machine Learning (ML) algorithms, sampling methods (Markov Chain Monte Carlo and rejection), and Structured Expert Judgment (SEJ) based on the Cooke's Model (CM). The tuned ANN model outperforms traditional simulation methods in terms of computational cost when predicting posterior probability distributions, while also displaying high accuracy as confirmed by R2-score and k-fold cross-validation analyses. Furthermore, using trained ANN models in conjunction with representative samples drawn from posterior probability distributions enables the rigorous prediction of reservoir production parameters. Finally, the CM is implemented to generate a robust model by combining the statistical quantiles, including P10–P25–P50–P75–P90 estimations obtained from MCMC and rejection sampling. The applicability of this methodology is demonstrated for Teal South reservoir with eight uncertain parameters.