Dynamic Markov Bases

This article presents a computational approach for generating Markov bases for multiway contingency tables whose cell counts might be constrained by fixed marginals and by lower and upper bounds. Our framework includes tables with structural zeros as a particular case. Instead of computing the entire Markov bases in an initial step, our framework finds sets of local moves that connect each table in the reference set with a set of neighbor tables. We construct a Markov chain on the reference set of tables that requires only a set of local moves at each iteration. The union of these sets of local moves forms a dynamic Markov basis. We illustrate the practicality of our algorithms in the estimation of exact p-values for a three-way table with structural zeros and a sparse eight-way table. This article has online supplementary materials.     

Minimal invariant Markov basis for sampling contingency tables with fixed marginals

In this paper we define an invariant Markov basis for a connected Markov chain over the set of contingency tables with fixed marginals and derive some characterizations of minimality of the invariant basis. We also give a necessary and sufficient condition for uniqueness of minimal invariant Markov bases. By considering the invariance, Markov bases can be presented very concisely. As an example, we present minimal invariant Markov bases for all 2 × 2 × 2 × 2 hierarchical models. The invariance here refers to permutation of indices of each axis of the contingency tables. If the categories of each axis do not have any order relations among them, it is natural to consider the action of the symmetric group on each axis of the contingency table. A general algebraic algorithm for obtaining a Markov basis was given by Diaconis and Sturmfels (The Annals of Statistics, 26, 363–397, 1998). Their algorithm is based on computing Gröbner basis of a well-specified polynomial ideal. However, the reduced Gröbner basis depends on the particular term order and is not symmetric. Therefore, it is of interest to consider the properties of invariant Markov basis.     

Nonparametric Bayesian Modeling for Multivariate Ordinal Data

This article proposes a probability model for k-dimensional ordinal outcomes, that is, it considers inference for data recorded in k-dimensional contingency tables with ordinal factors. The proposed approach is based on full posterior inference, assuming a flexible underlying prior probability model for the contingency table cell probabilities. We use a variation of the traditional multivariate probit model, with latent scores that determine the observed data. In our model, a mixture of normals prior replaces the usual single multivariate normal model for the latent variables. By augmenting the prior model to a mixture of normals we generalize inference in two important ways. First, we allow for varying local dependence structure across the contingency table. Second, inference in ordinal multivariate probit models is plagued by problems related to the choice and resampling of cutoffs defined for these latent variables. We show how the proposed mixture model approach entirely removes these problems. We illustrate the methodology with two examples, one simulated dataset and one dataset of interrater agreement.     

A dynamic programming approach to efficient sampling from Boltzmann distributions

Markov chain methods for Boltzmann sampling work in phases with decreasing temperatures. The number of transitions in each phase crucially affects terminal state distribution. We employ dynamic programming to allocate iterations to phases to improve guarantees on sample quality. Numerical experiments on the Ising model are presented.     

A markov chain sampler for contingency table exact inference

Summary  In the inference of contingency table, when the cell counts are not large enough for asymptotic approximation, conditioning exact method is used and often computationally impractical for large tables. Instead, various sampling methods can be used. Based on permutation, the Monte Carlo sampling may become again impractical for large tables. For this, existing the Markov chain method is to sample a few elements of the table at each iteration and is inefficient. Here we consider a Markov chain, in which a sub-table of user specified size is updated at each iteration, and it achieves high sampling efficiency. Some theoretical properties of the chain and its applications to some commonly used tables are discussed. As an illustration, this method is applied to the exact test of the Hardy-Weinberg equilibrium in the population genetics context.     

Statistically Efficient Thinning of a Markov Chain Sampler

It is common to subsample Markov chain output to reduce the storage burden. Geyer shows that discarding k ? 1 out of every k observations will not improve statistical efficiency, as quantified through variance in a given computational budget. That observation is often taken to mean that thinning Markov chain Monte Carlo (MCMC) output cannot improve statistical efficiency. Here, we suppose that it costs one unit of time to advance a Markov chain and then θ > 0 units of time to compute a sampled quantity of interest. For a thinned process, that cost θ is incurred less often, so it can be advanced through more stages. Here, we provide examples to show that thinning will improve statistical efficiency if θ is large and the sample autocorrelations decay slowly enough. If the lag ? ? 1 autocorrelations of a scalar measurement satisfy ρ_? > ρ_{? + 1} > 0, then there is always a θ < ∞ at which thinning becomes more efficient for averages of that scalar. Many sample autocorrelation functions resemble first order AR(1) processes with ρ_? = ρ^|?| for some ? 1 < ρ < 1. For an AR(1) process, it is possible to compute the most efficient subsampling frequency k. The optimal k grows rapidly as ρ increases toward 1. The resulting efficiency gain depends primarily on θ, not ρ. Taking k = 1 (no thinning) is optimal when ρ ? 0. For ρ > 0, it is optimal if and only if θ ? (1 ? ρ)²/(2ρ). This efficiency gain never exceeds 1 + θ. This article also gives efficiency bounds for autocorrelations bounded between those of two AR(1) processes. Supplementary materials for this article are available online.     

A finiteness theorem for Markov bases of hierarchical models

We show that the complexity of the Markov bases of multidimensional tables stabilizes eventually if a single table dimension is allowed to vary. In particular, if this table dimension is greater than a computable bound, the Markov bases consist of elements from Markov bases of smaller tables. We give an explicit formula for this bound in terms of Graver bases. We also compute these Markov and Graver complexities for all K×2×2×2 tables.     

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 Analysis of a Two-State Markov Modulated Poisson Process

Abstract

We postulate observations from a Poisson process whose rate parameter modulates between two values determined by an unobserved Markov chain. The theory switches from continuous to discrete time by considering the intervals between observations as a sequence of dependent random variables. A result from hidden Markov models allows us to sample from the posterior distribution of the model parameters given the observed event times using a Gibbs sampler with only two steps per iteration.     

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.     

On solving integral equations using Markov chain Monte Carlo methods

In this paper, we propose an original approach to the solution of Fredholm equations of the second kind. We interpret the standard Von Neumann expansion of the solution as an expectation with respect to a probability distribution defined on a union of subspaces of variable dimension. Based on this representation, it is possible to use trans-dimensional Markov chain Monte Carlo (MCMC) methods such as Reversible Jump MCMC to approximate the solution numerically. This can be an attractive alternative to standard Sequential Importance Sampling (SIS) methods routinely used in this context. To motivate our approach, we sketch an application to value function estimation for a Markov decision process. Two computational examples are also provided.     

Pattern Hit-and-Run for sampling efficiently on polytopes

Pattern Hit-and-Run (PHR) is a Markov chain Monte Carlo sampler for a target distribution that was originally designed for general sets embedded in a box. A specific set of interest to many applications is a polytope intersected with discrete or mixed continuous/discrete lattices. PHR requires an acceptance/rejection mechanism along a bidirectional walk to guarantee feasibility. We remove this inefficiency by utilizing the linearity of the constraints defining the polytope, so each iteration of PHR can be efficiently implemented even though the variables are allowed to be integer valued. Moreover, PHR converges to a uniform distribution in polynomial time for a class of discrete polytopes.     

Torpid mixing of the Wang–Swendsen–Kotecký algorithm for sampling colorings

We study the problem of sampling uniformly at random from the set of k-colorings of a graph with maximum degree Δ. We focus attention on the Markov chain Monte Carlo method, particularly on a popular Markov chain for this problem, the Wang–Swendsen–Kotecký (WSK) algorithm. The second author recently proved that the WSK algorithm quickly converges to the desired distribution when k11Δ/6. We study how far these positive results can be extended in general. In this note we prove the first non-trivial results on when the WSK algorithm takes exponentially long to reach the stationary distribution and is thus called torpidly mixing. In particular, we show that the WSK algorithm is torpidly mixing on a family of bipartite graphs when 3k<Δ/(20logΔ), and on a family of planar graphs for any number of colors. We also give a family of graphs for which, despite their small chromatic number, the WSK algorithm is not ergodic when kΔ/2, provided k is larger than some absolute constant k₀.     

Nearest neighbor conditional estimation for Harris recurrent Markov chains

This paper is concerned with consistent nearest neighbor time series estimation for data generated by a Harris recurrent Markov chain on a general state space. It is shown that nearest neighbor estimation is consistent in this general time series context, using simple and weak conditions. The results proved here, establish consistency, in a unified manner, for a large variety of problems, e.g. autoregression function estimation, and, more generally, extremum estimators as well as sequential forecasting. Finally, under additional conditions, it is also shown that the estimators are asymptotically normal.     

A Novel Approach for Markov Random Field With Intractable Normalizing Constant on Large Lattices

The pseudo likelihood method of Besag (1974) has remained a popular method for estimating Markov random field on a very large lattice, despite various documented deficiencies. This is partly because it remains the only computationally tractable method for large lattices. We introduce a novel method to estimate Markov random fields defined on a regular lattice. The method takes advantage of conditional independence structures and recursively decomposes a large lattice into smaller sublattices. An approximation is made at each decomposition. Doing so completely avoids the need to compute the troublesome normalizing constant. The computational complexity is O(N), where N is the number of pixels in the lattice, making it computationally attractive for very large lattices. We show through simulations, that the proposed method performs well, even when compared with methods using exact likelihoods. Supplementary material for this article is available online.     

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.     

Honest Importance Sampling With Multiple Markov Chains

Importance sampling is a classical Monte Carlo technique in which a random sample from one probability density, π₁, is used to estimate an expectation with respect to another, π. The importance sampling estimator is strongly consistent and, as long as two simple moment conditions are satisfied, it obeys a central limit theorem (CLT). Moreover, there is a simple consistent estimator for the asymptotic variance in the CLT, which makes for routine computation of standard errors. Importance sampling can also be used in the Markov chain Monte Carlo (MCMC) context. Indeed, if the random sample from π₁ is replaced by a Harris ergodic Markov chain with invariant density π₁, then the resulting estimator remains strongly consistent. There is a price to be paid, however, as the computation of standard errors becomes more complicated. First, the two simple moment conditions that guarantee a CLT in the iid case are not enough in the MCMC context. Second, even when a CLT does hold, the asymptotic variance has a complex form and is difficult to estimate consistently. In this article, we explain how to use regenerative simulation to overcome these problems. Actually, we consider a more general setup, where we assume that Markov chain samples from several probability densities, π₁, …, π_k, are available. We construct multiple-chain importance sampling estimators for which we obtain a CLT based on regeneration. We show that if the Markov chains converge to their respective target distributions at a geometric rate, then under moment conditions similar to those required in the iid case, the MCMC-based importance sampling estimator obeys a CLT. Furthermore, because the CLT is based on a regenerative process, there is a simple consistent estimator of the asymptotic variance. We illustrate the method with two applications in Bayesian sensitivity analysis. The first concerns one-way random effect models under different priors. The second involves Bayesian variable selection in linear regression, and for this application, importance sampling based on multiple chains enables an empirical Bayes approach to variable selection.     

Bayesian Inference in Hidden Markov Random Fields for Binary Data Defined on Large Lattices

Hidden Markov random fields represent a complex hierarchical model, where the hidden latent process is an undirected graphical structure. Performing inference for such models is difficult primarily because the likelihood of the hidden states is often unavailable. The main contribution of this article is to present approximate methods to calculate the likelihood for large lattices based on exact methods for smaller lattices. We introduce approximate likelihood methods by relaxing some of the dependencies in the latent model, and also by extending tractable approximations to the likelihood, the so-called pseudolikelihood approximations, for a large lattice partitioned into smaller sublattices. Results are presented based on simulated data as well as inference for the temporal-spatial structure of the interaction between up- and down-regulated states within the mitochondrial chromosome of the Plasmodium falciparum organism. Supplemental material for this article is available online.     

Adaptive Bayesian Nonstationary Modeling for Large Spatial Datasets Using Covariance Approximations

Gaussian process models have been widely used in spatial statistics but face tremendous modeling and computational challenges for very large nonstationary spatial datasets. To address these challenges, we develop a Bayesian modeling approach using a nonstationary covariance function constructed based on adaptively selected partitions. The partitioned nonstationary class allows one to knit together local covariance parameters into a valid global nonstationary covariance for prediction, where the local covariance parameters are allowed to be estimated within each partition to reduce computational cost. To further facilitate the computations in local covariance estimation and global prediction, we use the full-scale covariance approximation (FSA) approach for the Bayesian inference of our model. One of our contributions is to model the partitions stochastically by embedding a modified treed partitioning process into the hierarchical models that leads to automated partitioning and substantial computational benefits. We illustrate the utility of our method with simulation studies and the global Total Ozone Matrix Spectrometer (TOMS) data. Supplementary materials for this article are available online.     

Noisy Hamiltonian Monte Carlo for Doubly Intractable Distributions

Hamiltonian Monte Carlo (HMC) has been progressively incorporated within the statistician’s toolbox as an alternative sampling method in settings when standard Metropolis–Hastings is inefficient. HMC generates a Markov chain on an augmented state space with transitions based on a deterministic differential flow derived from Hamiltonian mechanics. In practice, the evolution of Hamiltonian systems cannot be solved analytically, requiring numerical integration schemes. Under numerical integration, the resulting approximate solution no longer preserves the measure of the target distribution, therefore an accept–reject step is used to correct the bias. For doubly intractable distributions—such as posterior distributions based on Gibbs random fields—HMC suffers from some computational difficulties: computation of gradients in the differential flow and computation of the accept–reject proposals poses difficulty. In this article, we study the behavior of HMC when these quantities are replaced by Monte Carlo estimates. Supplemental codes for implementing methods used in the article are available online.