MCMC Chains Research Articles

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.

A Bayesian proportional hazards mixture cure model for interval-censored data.

The proportional hazards mixture cure model is a popular analysis method for survival data where a subgroup of patients are cured. When the data are interval-censored, the estimation of this model is challenging due to its complex data structure. In this article, we propose a computationally efficient semiparametric Bayesian approach, facilitated by spline approximation and Poisson data augmentation, for model estimation and inference with interval-censored data and a cure rate. The spline approximation and Poisson data augmentation greatly simplify the MCMC algorithm and enhance the convergence of the MCMC chains. The empirical properties of the proposed method are examined through extensive simulation studies and also compared with the R package "GORCure". The use of the proposed method is illustrated through analyzing a data set from the Aerobics Center Longitudinal Study.

Phylogenetic Relationship and Molecular Divergence Dating Using SRY Gene Polymorphism about Four Ladoum Sheep Lineages in Senegal

Animal genetic resources are playing a vital role in livestock production and are essential to food security. The present study aims to contribute to a better understanding genetic local sheep breeds and to elucidate the phylogenetic relationships through the evolution of the SRY gene in four different lineages of Ladoum sheep raised in Senegal. After a brief analysis of genetic diversity, the phylogenetic relationships and molecular dating were inferred through haplotype networks and four phylogenetic reconstruction methods. The different haplotype networks are constructed with NETWORK ver. 5.0.0.0 using the Median-Joining method. Phylogenetic trees were reconstructed using neighbor-joining, maximum parsimony, maximum likelihood and Bayesian inference. The robustness of the nodes in phylogenetic trees of the three first methods was assessed by 1000 bootstraps. For Bayesian inference, the posterior probability distribution of the trees was estimated by 4 MCMC chains. 5,000,000 generations were performed for each of the chains by sampling the different parameters every 1000 generations. Results show a low polymorphism. Haplotypic diversity is much higher than the average nucleotide divergence between all pairs of haplotypes. The majority and central haplotype indicates a close relationship between “Batling” and “Tyson” individuals. “Birahim” lineage is very distinct from the rest. Phylogenetic trees confirm two genetically separate clades between “Birahim” and the other lineages. The period of divergence between “Birahim” lineage versus the common ancestor of the other three lineages was 2504 years ago. The polyphyly revealed in “Birahim” lindicates that this lineage does not contain the common ancestor of all individuals who compose it. It could therefore be derived from two or more sheep breeds with a common ancestor, Ovis aries. The monophyletic clade appears to be a group including a common ancestor and all of its genetic descendants. This group, bringing together the other three lineages, is in the process of being structured into sub-lineages. This study is the first to show that there are only two genetic lines within ladoum sheep in Senegal.

Dating Microbial Evolution with MCMCtree.

This protocol explains how to use the program MCMCtree to estimate divergence times in microbial phylogenies. The main advantage of MCMCtree is the implementation of an approximation to the molecular data likelihood that dramatically speeds up computation during Bayesian MCMC sampling of divergence times and evolutionary rates. The approximation allows the analysis of large phylogenies with hundreds of taxa and molecular alignments with thousands or millions of sites. Two examples are used to illustrate Bayesian clock dating with MCMCtree. The first is a phylogeny of (mostly) microbial eukaryotes and prokaryotes encompassing the major groups of life on Earth, and for which fossil information, to calibrate the nodes of the phylogeny, is available. The second is a phylogeny of influenza viruses with known sampling times. An overview of Bayesian MCMC sampling is given as well as practical advice on issues such as construction of the time and rate prior and assessment of convergence of MCMC chains. Strategies for estimating times in microbial phylogenies for which neither fossil information nor sampling times are known are discussed.

Anytime parallel tempering

Developing efficient MCMC algorithms is indispensable in Bayesian inference. In parallel tempering, multiple interacting MCMC chains run to more efficiently explore the state space and improve performance. The multiple chains advance independently through local moves, and the performance enhancement steps are exchange moves, where the chains pause to exchange their current sample amongst each other. To accelerate the independent local moves, they may be performed simultaneously on multiple processors. Another problem is then encountered: depending on the MCMC implementation and inference problem, local moves can take a varying and random amount of time to complete. There may also be infrastructure-induced variations, such as competing jobs on the same processors, which arises in cloud computing. Before exchanges can occur, all chains must complete the local moves they are engaged in to avoid introducing a potentially substantial bias (Proposition 1). To solve this issue of randomly varying local move completion times in multi-processor parallel tempering, we adopt the Anytime Monte Carlo framework of (Murray, L. M., Singh, S., Jacob, P. E., and Lee, A.: Anytime Monte Carlo. arXiv preprintarXiv:1612.03319, (2016): we impose real-time deadlines on the parallel local moves and perform exchanges at these deadlines without any processor idling. We show our methodology for exchanges at real-time deadlines does not introduce a bias and leads to significant performance enhancements over the naïve approach of idling until every processor’s local moves complete. The methodology is then applied in an ABC setting, where an Anytime ABC parallel tempering algorithm is derived for the difficult task of estimating the parameters of a Lotka–Volterra predator-prey model, and similar efficiency enhancements are observed.

Multimodal Bayesian registration of noisy functions using Hamiltonian Monte Carlo

Functional data registration is a necessary processing step for many applications. The observed data can be inherently noisy, often due to measurement error or natural process uncertainty; which most functional alignment methods cannot handle. A pair of functions can also have multiple optimal alignment solutions, which is not addressed in current literature. In this paper, a flexible Bayesian approach to functional alignment is presented, which appropriately accounts for noise in the data without any pre-smoothing required. Additionally, by running parallel MCMC chains, the method can account for multiple optimal alignments via the multi-modal posterior distribution of the warping functions. To most efficiently sample the warping functions, the approach relies on a modification of the standard Hamiltonian Monte Carlo to be well-defined on the infinite-dimensional Hilbert space. This flexible Bayesian alignment method is applied to both simulated data and real data sets to show its efficiency in handling noisy functions and successfully accounting for multiple optimal alignments in the posterior; characterizing the uncertainty surrounding the warping functions.

Performance of a method for weighting a range in the number of contributors in probabilistic genotyping

Uncertainty in the assignment of the number of contributors (NoC) can be encountered, particularly in higher-order mixtures, where alleles may be shared between contributors, may have dropped out, or may be masked by the stutter artefacts or allelic peaks of a more dominant contributor. Most probabilistic genotyping software requires the assignment of NoC prior to interpretation. NoC has been described as a nuisance parameter. Taylor et al. [1] describe a method to weigh the probability of the profile under different values of N and incorporate this into a likelihood ratio (LR). Within this paper we explore the performance of this variable number of contributors (varNoC) method programmed within the probabilistic genotyping software STRmix™. The desired combination of performance and runtime was obtained using the default STRmix™ version 2.7 MCMC settings in conjunction with a 2.5 % hyper-rectangle range, at least 10,000 naïve MC iterations and 8 MCMC chains. The varNoC LR demonstrated the typical sensitivity and specificity behaviour seen in previous studies, with a high level of reproducibility given repeat analyses. Profiles previously demonstrating ambiguity in the NoC assigned using conventional estimation methods, were able to be reliably interpreted and a varNoC LR assigned.

Corrigendum.

Corrigendum.

Hamiltonian Markov chain Monte Carlo for partitioned sample spaces with application to Bayesian deep neural nets

Allocating computation over multiple chains to reduce sampling time in MCMC is crucial in making MCMC more applicable in the state of the art models such as deep neural networks. One of the parallelization schemes for MCMC is partitioning the sample space to run different MCMC chains in each component of the partition (VanDerwerken and Schmidler in Parallel Markov chain Monte Carlo. arXiv:1312.7479, 2013; Basse et al. in Artificial intelligence and statistics, pp 1318–1327, 2016). In this work, we take Basse et al. (2016)’s bridge sampling approach and apply constrained Hamiltonian Monte Carlo on partitioned sample spaces. We propose a random dimension partition scheme that combines well with the constrained HMC. We empirically show that this approach can expedite MCMC sampling for any unnormalized target distribution such as Bayesian neural network in a high dimensional setting. Furthermore, in the presence of multi-modality, this algorithm is expected to be more efficient in mixing MCMC chains when proper partition elements are chosen.

Mixing of MCMC algorithms

ABSTRACTWe analyse MCMC chains focusing on how to find simulation parameters that give good mixing for discrete time, Harris ergodic Markov chains on a general state space X having invariant distribution π. The analysis uses an upper bound for the variance of the probability estimate. For each simulation parameter set, the bound is estimated from an MCMC chain using recurrence intervals. Recurrence intervals are a generalization of recurrence periods for discrete Markov chains. It is easy to compare the mixing properties for different simulation parameters. The paper gives general advice on how to improve the mixing of the MCMC chains and a new methodology for how to find an optimal acceptance rate for the Metropolis-Hastings algorithm. Several examples, both toy examples and large complex ones, illustrate how to apply the methodology in practice. We find that the optimal acceptance rate is smaller than the general recommendation in the literature in some of these examples.

Current constraints on anisotropic and isotropic dark energy models

We use Gaussian processes in combination with MCMC method to place constraints on cosmological parameters of three dark energy models including flat and curved FRW and Bianchi type I spacetimes. To do so, we use recently compiled 36 measurements of the Hubble parameter $H(z)$ in the redshifts intermediate $0.07\leqslant z \leqslant 2.36$. Moreover, we use these models to estimate the redshift of the deceleration-acceleration transition. We consider two Gaussian priors for current value of the Hubble constant i.e $H_{0}=73\pm1.74 (68\pm 2.8)$ km/s/Mpc to investigate the effect of the assumed $H_{0}$ on our parameters estimations. For statistical analysis we use NUTS sampler which is an extension of Hamiltonian Monte Carlo algorithm to generate MCMC chains for parameters of dark energy models. To compare the considered cosmologies, we perform Akaike information criterion (AIC) and Bayes factor ($\Psi$). In general, when we compared our results with 9 years WMAP as well as Planck 2015 Collaboration, we found that Bianchi type I model is slightly fits better to the observational Hubble data with respect to the non-flat FRW model.

Bayesian learning of weakly structural Markov graph laws using sequential Monte Carlo methods

We present a sequential sampling methodology for weakly structural Markov laws, arising naturally in a Bayesian structure learning context for decomposable graphical models. As a key component of our suggested approach, we show that the problem of graph estimation, which in general lacks natural sequential interpretation, can be recast into a sequential setting by proposing a recursive Feynman-Kac model that generates a flow of junction tree distributions over a space of increasing dimensions. We focus on particle McMC methods to provide samples on this space, in particular on particle Gibbs (PG), as it allows for generating McMC chains with global moves on an underlying space of decomposable graphs. To further improve the PG mixing properties, we incorporate a systematic refreshment step implemented through direct sampling from a backward kernel. The theoretical properties of the algorithm are investigated, showing that the proposed refreshment step improves the performance in terms of asymptotic variance of the estimated distribution. The suggested sampling methodology is illustrated through a collection of numerical examples demonstrating high accuracy in Bayesian graph structure learning in both discrete and continuous graphical models.

Prior Precision, Prior Accuracy, and the Estimation of Disease Prevalence Using Imperfect Diagnostic Tests.

Estimates of disease prevalence in any host population are complicated by uncertainty in the outcome of diagnostic tests on individuals. In the absence of gold standard diagnostics (tests that give neither false positives nor false negatives), Bayesian latent class inference can be applied to batteries of diagnostic tests, providing posterior estimates of the sensitivity and specificity of each test, alongside posterior estimates of disease prevalence. Here we explore the influence of precision and accuracy of prior information on the precision and accuracy of posterior estimates of these key parameters. Our simulations use three diagnostic tests, yielding eight possible diagnostic outcomes per individual. Seven degrees of freedom allow the estimation of seven parameters: sensitivity and specificity of each test, and disease prevalence. We show that prior precision begets posterior precision but only when priors are accurate. We also show that analyses without gold standard can use imprecise priors as long as they are initialised with accuracy. Imprecise priors risk the divergence of MCMC chains towards inaccurate posterior estimates, if inaccurate initial values are used. We note that inaccurate priors can yield inaccurate and imprecise inference. Bounded priors should certainly not be used unless their accuracy is well established. Inaccurate estimates of sensitivity or specificity can yield wildly inaccurate estimates of disease prevalence. Our analyses are motivated by studies of bovine tuberculosis in a wild badger population.

THE PAndAS VIEW OF THE ANDROMEDA SATELLITE SYSTEM. II. DETAILED PROPERTIES OF 23 M31 DWARF SPHEROIDAL GALAXIES

ABSTRACT We present a comprehensive analysis of the structural properties and luminosities of the 23 dwarf spheroidal galaxies that fall within the footprint of the Pan-Andromeda Archaeological Survey (PAndAS). These dwarf galaxies represent the large majority of Andromeda’s known satellite dwarf galaxies and cover a wide range in luminosity ( or ) and surface brightness ( mag arcsec−2). We confirm most previous measurements, but we find And XIX to be significantly larger than before ( , ) and cannot derive parameters for And XXVII as it is likely not a bound stellar system. We also significantly revise downward the luminosities of And XV and And XVI, which are now or . Finally, we provide the first detailed analysis of Cas II/And XXX, a fairly faint system ( ) of typical size ( ), located in close proximity to the two bright elliptical dwarf galaxies NGC 147 and NGC 185. Combined with the set of homogeneous distances published in an earlier contribution, our analysis dutifully tracks all relevant sources of uncertainty in the determination of the properties of the dwarf galaxies from the PAndAS photometric catalog. We further publish the posterior probability distribution functions of all the parameters we fit for in the form of MCMC chains available online; these inputs should be used in any analysis that aims to remain truthful to the data and properly account for covariance between parameters.

Bayesian EEG source localization using a structured sparsity prior

This paper deals with EEG source localization. The aim is to perform spatially coherent focal localization and recover temporal EEG waveforms, which can be useful in certain clinical applications. A new hierarchical Bayesian model is proposed with a multivariate Bernoulli Laplacian structured sparsity prior for brain activity. This distribution approximates a mixed ℓ20 pseudo norm regularization in a Bayesian framework. A partially collapsed Gibbs sampler is proposed to draw samples asymptotically distributed according to the posterior of the proposed Bayesian model. The generated samples are used to estimate the brain activity and the model hyperparameters jointly in an unsupervised framework. Two different kinds of Metropolis–Hastings moves are introduced to accelerate the convergence of the Gibbs sampler. The first move is based on multiple dipole shifts within each MCMC chain, whereas the second exploits proposals associated with different MCMC chains. Experiments with focal synthetic data shows that the proposed algorithm is more robust and has a higher recovery rate than the weighted ℓ21 mixed norm regularization. Using real data, the proposed algorithm finds sources that are spatially coherent with state of the art methods, namely a multiple sparse prior approach and the Champagne algorithm. In addition, the method estimates waveforms showing peaks at meaningful timestamps. This information can be valuable for activity spread characterization.

Taxonomic recommendations for British birds: seventh report

Taxonomic recommendations for British birds: seventh report

Using MCMC chain outputs to efficiently estimate Bayes factors

One of the most important methodological problems in psychological research is assessing the reasonableness of null models, which typically constrain a parameter to a specific value such as zero. Bayes factor has been recently advocated in the statistical and psychological literature as a principled means of measuring the evidence in data for various models, including those where parameters are set to specific values. Yet, it is rarely adopted in substantive research, perhaps because of the difficulties in computation. Fortunately, for this problem, the Savage–Dickey density ratio ( Dickey & Lientz, 1970) provides a conceptually simple approach to computing Bayes factor. Here, we review methods for computing the Savage–Dickey density ratio, and highlight an improved method, originally suggested by Gelfand and Smith (1990) and advocated by Chib (1995), that outperforms those currently discussed in the psychological literature. The improved method is based on conditional quantities, which may be integrated by Markov chain Monte Carlo sampling to estimate Bayes factors. These conditional quantities efficiently utilize all the information in the MCMC chains, leading to accurate estimation of Bayes factors. We demonstrate the method by computing Bayes factors in one-sample and one-way designs, and show how it may be implemented in WinBUGS.

A note on bimodality in the log-likelihood function for penalized spline mixed models

For a smoothing spline or general penalized spline model, the smoothing parameter can be estimated using residual maximum likelihood (REML) methods by expressing the spline in the form of a mixed model. The possibility of bimodality in the profile log-likelihood function for the smoothing parameter of these penalized spline mixed models is demonstrated. A canonical transformation into independent observations is used to provide efficient evaluation of the log-likelihood function and gives insight into the incompatibilities between the model and data that cause bimodality. This transformation can also be used to assess the influence of different frequency components in the data on the estimated smoothing parameter. It is demonstrated that, where bimodality occurs in the log-likelihood, Bayesian penalized spline models may show poor mixing in MCMC chains and be sensitive to the choice of prior distributions for variance components.

Does a Gibbs sampler approach to spatial Poisson regression models outperform a single site MH sampler?

A Gibbs sampler for a Poisson regression model including spatial effects is presented and evaluated. The approach is based on that a Poisson regression model can be transformed into an approximate normal linear model by data augmentation using the introduction of two sequences of latent variables. It is shown how this methodology can be extended to spatial Poisson regression models and details of the resulting Gibbs sampler are given. In particular, the influence of model parameterisation and different update strategies on the mixing of the MCMC chains is discussed. The developed Gibbs samplers are analysed in two simulation studies and applied to model the expected number of claims for policyholders of a German car insurance company. The mixing of the Gibbs samplers depends crucially on the model parameterisation and the update schemes. The best mixing is achieved when collapsed algorithms are used, reasonable low autocorrelations for the spatial effects are obtained in this case. For the regression effects however, autocorrelations are rather high, especially for data with very low heterogeneity. For comparison a single component Metropolis–Hastings algorithms is applied which displays very good mixing for all components. Although the Metropolis–Hastings sampler requires a higher computational effort, it outperforms the Gibbs samplers which would have to be run considerably longer in order to obtain the same precision of the parameters.

On learning strategies for evolutionary Monte Carlo

The real-parameter evolutionary Monte Carlo algorithm (EMC) has been proposed as an effective tool both for sampling from high-dimensional distributions and for stochastic optimization (Liang and Wong, 2001). EMC uses a temperature ladder similar to that in parallel tempering (PT; Geyer, 1991). In contrast with PT, EMC allows for crossover moves between parallel and tempered MCMC chains. In the context of EMC, we introduce four new moves, which enhance its efficiency as measured by the effective sample size. Secondly, we introduce a practical strategy for determining the temperature range and placing the temperatures in the ladder used in EMC and PT. Lastly, we prove the validity of the conditional sampling step of the snooker algorithm, a crossover move in EMC, which extends a result of Roberts and Gilks (1994).