On The Convergence Of Successive Substitution Sampling

The problem of finding marginal distributions of multidimensional random quantities has many applications in probability and statistics. Many of the solutions currently in use are very computationally intensive. For example, in a Bayesian inference problem with a hierarchical prior distribution, one is often driven to multidimensional numerical integration to obtain marginal posterior distributions of the model parameters of interest. Recently, however, a group of Monte Carlo integration techniques that fall under the general banner of successive substitution sampling (SSS) have proven to be powerful tools for obtaining approximate answers in a very wide variety of Bayesian modeling situations. Answers may also be obtained at low cost, both in terms of computer power and user sophistication. Important special cases of SSS include the “Gibbs sampler” described by Gelfand and Smith and the “IP algorithm” described by Tanner and Wong. The major problem plaguing users of SSS is the difficulty in ascertaining when “convergence” of the algorithm has been obtained. This problem is compounded by the fact that what is produced by the sampler is not the functional form of the desired marginal posterior distribution, but a random sample from this distribution. This article gives a general proof of the convergence of SSS and the sufficient conditions for both strong and weak convergence, as well as a convergence rate. We explore the connection between higher-order eigenfunctions of the transition operator and accelerated convergence via good initial distributions. We also provide asymptotic results for the sampling component of the error in estimating the distributions of interest. Finally, we give two detailed examples from familiar exponential family settings to illustrate the theory.