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

In this paper, we address the numerical posterior error control problem for the Bayesian approach to inverse problems or recently known as Bayesian Uncertainty Quantification (UQ). We generalize the results of Capistrán et al. (2016) to (a priori) expected Bayes factors (BF) and in a more general, infinite-dimensional setting. In this inverse problem, the unavoidable numerical approximation of the Forward Map (FM, i.e., the regressor function), arising from the numerical solution of a system of differential equations, demands error estimates of the corresponding approximate numerical posterior distribution. Our approach is to make such comparisons in the setting of Bayesian model selection and BFs. The main result of this paper is a bound on the absolute global error tolerated by the numerical solver of the FM in order to keep the BF of the numerical versus the theoretical posterior near one. For two examples, we provide a detailed analysis of the computation and implementation of the introduced bound. Furthermore, we show that the resulting numerical posterior turns out to be nearly identical from the theoretical posterior, given the control of the BF near one.

Women in Colombia face difficulties related to the patriarchal traits of their societies and well-known conflict afflicting the country since 1948. In this critical context, our aim is to study the relationship between baseline socio-demographic factors and variables associated to fertility, partnership patterns, and work activity. To best exploit the explanatory structure, we propose a Bayesian multivariate density regression model, which can accommodate mixed responses with censored, constrained, and binary traits. The flexible nature of the models allows for nonlinear regression functions and non-standard features in the errors, such as asymmetry or multi-modality. The model has interpretable covariate-dependent weights constructed through normalization, allowing for combinations of categorical and continuous covariates. Computational difficulties for inference are overcome through an adaptive truncation algorithm combining adaptive Metropolis-Hastings and sequential Monte Carlo to create a sequence of automatically truncated posterior mixtures. For our study on Colombian women’s life patterns, a variety of quantities are visualised and described, and in particular, our findings highlight the detrimental impact of family violence on women’s choices and behaviors.

We introduce a Bayesian non-parametric spatial factor analysis model with spatial dependency induced through a prior on factor loadings. For each column of the loadings matrix, spatial dependency is encoded using a probit stick-breaking process (PSBP) and a multiplicative gamma process shrinkage prior is used across columns to adaptively determine the number of latent factors. By encoding spatial information into the loadings matrix, meaningful factors are learned that respect the observed neighborhood dependencies, making them useful for assessing rates over space. Furthermore, the spatial PSBP prior can be used for clustering temporal trends, allowing users to identify regions within the spatial domain with similar temporal trajectories, an important task in many applied settings. In the manuscript, we illustrate the model’s performance in simulated data, but also in two real-world examples: longitudinal monitoring of glaucoma and malaria surveillance across the Peruvian Amazon. The R package spBFA, available on CRAN, implements the method.

The analysis of rank ordered data has a long history in the statistical literature across a diverse range of applications. In this paper we consider the Extended Plackett-Luce model that induces a flexible (discrete) distribution over permutations. The parameter space of this distribution is a combination of potentially high-dimensional discrete and continuous components and this presents challenges for parameter interpretability and also posterior computation. Particular emphasis is placed on the interpretation of the parameters in terms of observable quantities and we propose a general framework for preserving the mode of the prior predictive distribution. Posterior sampling is achieved using an effective simulation based approach that does not require imposing restrictions on the parameter space. Working in the Bayesian framework permits a natural representation of the posterior predictive distribution and we draw on this distribution to make probabilistic inferences and also to identify potential lack of model fit. The flexibility of the Extended Plackett-Luce model along with the effectiveness of the proposed sampling scheme are demonstrated using several simulation studies and real data examples.

Factor models aim to describe a dependence structure among high-dimensional random variables in terms of a low-dimensional unobserved random vector called a factor. One of the major practical issues of applying the factor model is to determine the factor dimensionality. In this paper, we propose a computationally feasible nonparametric prior distribution which achieves the posterior consistency of the factor dimensionality. We also derive the posterior contraction rate of the covariance matrix which is optimal when the factor dimensionality of the true covariance matrix is bounded. We conduct numerical studies that illustrate our theoretical results.

For the discovery of regression relationships between Y and a large set of p potential predictors $x_1,\dots ,x_p$, the flexible nonparametric nature of BART (Bayesian Additive Regression Trees) allows for a much richer set of possibilities than restrictive parametric approaches. However, subject matter considerations sometimes warrant a minimal assumption of monotonicity in at least some of the predictors. For such contexts, we introduce mBART, a constrained version of BART that can flexibly incorporate monotonicity in any predesignated subset of predictors using a multivariate basis of monotone trees, while avoiding the further confines of a full parametric form. For such monotone relationships, mBART provides (i) function estimates that are smoother and more interpretable, (ii) better out-of-sample predictive performance, and (iii) less post-data uncertainty. While many key aspects of the unconstrained BART model carry over directly to mBART, the introduction of monotonicity constraints necessitates a fundamental rethinking of how the model is implemented. In particular, the original BART Markov Chain Monte Carlo algorithm relied on a conditional conjugacy that is no longer available in a monotonically constrained space. Various simulated and real examples demonstrate the wide ranging potential of mBART.

Global-local shrinkage priors have been recognized as a useful class of priors that can strongly shrink small signals toward prior means while keeping large signals unshrunk. Although such priors have been extensively discussed under Gaussian responses, in practice, we often encounter count responses. Previous contributions on global-local shrinkage priors cannot be readily applied to count data. In this paper, we discuss global-local shrinkage priors for analyzing a sequence of counts. We provide sufficient conditions under which the posterior mean is unshrunk for very large signals, known as the tail robustness property. Then, we propose tractable priors to satisfy those conditions approximately or exactly and develop a custom posterior computation algorithm for Bayesian inference without tuning parameters. We demonstrate the proposed methods through simulation studies and an application to a real dataset.

We study frequentist properties of Bayesian and $L_0$ model selection, with a focus on (potentially non-linear) high-dimensional regression. We propose a construction to study how posterior probabilities and normalized $L_0$ criteria concentrate on the (Kullback-Leibler) optimal model and other subsets of the model space. When such concentration occurs, one also bounds the frequentist probabilities of selecting the correct model, type I and type II errors. These results hold generally, and help validate the use of posterior probabilities and $L_0$ criteria to control frequentist error probabilities associated to model selection and hypothesis tests. Regarding regression, we help understand the effect of the sparsity imposed by the prior or the $L_0$ penalty, and of problem characteristics such as the sample size, signal-to-noise, dimension and true sparsity. A particular finding is that one may use less sparse formulations than would be asymptotically optimal, but still attain consistency and often also significantly better finite-sample performance. We also prove new results related to misspecifying the mean or covariance structures, and give tighter rates for certain non-local priors than currently available.

Many applications in Bayesian statistics are extremely computationally intensive. However, they are often inherently parallel, making them prime targets for modern massively parallel processors. Multi-core and distributed computing is widely applied in the Bayesian community, however, very little attention has been given to fine-grain parallelisation using single instruction multiple data (SIMD) operations that are available on most modern CPUs. In this work, we practically demonstrate, using standard programming libraries, the utility of the SIMD approach for several topical Bayesian applications. Using the C programming language, we show that SIMD can improve the single-core floating point arithmetic performance by up to a factor of $6\times$ compared scalar C code and more than $25\times$ compared with optimised R code. Such improvements are multiplicative to any gains achieved through multi-core processing. We illustrate the potential of SIMD for accelerating Bayesian computations and provide the reader with techniques for exploiting modern massively parallel processing environments.

Non-linear hierarchical models are commonly used in many disciplines. However, inference in the presence of non-nested effects and on large datasets is challenging and computationally burdensome. This paper provides two contributions to scalable and accurate inference. First, I derive a new mean-field variational algorithm for estimating binomial logistic hierarchical models with an arbitrary number of non-nested random effects. Second, I propose “marginally augmented variational Bayes” (MAVB) that further improves the initial approximation through a step of Bayesian post-processing. I prove that MAVB provides a guaranteed improvement in the approximation quality at low computational cost and induces dependencies that were assumed away by the initial factorization assumptions.

I apply these techniques to a study of voter behavior using a high-dimensional application of the popular approach of multilevel regression and post-stratification (MRP). Existing estimation took hours whereas the algorithms proposed run in minutes. The posterior means are well-recovered even under strong factorization assumptions. Applying MAVB further improves the approximation by partially correcting the under-estimated variance. The proposed methodology is implemented in an open source software package.

We take a new look at the problem of disentangling the volatility and jumps processes of daily stock returns. We first provide a computational framework for the univariate stochastic volatility model with Poisson-driven jumps that offers a competitive inference alternative to the existing tools. This methodology is then extended to a large set of stocks for which we assume that their unobserved jump intensities co-evolve in time through a dynamic factor model. To evaluate the proposed modelling approach we conduct out-of-sample forecasts and we compare the posterior predictive distributions obtained from the different models. We provide evidence that joint modelling of jumps improves the predictive ability of the stochastic volatility models.