Conditioning in Markov Chain Monte Carlo

Abstract

The so-called “Rao-Blackwellized” estimators proposed by Gelfand and Smith do not always reduce variance in Markov chain Monte Carlo when the dependence in the Markov chain is taken into account. An illustrative example is given, and a theorem characterizing the necessary and sufficient condition for such an estimator to always reduce variance is proved.     

Bayesian Monte Carlo estimation for profile hidden Markov models

Hidden Markov models are used as tools for pattern recognition in a number of areas, ranging from speech processing to biological sequence analysis. Profile hidden Markov models represent a class of so-called “left–right” models that have an architecture that is specifically relevant to classification of proteins into structural families based on their amino acid sequences. Standard learning methods for such models employ a variety of heuristics applied to the expectation-maximization implementation of the maximum likelihood estimation procedure in order to find the global maximum of the likelihood function. Here, we compare maximum likelihood estimation to fully Bayesian estimation of parameters for profile hidden Markov models with a small number of parameters. We find that, relative to maximum likelihood methods, Bayesian methods assign higher scores to data sequences that are distantly related to the pattern consensus, show better performance in classifying these sequences correctly, and continue to perform robustly with regard to misspecification of the number of model parameters. Though our study is limited in scope, we expect our results to remain relevant for models with a large number of parameters and other types of left–right hidden Markov models.     

Importance sampling correction versus standard averages of reversible MCMCs in terms of the asymptotic variance

We establish an ordering criterion for the asymptotic variances of two consistent Markov chain Monte Carlo (MCMC) estimators: an importance sampling (IS) estimator, based on an approximate reversible chain and subsequent IS weighting, and a standard MCMC estimator, based on an exact reversible chain. Essentially, we relax the criterion of the Peskun type covariance ordering by considering two different invariant probabilities, and obtain, in place of a strict ordering of asymptotic variances, a bound of the asymptotic variance of IS by that of the direct MCMC. Simple examples show that IS can have arbitrarily better or worse asymptotic variance than Metropolis–Hastings and delayed-acceptance (DA) MCMC. Our ordering implies that IS is guaranteed to be competitive up to a factor depending on the supremum of the (marginal) IS weight. We elaborate upon the criterion in case of unbiased estimators as part of an auxiliary variable framework. We show how the criterion implies asymptotic variance guarantees for IS in terms of pseudo-marginal (PM) and DA corrections, essentially if the ratio of exact and approximate likelihoods is bounded. We also show that convergence of the IS chain can be less affected by unbounded high-variance unbiased estimators than PM and DA chains.     

A Practical Sequential Stopping Rule for High-Dimensional Markov Chain Monte Carlo

A current challenge for many Bayesian analyses is determining when to terminate high-dimensional Markov chain Monte Carlo simulations. To this end, we propose using an automated sequential stopping procedure that terminates the simulation when the computational uncertainty is small relative to the posterior uncertainty. Further, we show this stopping rule is equivalent to stopping when the effective sample size is sufficiently large. Such a stopping rule has previously been shown to work well in settings with posteriors of moderate dimension. In this article, we illustrate its utility in high-dimensional simulations while overcoming some current computational issues. As examples, we consider two complex Bayesian analyses on spatially and temporally correlated datasets. The first involves a dynamic space-time model on weather station data and the second a spatial variable selection model on fMRI brain imaging data. Our results show the sequential stopping rule is easy to implement, provides uncertainty estimates, and performs well in high-dimensional settings. Supplementary materials for this article are available online.     

Monte carlo estimation of the sampling distribution of nonlinear model parameter estimators

The sampling distribution of parameter estimators can be summarized by moments, fractiles or quantiles. For nonlinear models, these quantities are often approximated by power series, approximated by transformed systems, or estimated by Monte Carlo sampling. A control variate approach based on a linear approximation of the nonlinear model is introduced here to reduce the Monte Carlo sampling necessary to achieve a given accuracy. The particular linear approximation chosen has several advantages: its moments and other properties are known, it is easy to implement, and there is a correspondence to asymptotic results that permits assessment of control variate effectiveness prior to sampling via measures of nonlinearity. Empirical results for several nonlinear problems are presented.This research was supported in part by the Office of Naval Research under Contract N00014-79-C-0832.     

On Markov Chain Monte Carlo Acceleration

Abstract

Markov chain Monte Carlo (MCMC) methods are currently enjoying a surge of interest within the statistical community. The goal of this work is to formalize and support two distinct adaptive strategies that typically accelerate the convergence of an MCMC algorithm. One approach is through resampling; the other incorporates adaptive switching of the transition kernel. Support is both by analytic arguments and simulation study. Application is envisioned in low-dimensional but nontrivial problems. Two pathological illustrations are presented. Connections with reparameterization are discussed as well as possible difficulties with infinitely often adaptation.     

Partition-Weighted Monte Carlo Estimation

Although various efficient and sophisticated Markov chain Monte Carlo sampling methods have been developed during the last decade, the sample mean is still a dominant in computing Bayesian posterior quantities. The sample mean is simple, but may not be efficient. The weighted sample mean is a natural generalization of the sample mean. In this paper, a new weighted sample mean is proposed by partitioning the support of posterior distribution, so that the same weight is assigned to observations that belong to the same subset in the partition. A novel application of this new weighted sample mean in computing ratios of normalizing constants and necessary theory are provided. Illustrative examples are given to demonstrate the methodology.     

Discrete stochastic consistent estimators of the Monte Carlo method

The efficiency of discrete stochastic consistent estimators (the weighted uniform sampling and estimator with a correcting multiplier) of the Monte Carlo method is investigated. Confidence intervals and upper bounds on the variances are obtained, and the computational cost of the corresponding discrete stochastic numerical scheme is estimated.     

The Alive Particle Filter and Its Use in Particle Markov Chain Monte Carlo

In the following article, we investigate a particle filter for approximating Feynman–Kac models with indicator potentials and we use this algorithm within Markov chain Monte Carlo (MCMC) to learn static parameters of the model. Examples of such models include approximate Bayesian computation (ABC) posteriors associated with hidden Markov models (HMMs) or rare-event problems. Such models require the use of advanced particle filter or MCMC algorithms to perform estimation. One of the drawbacks of existing particle filters is that they may “collapse,” in that the algorithm may terminate early, due to the indicator potentials. In this article, using a newly developed special case of the locally adaptive particle filter, we use an algorithm that can deal with this latter problem, while introducing a random cost per-time step. In particular, we show how this algorithm can be used within MCMC, using particle MCMC. It is established that, when not taking into account computational time, when the new MCMC algorithm is applied to a simplified model it has a lower asymptotic variance in comparison to a standard particle MCMC algorithm. Numerical examples are presented for ABC approximations of HMMs.     

基于贝叶斯MCMC算法的美式期权定价

鉴于美式期权的定价具有后向迭代搜索特征,本文结合Longstaff和Schwartz提出的美式期权定价的最小二乘模拟方法,研究基于马尔科夫链蒙特卡洛算法对回归方程系数的估计,实现对美式期权的双重模拟定价.通过对无红利美式看跌股票期权定价进行大量实证模拟,从期权价值定价误差等方面同著名的最小二乘蒙特卡洛模拟方法进行对比分析,结果表明基于MCMC回归算法给出的美式期权定价具有更高的精确度.模拟实证结果表明本文提出的对美式期权定价方法具有较好的可行性、有效性与广泛的适用性.该方法的不足之处就是类似于一般的蒙特卡洛方法,会使得求解的计算量有所加大.     

Markov Chain Monte Carlo Convergence Assessment via Two-Way Analysis of Variance

Abstract

In this article we discuss the problem of assessing the performance of Markov chain Monte Carlo (MCMC) algorithms on the basis of simulation output. In essence, we extend the original ideas of Gelman and Rubin and, more recently, Brooks and Gelman, to problems where we are able to split the variation inherent within the MCMC simulation output into two distinct groups. We show how such a diagnostic may be useful in assessing the performance of MCMC samplers addressing model choice problems, such as the reversible jump MCMC algorithm. In the model choice context, we show how the reversible jump MCMC simulation output for parameters that retain a coherent interpretation throughout the simulation, can be used to assess convergence. By considering various decompositions of the sampling variance of this parameter, we can assess the performance of our MCMC sampler in terms of its mixing properties both within and between models and we illustrate our approach in both the graphical Gaussian models and normal mixtures context. Finally, we provide an example of the application of our diagnostic to the assessment of the influence of different starting values on MCMC simulation output, thereby illustrating the wider utility of our method beyond the Bayesian model choice and reversible jump MCMC context.     

Resampling Markov Chain Monte Carlo Algorithms: Basic Analysis and Empirical Comparisons

Sampling from complex distributions is an important but challenging topic in scientific and statistical computation. We synthesize three ideas, tempering, resampling, and Markov moving, and propose a general framework of resampling Markov chain Monte Carlo (MCMC). This framework not only accommodates various existing algorithms, including resample-move, importance resampling MCMC, and equi-energy sampling, but also leads to a generalized resample-move algorithm. We provide some basic analysis of these algorithms within the general framework, and present three simulation studies to compare these algorithms together with parallel tempering in the difficult situation where new modes emerge in the tails of previous tempering distributions. Our analysis and empirical results suggest that generalized resample-move tends to perform the best among all the algorithms studied when the Markov kernels lead to fast mixing or even locally so toward restricted distributions, whereas parallel tempering tends to perform the best when the Markov kernels lead to slow mixing, without even converging fast to restricted distributions. Moreover, importance resampling MCMC and equi-energy sampling perform similarly to each other, often worse than independence Metropolis resampling MCMC. Therefore, different algorithms seem to have advantages in different settings.     

Robustness and applicability of Markov chain Monte Carlo algorithms for eigenvalue problems

In this paper we analyse applicability and robustness of Markov chain Monte Carlo algorithms for eigenvalue problems. We restrict our consideration to real symmetric matrices.

Almost Optimal Monte Carlo (MAO) algorithms for solving eigenvalue problems are formulated. Results for the structure of both – systematic and probability error are presented. It is shown that the values of both errors can be controlled independently by different algorithmic parameters. The results present how the systematic error depends on the matrix spectrum. The analysis of the probability error is presented. It shows that the close (in some sense) the matrix under consideration is to the stochastic matrix the smaller is this error. Sufficient conditions for constructing robust and interpolation Monte Carlo algorithms are obtained. For stochastic matrices an interpolation Monte Carlo algorithm is constructed.

A number of numerical tests for large symmetric dense matrices are performed in order to study experimentally the dependence of the systematic error from the structure of matrix spectrum. We also study how the probability error depends on the balancing of the matrix.     

User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient

In this paper, we study the problem of sampling from a given probability density function that is known to be smooth and strongly log-concave. We analyze several methods of approximate sampling based on discretizations of the (highly overdamped) Langevin diffusion and establish guarantees on its error measured in the Wasserstein-2 distance. Our guarantees improve or extend the state-of-the-art results in three directions. First, we provide an upper bound on the error of the first-order Langevin Monte Carlo (LMC) algorithm with optimized varying step-size. This result has the advantage of being horizon free (we do not need to know in advance the target precision) and to improve by a logarithmic factor the corresponding result for the constant step-size. Second, we study the case where accurate evaluations of the gradient of the log-density are unavailable, but one can have access to approximations of the aforementioned gradient. In such a situation, we consider both deterministic and stochastic approximations of the gradient and provide an upper bound on the sampling error of the first-order LMC that quantifies the impact of the gradient evaluation inaccuracies. Third, we establish upper bounds for two versions of the second-order LMC, which leverage the Hessian of the log-density. We provide non asymptotic guarantees on the sampling error of these second-order LMCs. These guarantees reveal that the second-order LMC algorithms improve on the first-order LMC in ill-conditioned settings.     

Bayesian and non-bayesian analysis of gamma stochastic frontier models by Markov Chain Monte Carlo methods

Summary  This paper considers simulation-based approaches for the gamma stochastic frontier model. Efficient Markov chain Monte Carlo methods are proposed for sampling the posterior distribution of the parameters. Maximum likelihood estimation is also discussed based on the stochastic approximation algorithm. The methods are applied to a data set of the U.S. electric utility industry. The authors are grateful to two anonymous referees for their useful comments, which improved an earlier version of the paper. The first author also thanks the financial support by the Japanese Ministry of Education, Culture, Sports, Science and Technology under the Grant-in-Aid for Scientific Research No.14730022.     

Bayesian estimation of NIG models via Markov chain Monte Carlo methods

The normal inverse Gaussian (NIG) distribution is a promising alternative for modelling financial data since it is a continuous distribution that allows for skewness and fat tails. There is an increasing number of applications of the NIG distribution to financial problems. Due to the complicated nature of its density, estimation procedures are not simple. In this paper we propose Bayesian estimation for the parameters of the NIG distribution via an MCMC scheme based on the Gibbs sampler. Our approach makes use of the data augmentation provided by the mixture representation of the distribution. We also extend the model to allow for modelling heteroscedastic regression situations. Examples with financial and simulated data are provided. Copyright © 2004 John Wiley & Sons, Ltd.     

Smoothed Monte Carlo estimators for the time-in-the-red in risk processes

We consider a modified version of the de Finetti model in insurance risk theory in which, when surpluses become negative the company has the possibility of borrowing, and thus continue its operation. For this model we examine the problem of estimating the time-in-the red over a finite horizon via simulation. We propose a smoothed estimator based on a conditioning argument which is very simple to implement as well as particularly efficient, especially when the claim distribution is heavy tailed. We establish unbiasedness for this estimator and show that its variance is lower than the naïve estimator based on counts. Finally we present a number of simulation results showing that the smoothed estimator has variance which is often significantly lower than that of the naïve Monte-Carlo estimator.     

Explicit error bounds for lazy reversible Markov chain Monte Carlo

We prove explicit, i.e., non-asymptotic, error bounds for Markov Chain Monte Carlo methods, such as the Metropolis algorithm. The problem is to compute the expectation (or integral) of f$f$ with respect to a measure π$\pi$ which can be given by a density ?$?$ with respect to another measure. A straight simulation of the desired distribution by a random number generator is in general not possible. Thus it is reasonable to use Markov chain sampling with a burn-in. We study such an algorithm and extend the analysis of Lovasz and Simonovits [L. Lovász, M. Simonovits, Random walks in a convex body and an improved volume algorithm, Random Structures Algorithms 4 (4) (1993) 359–412] to obtain an explicit error bound.     

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.