Inference from Iterative Simulation Using Multiple Sequences

TLDR

Iterative simulation methods such as the Gibbs sampler can efficiently summarize multivariate distributions, yet naïve use may produce misleading results, motivating the need for reliable inference techniques for Bayesian posteriors. The study aims to provide simple, broadly applicable methods for researchers focused on scientific inference from iterative simulation outputs rather than on the underlying probability theory. The authors recommend running multiple independent, overdispersed chains, collecting at each iteration distributional estimates and convergence diagnostics, then applying normal‑theory approximations to infer Bayesian posteriors, demonstrated on a random‑effects mixture model of reaction times.

Abstract

The Gibbs sampler, the algorithm of Metropolis and similar iterative simulation methods are potentially very helpful for summarizing multivariate distributions. Used naively, however, iterative simulation can give misleading answers. Our methods are simple and generally applicable to the output of any iterative simulation; they are designed for researchers primarily interested in the science underlying the data and models they are analyzing, rather than for researchers interested in the probability theory underlying the iterative simulations themselves. Our recommended strategy is to use several independent sequences, with starting points sampled from an overdispersed distribution. At each step of the iterative simulation, we obtain, for each univariate estimand of interest, a distributional estimate and an estimate of how much sharper the distributional estimate might become if the simulations were continued indefinitely. Because our focus is on applied inference for Bayesian posterior distributions in real problems, which often tend toward normality after transformations and marginalization, we derive our results as normal-theory approximations to exact Bayesian inference, conditional on the observed simulations. The methods are illustrated on a random-effects mixture model applied to experimental measurements of reaction times of normal and schizophrenic patients.