Bayes and MCMC for Undergraduates

Students of statistics should be taught the ideas and methods that are widely used in practice and that will help them understand the world of statistics. Today, this means teaching them about Bayesian methods. In this article, I present ideas on teaching an undergraduate Bayesian course that uses Markov chain Monte Carlo and that can be a second course or, for strong students, a first course in statistics.     

Regenerative Markov Chain Monte Carlo for Any Distribution

While Markov chain Monte Carlo (MCMC) methods are frequently used for difficult calculations in a wide range of scientific disciplines, they suffer from a serious limitation: their samples are not independent and identically distributed. Consequently, estimates of expectations are biased if the initial value of the chain is not drawn from the target distribution. Regenerative simulation provides an elegant solution to this problem. In this article, we propose a simple regenerative MCMC algorithm to generate variates for any distribution.     

Markov Chains and De-initializing Processes

We define a notion of de-initializing Markov chains. We prove that to analyse convergence of Markov chains to stationarity, it suffices to analyse convergence of a de-initializing chain. Applications are given to Markov chain Monte Carlo algorithms and to convergence diagnostics.     

Bounding convergence rates for Markov chains: An example of the use of computer algebra

Kolassa and Tanner (J. Am. Stat. Assoc. (1994) 89, 697–702) present the Gibbs-Skovgaard algorithm for approximate conditional inference. Kolassa (Ann Statist. (1999), 27, 129–142) gives conditions under which their Markov chain is known to converge. This paper calculates explicity bounds on convergence rates in terms calculable directly from chain transition operators. These results are useful in cases like those considered by Kolassa (1999).     

LONG-RANGE DEPENDENT COMMON FACTOR MODELS: A BAYESIAN APPROACH

We propose a simulation-based Bayesian approach to analyze multivariate time series with possible common long-range dependent factors. A state-space approach is used to represent the likelihood function in a tractable manner. The approach taken here allows for extension to fit a non-Gaussian multivariate stochastic volatility (MVSV) model with common long-range dependent components. The method is illustrated for a set of stock returns for companies having similar annual sales.     

On the Value of derivative evaluations and random walk suppression in Markov Chain Monte Carlo algorithms

Two strategies that can potentially improve Markov Chain Monte Carlo algorithms are to use derivative evaluations of the target density, and to suppress random walk behaviour in the chain. The use of one or both of these strategies has been investigated in a few specific applications, but neither is used routinely. We undertake a broader evaluation of these techniques, with a view to assessing their utility for routine use. In addition to comparing different algorithms, we also compare two different ways in which the algorithms can be applied to a multivariate target distribution. Specifically, the univariate version of an algorithm can be applied repeatedly to one-dimensional conditional distributions, or the multivariate version can be applied directly to the target distribution.     

Bayesian inference in hidden Markov models through the reversible jump Markov chain Monte Carlo method

Hidden Markov models form an extension of mixture models which provides a flexible class of models exhibiting dependence and a possibly large degree of variability. We show how reversible jump Markov chain Monte Carlo techniques can be used to estimate the parameters as well as the number of components of a hidden Markov model in a Bayesian framework. We employ a mixture of zero-mean normal distributions as our main example and apply this model to three sets of data from finance, meteorology and geomagnetism.     

Mixed Marginal Copula Modeling

ABSTRACT

This article extends the literature on copulas with discrete or continuous marginals to the case where some of the marginals are a mixture of discrete and continuous components. We do so by carefully defining the likelihood as the density of the observations with respect to a mixed measure. The treatment is quite general, although we focus on mixtures of Gaussian and Archimedean copulas. The inference is Bayesian with the estimation carried out by Markov chain Monte Carlo. We illustrate the methodology and algorithms by applying them to estimate a multivariate income dynamics model. Supplementary materials for this article are available online.     

Approximating hidden Gaussian Markov random fields

Summary.  Gaussian Markov random-field (GMRF) models are frequently used in a wide variety of applications. In most cases parts of the GMRF are observed through mutually independent data; hence the full conditional of the GMRF, a hidden GMRF (HGMRF), is of interest. We are concerned with the case where the likelihood is non-Gaussian, leading to non-Gaussian HGMRF models. Several researchers have constructed block sampling Markov chain Monte Carlo schemes based on approximations of the HGMRF by a GMRF, using a second-order expansion of the log-density at or near the mode. This is possible as the GMRF approximation can be sampled exactly with a known normalizing constant. The Markov property of the GMRF approximation yields computational efficiency.The main contribution in the paper is to go beyond the GMRF approximation and to construct a class of non-Gaussian approximations which adapt automatically to the particular HGMRF that is under study. The accuracy can be tuned by intuitive parameters to nearly any precision. These non-Gaussian approximations share the same computational complexity as those which are based on GMRFs and can be sampled exactly with computable normalizing constants. We apply our approximations in spatial disease mapping and model-based geostatistical models with different likelihoods, obtain procedures for block updating and construct Metropolized independence samplers.     

Problem Solving is Often a Matter of Cooking Up an Appropriate Markov Chain*

Abstract. By means of a series of examples, from classic contributions to probability theory as well as the author's own, an attempt is made to convince the reader that problem solving is often a matter of cooking up an appropriate Markov chain. Topics touched upon along the way include coupling, correlation inequalities, and percolation.     

Markov Chain Monte Carlo Methods for Fitting Spatiotemporal Stochastic Models in Plant Epidemiology

Strategies for controlling plant epidemics are investigated by fitting continuous time spatiotemporal stochastic models to data consisting of maps of disease incidence observed at discrete times. Markov chain Monte Carlo methods are used for fitting two such models to data describing the spread of citrus tristeza virus (CTV) in an orchard. The approach overcomes some of the difficulties encountered when fitting stochastic models to infrequent observations of a continuous process. The results of the analysis cast doubt on the effectiveness of a strategy identified from a previous spatial analysis of the CTV data. Extensions of the approaches to more general models and other problems are also considered.     

A distance-based diagnostic for trans-dimensional Markov chains

Over the last decade the use of trans-dimensional sampling algorithms has become endemic in the statistical literature. In spite of their application however, there are few reliable methods to assess whether the underlying Markov chains have reached their stationary distribution. In this article we present a distance-based method for the comparison of trans-dimensional Markov chain sample output for a broad class of models. This diagnostic will simultaneously assess deviations between and within chains. Illustration of the analysis of Markov chain sample-paths is presented in simulated examples and in two common modelling situations: a finite mixture analysis and a change-point problem.     

Minimal Basis for a Connected Markov Chain over 3 × 3 ×K Contingency Tables with Fixed Two-Dimensional Marginals

This paper considers a connected Markov chain for sampling 3 × 3 ×K contingency tables having fixed two‐dimensional marginal totals. Such sampling arises in performing various tests of the hypothesis of no three‐factor interactions. A Markov chain algorithm is a valuable tool for evaluating P‐values, especially for sparse datasets where large‐sample theory does not work well. To construct a connected Markov chain over high‐dimensional contingency tables with fixed marginals, algebraic algorithms have been proposed. These algorithms involve computations in polynomial rings using Gröbner bases. However, algorithms based on Gröbner bases do not incorporate symmetry among variables and are very time‐consuming when the contingency tables are large. We construct a minimal basis for a connected Markov chain over 3 × 3 ×K contingency tables. The minimal basis is unique. Some numerical examples illustrate the practicality of our algorithms.     

Importance sampling type estimators based on approximate marginal Markov chain Monte Carlo

We consider importance sampling (IS) type weighted estimators based on Markov chain Monte Carlo (MCMC) targeting an approximate marginal of the target distribution. In the context of Bayesian latent variable models, the MCMC typically operates on the hyperparameters, and the subsequent weighting may be based on IS or sequential Monte Carlo (SMC), but allows for multilevel techniques as well. The IS approach provides a natural alternative to delayed acceptance (DA) pseudo-marginal/particle MCMC, and has many advantages over DA, including a straightforward parallelization and additional flexibility in MCMC implementation. We detail minimal conditions which ensure strong consistency of the suggested estimators, and provide central limit theorems with expressions for asymptotic variances. We demonstrate how our method can make use of SMC in the state space models context, using Laplace approximations and time-discretized diffusions. Our experimental results are promising and show that the IS-type approach can provide substantial gains relative to an analogous DA scheme, and is often competitive even without parallelization.     

Statistical inference for discretely observed Markov jump processes

Summary.  Likelihood inference for discretely observed Markov jump processes with finite state space is investigated. The existence and uniqueness of the maximum likelihood estimator of the intensity matrix are investigated. This topic is closely related to the imbedding problem for Markov chains. It is demonstrated that the maximum likelihood estimator can be found either by the EM algorithm or by a Markov chain Monte Carlo procedure. When the maximum likelihood estimator does not exist, an estimator can be obtained by using a penalized likelihood function or by the Markov chain Monte Carlo procedure with a suitable prior. The methodology and its implementation are illustrated by examples and simulation studies.     

Markov chain Monte Carlo tests for designed experiments

We consider conditional exact tests of factor effects in designed experiments for discrete response variables. Similarly to the analysis of contingency tables, a Markov chain Monte Carlo method can be used for performing exact tests, when large-sample approximations are poor and the enumeration of the conditional sample space is infeasible. For designed experiments with a single observation for each run, we formulate log-linear or logistic models and consider a connected Markov chain over an appropriate sample space. In particular, we investigate fractional factorial designs with 2^p-q$2^{p-q}$ runs, noting correspondences to the models for 2^p-q$2^{p-q}$ contingency tables.     

Markov chain Monte Carlo methods for high dimensional inversion in remote sensing

Summary.  We discuss the inversion of the gas profiles (ozone, NO₃, NO₂, aerosols and neutral density) in the upper atmosphere from the spectral occultation measurements. The data are produced by the 'Global ozone monitoring of occultation of stars' instrument on board the Envisat satellite that was launched in March 2002. The instrument measures the attenuation of light spectra at various horizontal paths from about 100 km down to 10–20 km. The new feature is that these data allow the inversion of the gas concentration height profiles. A short introduction is given to the present operational data management procedure with examples of the first real data inversion. Several solution options for a more comprehensive statistical inversion are presented. A direct inversion leads to a non-linear model with hundreds of parameters to be estimated. The problem is solved with an adaptive single-step Markov chain Monte Carlo algorithm. Another approach is to divide the problem into several non-linear smaller dimensional problems, to run parallel adaptive Markov chain Monte Carlo chains for them and to solve the gas profiles in repetitive linear steps. The effect of grid size is discussed, and we present how the prior regularization takes the grid size into account in a way that effectively leads to a grid-independent inversion.     

Fast sampling of Gaussian Markov random fields

This paper demonstrates how Gaussian Markov random fields (conditional autoregressions) can be sampled quickly by using numerical techniques for sparse matrices. The algorithm is general and efficient, and expands easily to various forms for conditional simulation and evaluation of normalization constants. We demonstrate its use by constructing efficient block updates in Markov chain Monte Carlo algorithms for disease mapping.     

Convergence of Slice Sampler Markov Chains

We analyse theoretical properties of the slice sampler. We find that the algorithm has extremely robust geometric ergodicity properties. For the case of just one auxiliary variable, we demonstrate that the algorithm is stochastically monotone, and we deduce analytic bounds on the total variation distance from stationarity of the method by using Foster–Lyapunov drift condition methodology.     

Bayesian Inference and Forecasting in Dynamic Neural Networks with Fully Markov Switching ARCH Noises

We deal with one-layer feed-forward neural network for the Bayesian analysis of nonlinear time series. Noises are modeled nonlinearly and nonnormally, by means of ARCH models whose parameters are all dependent on a hidden Markov chain. Parameter estimation is performed by sampling from the posterior distribution via Evolutionary Monte Carlo algorithm, in which two new crossover operators have been introduced. Unknown parameters of the model also include the missing values which can occur within the observed series, so, considering future values as missing, it is also possible to compute point and interval multi-step-ahead predictions.