MCMC Sampling Algorithm Research Articles

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.

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.

Joint classification and prediction of random curves using heavy‐tailed process functional regression

We propose a heavy-tailed process functional regression to jointly perform classification and prediction of time-varying functional data. We use two independent scale mixtures of Gaussian Processes to respectively model random effects and random errors, yielding robust inferences against both magnitude and shape outliers. We classify random curves by posterior predictive probabilities of class labels and offer a weighted prediction of future curve trends. A Bayesian estimation procedure is implemented through an MCMC sampling algorithm. The performance of classification and prediction of the proposed model is evaluated using simulated studies and some real data sets.

Direct Estimation of Diagnostic Classification Model Attribute Mastery Profiles via a Collapsed Gibbs Sampling Algorithm.

This paper proposes a novel collapsed Gibbs sampling algorithm that marginalizes model parameters and directly samples latent attribute mastery patterns in diagnostic classification models. This estimation method makes it possible to avoid boundary problems in the estimation of model item parameters by eliminating the need to estimate such parameters. A simulation study showed the collapsed Gibbs sampling algorithm can accurately recover the true attribute mastery status in various conditions. A second simulation showed the collapsed Gibbs sampling algorithm was computationally more efficient than another MCMC sampling algorithm, implemented by JAGS. In an analysis of real data, the collapsed Gibbs sampling algorithm indicated good classification agreement with results from a previous study.

Bayesian Estimation of Adsorption and Desorption Parameters for Pore Scale Transport

Stochastic parameter estimation and inversion have become increasingly popular in recent years. Nowadays, it is computationally reasonable and regular to solve complex inverse problems within the Bayesian framework. Applications of Bayesian inferences for inverse problems require investigation of the posterior distribution, which usually has a complex landscape and is highly dimensional. In these cases, Markov chain Monte Carlo methods (MCMC) are often used. This paper discusses a Bayesian approach for identifying adsorption and desorption rates in combination with a pore-scale reactive flow. Markov chain Monte Carlo sampling is used to estimate adsorption and desorption rates. The reactive transport in porous media is governed by incompressible Stokes equations, coupled with convection–diffusion equation for species’ transport. Adsorption and desorption are accounted via Robin boundary conditions. The Henry isotherm is considered for describing the reaction terms. The measured concentration at the outlet boundary is provided as additional information for the identification procedure. Metropolis–Hastings and Adaptive Metropolis algorithms are implemented. Credible intervals have been plotted from sampled posterior distributions for both algorithms. The impact of the noise in the measurements and influence of several measurements for Bayesian identification procedure is studied. Sample analysis using the autocorrelation function and acceptance rate is performed to estimate mixing of the Markov chain. As result, we conclude that MCMC sampling algorithm within the Bayesian framework is good enough to determine an admissible set of parameters via credible intervals.

Bayesian auxiliary variable model for birth records data with qualitative and quantitative responses

Many applications involve data with qualitative and quantitative responses. When there is an association between the two responses, a joint model will produce improved results than fitting them separately. In this paper, a Bayesian method is proposed to jointly model such data. The joint model links the qualitative and quantitative responses and can assess their dependency strength via a latent variable. The posterior distributions of parameters are obtained through an efficient MCMC sampling algorithm. The simulation is conducted to show that the proposed method improves the prediction capacity for both responses. Further, the proposed joint model is applied to the birth records data acquired by the Virginia Department of Health for studying the mutual dependence between preterm birth of infants and their birth weights.

An Almost Constant Lower Bound of the Isoperimetric Coefficient in the KLS Conjecture

We prove an almost constant lower bound of the isoperimetric coefficient in the KLS conjecture. The lower bound has the dimension dependency d^{-o_d(1)}. When the dimension is large enough, our lower bound is tighter than the previous best bound which has the dimension dependency d^{-1/4}. Improving the current best lower bound of the isoperimetric coefficient in the KLS conjecture has many implications, including improvements of the current best bounds in Bourgain’s slicing conjecture and in the thin-shell conjecture, better concentration inequalities for Lipschitz functions of log-concave measures and better mixing time bounds for MCMC sampling algorithms on log-concave measures.

Point-Collocation NIPC 기법 기반 천이 난류 모델 계수의 베이지안 추론 연구

In the present work, Bayesian inference based on the point-collocation Non-intrusive Polynomial Chaos was conducted in transitional flow around the flat plate using Menter’s γ-Re<SUB>θ</SUB> transition model. Three model coefficients, C<SUB>a2</SUB>, C<SUB>e1</SUB> and C<SUB>e2</SUB> were considered for random variables based on the sensitivity analysis of the present turbulence model and quantity of interest was set to the drag coefficient, C<SUB>d</SUB>. With the assumption of uniform distribution of three model coefficients within ±10%, the surrogate model was obtained based on the general Polynomial Chaos expansion. The simulation results were used for observation to calculate the likelihood function and MCMC sampling algorithm was adopted for Bayesian inference. The correlation between C<SUB>e1</SUB> and C<SUB>d</SUB> was predicted with highest value than other two model coefficients. Based on the original Bayesian inference result, the effect of number of observation data was studied. Also the 2nd order of gPC for surrogate model was adequate through comparison with results of 3rd order gPC. Finally, two prior distributions of input random variables for Bayesian inference were considered and the differences of posterior distributions of input and output were investigated.

Analysis of a multilevel Markov chain Monte Carlo finite element method for Bayesian inversion of log-normal diffusions

We develop the multilevel Markov chain Monte Carlo finite element method (MLMCMC-FEM) to sample from the posterior density of the Bayesian inverse problems. The unknown is the diffusion coefficient of a linear, second-order divergence form, elliptic equation in a bounded, polytopal subdomain of . We provide a convergence analysis with absolute mean convergence rate estimates for the proposed modified MLMCMC-FEM showing in particular error versus work bounds, which are explicit in the discretization parameters. This work generalizes the MLMCMC-FEM algorithm and the error versus work analysis for the uniform prior measure from Hoang et al (2013 Inverse Problems 29), which we also review here, to linear, elliptic, divergence-form PDEs with a log-Gaussian uncertain coefficient and Gaussian prior measure. In comparison to Hoang et al (2013 Inverse Problems 29), we show by mathematical proofs and numerical examples that the unboundedness of the parameter range under Gaussian prior and the non-uniform ellipticity of the forward model require essential modifications in the MCMC sampling algorithm and in the error analysis. The proposed novel multilevel MCMC sampler applies to general Bayesian inverse problems for linear, second order elliptic divergence-form PDEs with log-Gaussian coefficients. It only requires a numerical forward solver with essentially optimal complexity for producing an approximation of the posterior expectation of a quantity of interest within a prescribed accuracy. Numerical examples using independence and pCN samplers are in agreement with our error versus work analysis.

Fast MCMC sampling algorithms on polytopes

We propose and analyze two new MCMC sampling algorithms, the Vaidya walk and the John walk, for generating samples from the uniform distribution over a polytope. Both random walks are sampling algorithms derived from interior point methods. The former is based on volumetric-logarithmic barrier introduced by Vaidya whereas the latter uses John's ellipsoids. We show that the Vaidya walk mixes in significantly fewer steps than the logarithmic-barrier based Dikin walk studied in past work. For a polytope in ℝd defined by n ≥ d linear constraints, we show that the mixing time from a warm start is bounded as O(n0.5d1.5), compared to the O(nd) mixing time bound for the Dikin walk. The cost of each step of the Vaidya walk is of the same order as the Dikin walk, and at most twice as large in terms of constant pre-factors. For the John walk, we prove an O(d2.5 log4(n/d) bound on its mixing time and conjecture that an improved variant of it could achieve a mixing time of O(d2 poly-log(n/d)). Additionally, we propose variants of the Vaidya and John walks that mix in polynomial time from a deterministic starting point. The speed-up of the Vaidya walk over the Dikin walk are illustrated in numerical examples.

Bayesian estimation and hypothesis tests for a circular Generalized Linear Model

Motivated by a study from cognitive psychology, we develop a Generalized Linear Model for circular data within the Bayesian framework, using the von Mises distribution. Although circular data arise in a wide variety of scientific fields, the number of methods for their analysis is limited. Our model allows inclusion of both continuous and categorical covariates. In a frequentist setting, this type of model is plagued by the likelihood surface of its regression coefficients, which is not logarithmically concave. In a Bayesian context, a weakly informative prior solves this issue, while for other parametersnoninformative priors are available. In addition to an MCMC sampling algorithm, we develop Bayesian hypothesis tests based on the Bayes factor for both equality and inequality constrained hypotheses. In a simulation study, it can be seen that our method performs well. The analyses are available in the package CircGLMBayes. Finally, we apply this model to a dataset from experimental psychology, and show that it provides valuable insight for applied researchers. Extensions to dependent observations are within reach by means of the multivariate von Mises distribution.

Efficient orbit-aware triad and quad census in directed and undirected graphs

The prevalence of select substructures is an indicator of network effects in applications such as social network analysis and systems biology. Moreover, subgraph statistics are pervasive in stochastic network models, and they need to be assessed repeatedly in MCMC sampling and estimation algorithms. We present a new approach to count all induced and non-induced four-node subgraphs (the quad census) on a per-node and per-edge basis, complete with a separation into their non-automorphic roles in these subgraphs. It is the first approach to do so in a unified manner, and is based on only a clique-listing subroutine. Computational experiments indicate that, despite its simplicity, the approach outperforms previous, less general approaches.By way of the more presentable triad census, we additionally show how to extend the quad census to directed graphs. As a byproduct we obtain the asymptotically fastest triad census algorithm to date.

Dynamic Quantile Function Models

We offer a novel way of thinking about the modelling of the time-varying distributions of financial asset returns. Borrowing ideas from symbolic data analysis, we consider data representations beyond scalars and vectors. Specifically, we consider a quantile function as an observation, and develop a new class of dynamic models for quantile-function-valued (QF-valued) time series. In order to make statistical inferences and account for parameter uncertainty, we propose a method whereby a likelihood function can be constructed for QF-valued data, and develop an adaptive MCMC sampling algorithm for simulating from the posterior distribution. Compared to modelling realised measures, modelling the entire quantile functions of intra-daily returns allows one to gain more insight into the dynamic structure of price movements. Via simulations, we show that the proposed MCMC algorithm is effective in recovering the posterior distribution, and that the posterior means are reasonable point estimates of the model parameters. For empirical studies, the new model is applied to analysing one-minute returns of major international stock indices. Through quantile scaling, we further demonstrate the usefulness of our method by forecasting one-step-ahead the Value-at-Risk of daily returns.

On expressiveness of the chain graph interpretations

In this article we study the expressiveness of the different chain graph interpretations. Chain graphs is a class of probabilistic graphical models that can contain two types of edges, representing different types of relationships between the variables in question. Chain graphs is also a superclass of directed acyclic graphs, i.e. Bayesian networks, and can thereby represent systems more accurately than this less expressive class of models. Today there do however exist several different ways of interpreting chain graphs and what conditional independences they encode, giving rise to different so-called chain graph interpretations. Previous research has approximated the number of representable independence models for the Lauritzen–Wermuth–Frydenberg and the multivariate regression chain graph interpretations using an MCMC based approach. In this article we use a similar approach to approximate the number of models representable by the latest chain graph interpretation in research, the Andersson–Madigan–Perlman interpretation. Moreover we summarize and compare the different chain graph interpretations with each other. Our results confirm previous results that directed acyclic graphs only can represent a small fraction of the models representable by chain graphs, even for a low number of nodes. The results also show that the Andersson–Madigan–Perlman and multivariate regression interpretations can represent about the same amount of models and twice the amount of models compared to the Lauritzen–Wermuth–Frydenberg interpretation. However, at the same time almost all models representable by the latter interpretation can only be represented by that interpretation while the former two have a large intersection in terms of representable models. • We study the expressiveness of the different chain graph interpretations in terms of representable independence models. • We show how an MCMC sampling algorithm can be used to sample chain graph models. • We show that chain graphs can represent exponentially many more independence models compared to Bayesian networks. • We show (numerically) that the average number of chain graphs per chain graph model appears to converge. • We show that AMP and MVR chain graphs are more expressive than LWF chain graphs.

Flexible paleoclimate age-depth models using an autoregressive gamma process

Radiocarbon dating is routinely used in paleoecology to build chronolo- gies of lake and peat sediments, aiming at inferring a model that would relate the sediment depth with its age. We present a new approach for chronology building (called \Bacon) that has received enthusiastic attention by paleoecologists. Our methodology is based on controlling core accumulation rates using a gamma au- toregressive semiparametric model with an arbitrary number of subdivisions along the sediment. Using prior knowledge about accumulation rates is crucial and in- formative priors are routinely used. Since many sediment cores are currently ana- lyzed, using difierent data sets and prior distributions, a robust (adaptive) MCMC is very useful. We use the t-walk (Christen and Fox, 2010), a self adjusting, robust MCMC sampling algorithm, that works acceptably well in many situations. Out- liers are also addressed using a recent approach that considers a Student-t model for radiocarbon data. Two examples are presented here, that of a peat core and a core from a lake, and our results are compared with other approaches. Past climates and environments can be reconstructed from deposits such as ocean or lake sediments, ice sheets and peat bogs. Within a vertical sediment proflle (core), mea- surements of microfossils, macrofossils, isotopes and other variables at a range of depths serve as proxy estimates or \proxies of climate and environmental conditions when the sediment of those depths was deposited. It is crucial to establish reliable relationships between these depths and their ages. Age-depth relationships are used to study the evolution of climate/environmental proxies along sediment depth and therefore through time (e.g., Lowe and Walker 1997). Age-depth models are constructed in various ways. For sediment depths containing organic matter, and for ages younger than c. 50,000 years, radiocarbon dating is often used to create an age-depth model. Cores are divided into slices and some of these are radiocarbon dated. A curve is fltted to the radiocarbon data and interpolated to obtain an age estimate for every depth of the core. The flrst restriction to be considered is that age should be increasing monotonically with depth, because sediment can never have accumulated backwards in time (extraordinary events leading to mixed or reversed sediments are, most of the time, noticeable in the stratigraphy and therefore such cores are ruled out from further analyses). Moreover, cores may have missing sections, leading to ∞at parts in the age depth models.

Nonparametric function estimation subject to monotonicity, convexity and other shape constraints

This paper uses free-knot and fixed-knot regression splines in a Bayesian context to develop methods for the nonparametric estimation of functions subject to shape constraints in models with log-concave likelihood functions. The shape constraints we consider include monotonicity, convexity and functions with a single minimum. A computationally efficient MCMC sampling algorithm is developed that converges faster than previous methods for non-Gaussian models. Simulation results indicate the monotonically constrained function estimates have good small sample properties relative to (i) unconstrained function estimates, and (ii) function estimates obtained from other constrained estimation methods when such methods exist. Also, asymptotic results show the methodology provides consistent estimates for a large class of smooth functions. Two detailed illustrations exemplify the ideas.

Testing forecasts of a new Bayesian time-predictable model of eruption occurrence

In this paper we propose a model to forecast eruptions in a real forward perspective. Specifically, the model provides a forecast of the next eruption after the end of the last one, using only the data available up to that time. We focus our attention on volcanoes with open conduit regime and high eruption frequency. We assume a generalization of the classical time predictable model to describe the eruptive behavior of open conduit volcanoes and we use a Bayesian hierarchical model to make probabilistic forecasts. We apply the model to Kilauea volcano eruptive data and Mount Etna volcano flank eruption data. The aims of the proposed model are: (1) to test whether or not the Kilauea and Mount Etna volcanoes follow a time predictable behavior; (2) to discuss the volcanological implications of the time predictable model parameters inferred; (3) to compare the forecast capabilities of this model with other models present in literature. The results obtained using the MCMC sampling algorithm show that both volcanoes follow a time predictable behavior. The numerical values inferred for the parameters of the time predictable model suggest that the amount of the erupted volume could change the dynamics of the magma chamber refilling process during the repose period. The probability gain of this model compared with other models already present in literature is appreciably greater than zero. This means that our model provides better forecast than previous models and it could be used in a probabilistic volcanic hazard assessment scheme.

Estimation of fish and wildlife disease prevalence from imperfect diagnostic tests on pooled samples with varying pool sizes

Methods of estimating disease or parasite prevalence in free-ranging and some captive fish and wildlife populations are frequently lacking in precision due to limited numbers of observations and different assay procedures. Recently statistical methods and software programs have been developed to use Bayesian and other methods to obtain estimates of disease prevalence from diagnostic tests in which sensitivity and/or specificity is not perfect (imperfect) and with sampling schemes using pooled samples. However, these published methods and software programs that consider pooled data sampling have generally considered the case of one uniform pool size for all samples. We present a method for estimating disease prevalence from imperfect diagnostic tests with pooled data collected from a variety of pool sizes. We use a Bayesian approach and obtain a sample from the posterior distribution of prevalence, sensitivity, and specificity, using an MCMC sampling algorithm implemented in the WINBUGS statistical package. We illustrate the use of these methods with three examples and perform efficiency calculations to investigate the performance of these estimators relative to maximum likelihood estimators that assume perfect diagnostic tests. Our results illustrate that the estimates produced from these methods adjust for imperfect tests, and are often more efficient than estimates assuming perfect tests, except in some situations when there is not much prior information on diagnostic test sensitivity and specificity.