A sharp first order analysis of Feynman–Kac particle models,Part II: Particle Gibbs samplers

This article provides a new theory for the analysis of the particle Gibbs (PG) sampler (Andrieu et al., 2010). Following the work of Del Moral and Jasra (2017) we provide some analysis of the particle Gibbs sampler, giving first order expansions of the kernel and minorization estimates. In addition, first order propagation of chaos estimates are derived for empirical measures of the dual particle model with a frozen path, also known as the conditional sequential Monte Carlo (SMC) update of the PG sampler. Backward and forward PG samplers are discussed, including a first comparison of the contraction estimates obtained by first order estimates. We illustrate our results with an example of fixed parameter estimation arising in hidden Markov models.     

Performance of the Gibbs,Hit-and-Run,and Metropolis Samplers

Abstract

We consider the performance of three Monte Carlo Markov-chain samplers—the Gibbs sampler, which cycles through coordinate directions; the Hit-and-Run (H&R) sampler, which randomly moves in any direction; and the Metropolis sampler, which moves with a probability that is a ratio of likelihoods. We obtain several analytical results. We provide a sufficient condition of the geometric convergence on a bounded region S for the H&R sampler. For a general region S, we review the Schervish and Carlin sufficient geometric convergence condition for the Gibbs sampler. We show that for a multivariate normal distribution this Gibbs sufficient condition holds and for a bivariate normal distribution the Gibbs marginal sample paths are each an AR(1) process, and we obtain the standard errors of sample means and sample variances, which we later use to verify empirical Monte Carlo results. We empirically compare the Gibbs and H&R samplers on bivariate normal examples. For zero correlation, the Gibbs sampler provides independent data, resulting in better performance than H&R. As the absolute value of the correlation increases, H&R performance improves, with H&R substantially better for correlations above .9. We also suggest and study methods for choosing the number of replications, for estimating the standard error of point estimators and for reducing point-estimator variance. We suggest using a single long run instead of using multiple iid separate runs. We suggest using overlapping batch statistics (obs) to get the standard errors of estimates; additional empirical results show that obs is accurate. Finally, we review the geometric convergence of the Metropolis algorithm and develop a Metropolisized H&R sampler. This sampler works well for high-dimensional and complicated integrands or Bayesian posterior densities.     

A Generalized Markov Sampler

A recent development of the Markov chain Monte Carlo (MCMC) technique is the emergence of MCMC samplers that allow transitions between different models. Such samplers make possible a range of computational tasks involving models, including model selection, model evaluation, model averaging and hypothesis testing. An example of this type of sampler is the reversible jump MCMC sampler, which is a generalization of the Metropolis–Hastings algorithm. Here, we present a new MCMC sampler of this type. The new sampler is a generalization of the Gibbs sampler, but somewhat surprisingly, it also turns out to encompass as particular cases all of the well-known MCMC samplers, including those of Metropolis, Barker, and Hastings. Moreover, the new sampler generalizes the reversible jump MCMC. It therefore appears to be a very general framework for MCMC sampling. This paper describes the new sampler and illustrates its use in three applications in Computational Biology, specifically determination of consensus sequences, phylogenetic inference and delineation of isochores via multiple change-point analysis.     

Geometric Ergodicity of Gibbs and Block Gibbs Samplers for a Hierarchical Random Effects Model

We consider fixed scan Gibbs and block Gibbs samplers for a Bayesian hierarchical random effects model with proper conjugate priors. A drift condition given in Meyn and Tweedie (1993, Chapter 15) is used to show that these Markov chains are geometrically ergodic. Showing that a Gibbs sampler is geometrically ergodic is the first step toward establishing central limit theorems, which can be used to approximate the error associated with Monte Carlo estimates of posterior quantities of interest. Thus, our results will be of practical interest to researchers using these Gibbs samplers for Bayesian data analysis.     

Sparse Partially Collapsed MCMC for Parallel Inference in Topic Models

Topic models, and more specifically the class of latent Dirichlet allocation (LDA), are widely used for probabilistic modeling of text. Markov chain Monte Carlo (MCMC) sampling from the posterior distribution is typically performed using a collapsed Gibbs sampler. We propose a parallel sparse partially collapsed Gibbs sampler and compare its speed and efficiency to state-of-the-art samplers for topic models on five well-known text corpora of differing sizes and properties. In particular, we propose and compare two different strategies for sampling the parameter block with latent topic indicators. The experiments show that the increase in statistical inefficiency from only partial collapsing is smaller than commonly assumed, and can be more than compensated by the speedup from parallelization and sparsity on larger corpora. We also prove that the partially collapsed samplers scale well with the size of the corpus. The proposed algorithm is fast, efficient, exact, and can be used in more modeling situations than the ordinary collapsed sampler. Supplementary materials for this article are available online.     

Stability Relationships Among the Gibbs Sampler and its Subchains

The use of Gibbs samplers driven by improper posteriors has been a controversial issue in the statistics literature over the last few years. It has recently been demonstrated that it is possible to make valid statistical inferences through such Gibbs samplers. Furthermore, theoretical and empirical evidence has been given to support the idea that there are actually computational advantages to using these nonpositive recurrent Markov chains rather than more standard positive recurrent chains. These results provide motivation for a general study of the behavior of the Gibbs Markov chain when it is not positive recurrent. This article concerns stability relationships among the two-variable Gibbs sampler and its subchains. We show that these three Markov chains always share the same stability; that is, they are either all positive recurrent, all null recurrent, or all transient. In addition, we establish general results concerning the ways in which positive recurrent Markov chains can arise from null recurrent and transient Gibbs chains. Six examples of varying complexity are used to illustrate the results.     

Uncertainty Quantification for Modern High-Dimensional Regression via Scalable Bayesian Methods

Tremendous progress has been made in the last two decades in the area of high-dimensional regression, especially in the “large p, small n” setting. Such sample starved settings inevitably lead to models which are potentially very unstable and hence quite unreliable. To this end, Bayesian shrinkage methods have generated a lot of recent interest in the modern high-dimensional regression and model selection context. Such methods span the wide spectrum of modern regression approaches and include among others, spike-and-slab priors, the Bayesian lasso, ridge regression, and global-local shrinkage priors such as the Horseshoe prior and the Dirichlet–Laplace prior. These methods naturally facilitate tractable uncertainty quantification and have thus been used extensively across diverse applications. A common unifying feature of these models is that the corresponding priors on the regression coefficients can be expressed as a scale mixture of normals. This property has been leveraged extensively to develop various three-step Gibbs samplers to explore the corresponding intractable posteriors. The convergence of such samplers however is very slow in high dimensions settings, making them disconnected to the very setting that they are intended to work in. To address this challenge, we propose a comprehensive and unifying framework to draw from the same family of posteriors via a class of tractable and scalable two-step blocked Gibbs samplers. We demonstrate that our proposed class of two-step blocked samplers exhibits vastly superior convergence behavior compared to the original three-step sampler in high-dimensional regimes on simulated data as well as data from a variety of applications including gene expression data, infrared spectroscopy data, and socio-economic/law enforcement data. We also provide a detailed theoretical underpinning to the new method by deriving explicit upper bounds for the (geometric) rate of convergence, and by proving that the proposed two-step sampler has superior spectral properties. Supplementary material for this article is available online.     

Twisting the Alive Particle Filter

This work focuses on sampling from hidden Markov models (Cappe et al. 2005) whose observations have intractable density functions. We develop a new sequential Monte Carlo (e.g. Doucet, 2011) algorithm and a new particle marginal Metropolis-Hastings (Andrieu et al J R Statist Soc Ser B 72:269-342, 2010) algorithm for these purposes. We build from Jasra et al (2013) and Whiteley and Lee (Ann Statist 42:115-141, 2014) to construct the sequential Monte Carlo (SMC) algorithm, which we call the alive twisted particle filter. Like the alive particle filter (Amrein and Künsch, 2011, Jasra et al, 2013), our new SMC algorithm adopts an approximate Bayesian computation (Tavare et al. Genetics 145:505-518, 1997) estimate of the HMM. Our alive twisted particle filter also uses a twisted proposal as in Whiteley and Lee (Ann Statist 42:115-141, 2014) to obtain a low-variance estimate of the HMM normalising constant. We demonstrate via numerical examples that, in some scenarios, this estimate has a much lower variance than that of the estimate obtained via the alive particle filter. The low variance of this normalising constant estimate encourages the implementation of our SMC algorithm within a particle marginal Metropolis-Hastings (PMMH) scheme, and we call the resulting methodology “alive twisted PMMH”. We numerically demonstrate, on a stochastic volatility model, how our alive twisted PMMH can converge faster than the standard alive PMMH of Jasra et al (2013).     

Bayesian Variable Selection via Particle Stochastic Search

We focus on Bayesian variable selection in regression models. One challenge is to search the huge model space adequately, while identifying high posterior probability regions. In the past decades, the main focus has been on the use of Markov chain Monte Carlo (MCMC) algorithms for these purposes. In this article, we propose a new computational approach based on sequential Monte Carlo (SMC), which we refer to as particle stochastic search (PSS). We illustrate PSS through applications to linear regression and probit models.     

Metropolis-Hastings Within Partially Collapsed Gibbs Samplers

The partially collapsed Gibbs (PCG) sampler offers a new strategy for improving the convergence of a Gibbs sampler. PCG achieves faster convergence by reducing the conditioning in some of the draws of its parent Gibbs sampler. Although this can significantly improve convergence, care must be taken to ensure that the stationary distribution is preserved. The conditional distributions sampled in a PCG sampler may be incompatible and permuting their order may upset the stationary distribution of the chain. Extra care must be taken when Metropolis-Hastings (MH) updates are used in some or all of the updates. Reducing the conditioning in an MH within Gibbs sampler can change the stationary distribution, even when the PCG sampler would work perfectly if MH were not used. In fact, a number of samplers of this sort that have been advocated in the literature do not actually have the target stationary distributions. In this article, we illustrate the challenges that may arise when using MH within a PCG sampler and develop a general strategy for using such updates while maintaining the desired stationary distribution. Theoretical arguments provide guidance when choosing between different MH within PCG sampling schemes. Finally, we illustrate the MH within PCG sampler and its computational advantage using several examples from our applied work.     

Pattern discrete and mixed Hit-and-Run for global optimization

We develop new Markov chain Monte Carlo samplers for neighborhood generation in global optimization algorithms based on Hit-and-Run. The success of Hit-and-Run as a sampler on continuous domains motivated Discrete Hit-and-Run with random biwalk for discrete domains. However, the potential for efficiencies in the implementation, which requires a randomization at each move to create the biwalk, lead us to a different approach that uses fixed patterns in generating the biwalks. We define Sphere and Box Biwalks that are pattern-based and easily implemented for discrete and mixed continuous/discrete domains. The pattern-based Hit-and-Run Markov chains preserve the convergence properties of Hit-and-Run to a target distribution. They also converge to continuous Hit-and-Run as the mesh of the discretized variables becomes finer, approaching a continuum. Moreover, we provide bounds on the finite time performance for the discrete cases of Sphere and Box Biwalks. We embed our samplers in an Improving Hit-and-Run global optimization algorithm and test their performance on a number of global optimization test problems.     

On the Convergence Complexity of Gibbs Samplers for a Family of Simple Bayesian Random Effects Models

The emergence of big data has led to so-called convergence complexity analysis, which is the study of how Markov chain Monte Carlo (MCMC) algorithms behave as the sample size, n, and/or the number of parameters, p, in the underlying data set increase. This type of analysis is often quite challenging, in part because existing results for fixed n and p are simply not sharp enough to yield good asymptotic results. One of the first convergence complexity results for an MCMC algorithm on a continuous state space is due to Yang and Rosenthal (2019), who established a mixing time result for a Gibbs sampler (for a simple Bayesian random effects model) that was introduced and studied by Rosenthal (Stat Comput 6:269–275, 1996). The asymptotic behavior of the spectral gap of this Gibbs sampler is, however, still unknown. We use a recently developed simulation technique (Qin et al. Electron J Stat 13:1790–1812, 2019) to provide substantial numerical evidence that the gap is bounded away from 0 as n → ∞. We also establish a pair of rigorous convergence complexity results for two different Gibbs samplers associated with a generalization of the random effects model considered by Rosenthal (Stat Comput 6:269–275, 1996). Our results show that, under a strong growth condition, the spectral gaps of these Gibbs samplers converge to 1 as the sample size increases.

     

Coupling stochastic EM and approximate Bayesian computation for parameter inference in state-space models

We study the class of state-space models and perform maximum likelihood estimation for the model parameters. We consider a stochastic approximation expectation–maximization (SAEM) algorithm to maximize the likelihood function with the novelty of using approximate Bayesian computation (ABC) within SAEM. The task is to provide each iteration of SAEM with a filtered state of the system, and this is achieved using an ABC sampler for the hidden state, based on sequential Monte Carlo methodology. It is shown that the resulting SAEM-ABC algorithm can be calibrated to return accurate inference, and in some situations it can outperform a version of SAEM incorporating the bootstrap filter. Two simulation studies are presented, first a nonlinear Gaussian state-space model then a state-space model having dynamics expressed by a stochastic differential equation. Comparisons with iterated filtering for maximum likelihood inference, and Gibbs sampling and particle marginal methods for Bayesian inference are presented.     

Parallel Resampling in the Particle Filter

Modern parallel computing devices, such as the graphics processing unit (GPU), have gained significant traction in scientific and statistical computing. They are particularly well-suited to data-parallel algorithms such as the particle filter, or more generally sequential Monte Carlo (SMC), which are increasingly used in statistical inference. SMC methods carry a set of weighted particles through repeated propagation, weighting, and resampling steps. The propagation and weighting steps are straightforward to parallelize, as they require only independent operations on each particle. The resampling step is more difficult, as standard schemes require a collective operation, such as a sum, across particle weights. Focusing on this resampling step, we analyze two alternative schemes that do not involve a collective operation (Metropolis and rejection resamplers), and compare them to standard schemes (multinomial, stratified, and systematic resamplers). We find that, in certain circumstances, the alternative resamplers can perform significantly faster on a GPU, and to a lesser extent on a CPU, than the standard approaches. Moreover, in single precision, the standard approaches are numerically biased for upward of hundreds of thousands of particles, while the alternatives are not. This is particularly important given greater single- than double-precision throughput on modern devices, and the consequent temptation to use single precision with a greater number of particles. Finally, we provide auxiliary functions useful for implementation, such as for the permutation of ancestry vectors to enable in-place propagation. Supplementary materials are available online.     

Convergence rate of gibbs sampler and its application

Based on the convergence rate defined by the Pearson-χ~2 distance,this pa- per discusses properties of different Gibbs sampling schemes.Under a set of regularity conditions,it is proved in this paper that the rate of convergence on systematic scan Gibbs samplers is the norm of a forward operator.We also discuss that the collapsed Gibbs sam- pler has a faster convergence rate than the systematic scan Gibbs sampler as proposed by Liu et al.Based on the definition of convergence rate of the Pearson-χ~2 distance, this paper proved this result quantitatively.According to Theorem 2,we also proved that the convergence rate defined with the spectral radius of matrix by Robert and Shau is equivalent to the corresponding radius of the forward operator.     

Quasi-Conjugate Bayes Estimates for GPD Parameters and Application to Heavy Tails Modelling

We present a quasi-conjugate Bayes approach for estimating Generalized Pareto Distribution (GPD) parameters, distribution tails and extreme quantiles within the Peaks-Over-Threshold framework. Damsleth conjugate Bayes structure on Gamma distributions is transfered to GPD. Posterior estimates are then computed by Gibbs samplers with Hastings-Metropolis steps. Accurate Bayes credibility intervals are also defined, they provide assessment of the quality of the extreme events estimates. An empirical Bayesian method is used in this work, but the suggested approach could incorporate prior information. It is shown that the obtained quasi-conjugate Bayes estimators compare well with the GPD standard estimators when simulated and real data sets are studied. AMS 2000 Subject Classification Primary—62G32, 62F15, 62G09     

An Online Expectation–Maximization Algorithm for Changepoint Models

Changepoint models are widely used to model the heterogeneity of sequential data. We present a novel sequential Monte Carlo (SMC) online expectation–maximization (EM) algorithm for estimating the static parameters of such models. The SMC online EM algorithm has a cost per time which is linear in the number of particles and could be particularly important when the data is representable as a long sequence of observations, since it drastically reduces the computational requirements for implementation. We present an asymptotic analysis for the stability of the SMC estimates used in the online EM algorithm and demonstrate the performance of this scheme by using both simulated and real data originating from DNA analysis. The supplementary materials for the article are available online.     

On reparametrization and the Gibbs sampler

Gibbs samplers derived under different parametrizations of the target density can have radically different rates of convergence. In this article, we specify conditions under which reparametrization leaves the convergence rate of a Gibbs chain unchanged. An example illustrates how these results can be exploited in convergence rate analyses.     

Implementations of the Monte Carlo EM Algorithm

The Monte Carlo EM (MCEM) algorithm is a modification of the EM algorithm where the expectation in the E-step is computed numerically through Monte Carlo simulations. The most exible and generally applicable approach to obtaining a Monte Carlo sample in each iteration of an MCEM algorithm is through Markov chain Monte Carlo (MCMC) routines such as the Gibbs and Metropolis–Hastings samplers. Although MCMC estimation presents a tractable solution to problems where the E-step is not available in closed form, two issues arise when implementing this MCEM routine: (1) how do we minimize the computational cost in obtaining an MCMC sample? and (2) how do we choose the Monte Carlo sample size? We address the first question through an application of importance sampling whereby samples drawn during previous EM iterations are recycled rather than running an MCMC sampler each MCEM iteration. The second question is addressed through an application of regenerative simulation. We obtain approximate independent and identical samples by subsampling the generated MCMC sample during different renewal periods. Standard central limit theorems may thus be used to gauge Monte Carlo error. In particular, we apply an automated rule for increasing the Monte Carlo sample size when the Monte Carlo error overwhelms the EM estimate at any given iteration. We illustrate our MCEM algorithm through analyses of two datasets fit by generalized linear mixed models. As a part of these applications, we demonstrate the improvement in computational cost and efficiency of our routine over alternative MCEM strategies.     

Optimizing random scan Gibbs samplers

The Gibbs sampler is a popular Markov chain Monte Carlo routine for generating random variates from distributions otherwise difficult to sample. A number of implementations are available for running a Gibbs sampler varying in the order through which the full conditional distributions used by the Gibbs sampler are cycled or visited. A common, and in fact the original, implementation is the random scan strategy, whereby the full conditional distributions are updated in a randomly selected order each iteration. In this paper, we introduce a random scan Gibbs sampler which adaptively updates the selection probabilities or “learns” from all previous random variates generated during the Gibbs sampling. In the process, we outline a number of variations on the random scan Gibbs sampler which allows the practitioner many choices for setting the selection probabilities and prove convergence of the induced (Markov) chain to the stationary distribution of interest. Though we emphasize flexibility in user choice and specification of these random scan algorithms, we present a minimax random scan which determines the selection probabilities through decision theoretic considerations on the precision of estimators of interest. We illustrate and apply the results presented by using the adaptive random scan Gibbs sampler developed to sample from multivariate Gaussian target distributions, to automate samplers for posterior simulation under Dirichlet process mixture models, and to fit mixtures of distributions.