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.     

Computing p-values in conditional independence models for a contingency table

We present a Markov chain Monte Carlo (MCMC) method for generating Markov chains using Markov bases for conditional independence models for a four-way contingency table. We then describe a Markov basis characterized by Markov properties associated with a given conditional independence model and show how to use the Markov basis to generate random tables of a Markov chain. The estimates of exact p-values can be obtained from random tables generated by the MCMC method. Numerical experiments examine the performance of the proposed MCMC method in comparison with the χ ² approximation using large sparse contingency tables.     

Implementing pure adaptive search for global optimization using Markov chain sampling

The Pure Adaptive Search (PAS) algorithm for global optimization yields a sequence of points, each of which is uniformly distributed in the level set corresponding to its predecessor. This algorithm has the highly desirable property of solving a large class of global optimization problems using a number of iterations that increases at most linearly in the dimension of the problem. Unfortunately, PAS has remained of mostly theoretical interest due to the difficulty of generating, in each iteration, a point uniformly distributed in the improving feasible region. In this article, we derive a coupling equivalence between generating an approximately uniformly distributed point using Markov chain sampling, and generating an exactly uniformly distributed point with a certain probability. This result is used to characterize the complexity of a PAS-implementation as a function of (a) the number of iterations required by PAS to achieve a certain solution quality guarantee, and (b) the complexity of the sampling algorithm used. As an application, we use this equivalence to show that PAS, using the so-called Random ball walk Markov chain sampling method for generating nearly uniform points in a convex region, can be used to solve most convex programming problems in polynomial time.     

Exact sampling with highly uniform point sets

In 1996, Propp and Wilson came up with a remarkably clever method for generating exact samples from the stationary distribution of a Markov chain [J.G. Propp, D.B. Wilson, Exact sampling with coupled Markov chains and applications to statistical mechanics, Random Structures and Algorithms 9 (1–2) (1996) 223–252]. Their method, called “perfect sampling” or “exact sampling” avoids the inherent bias of samples that are generated by running the chain for a large but fixed number of steps. It does so by using a strategy called “coupling from the past”. Although the sampling mechanism used in their method is typically driven by independent random points, more structured sampling can also be used. Recently, Craiu and Meng [R.V. Craiu, X.-L. Meng, Antithetic coupling for perfect sampling, in: E.I. George (Ed.), Bayesian Methods with Applications to Science, Policy, and Official Statistics (Selected Papers from ISBA 2000), 2000, pp. 99–108; R.V. Craiu, X.-L. Meng, Multi-process parallel antithetic coupling for forward and backward Markov Chain Monte Carlo, Annals of Statistics 33 (2005) 661–697] suggested using different forms of antithetic coupling for that purpose. In this paper, we consider the use of highly uniform point sets to drive the exact sampling in Propp and Wilson’s method, and illustrate the effectiveness of the proposed method with a few numerical examples.     

On Markov chain Monte Carlo Algorithms for Computing Conditional Expectations Based on Sufficient Statistics

Much work has focused on developing exact tests for the analysis of discrete data using log linear or logistic regression models. A parametric model is tested for a dataset by conditioning on the value of a sufficient statistic and determining the probability of obtaining another dataset as extreme or more extreme relative to the general model, where extremeness is determined by the value of a test statistic such as the chi-square or the log-likelihood ratio. Exact determination of these probabilities can be infeasible for high dimensional problems, and asymptotic approximations to them are often inaccurate when there are small data entries and/or there are many nuisance parameters. In these cases Monte Carlo methods can be used to estimate exact probabilities by randomly generating datasets (tables) that match the sufficient statistic of the original table. However, naive Monte Carlo methods produce tables that are usually far from matching the sufficient statistic. The Markov chain Monte Carlo method used in this work (the regression/attraction approach) uses attraction to concentrate the distribution around the set of tables that match the sufficient statistic, and uses regression to take advantage of information in tables that “almost” match. It is also more general than others in that it does not require the sufficient statistic to be linear, and it can be adapted to problems involving continuous variables. The method is applied to several high dimensional settings including four-way tables with a model of no four-way interaction, and a table of continuous data based on beta distributions. It is powerful enough to deal with the difficult problem of four-way tables and flexible enough to handle continuous data with a nonlinear sufficient statistic.     

Markov chain analysis of genetic algorithms applied to?fitness functions perturbed concurrently by additive and multiplicative noise

We analyze the transition and convergence properties of genetic algorithms (GAs) applied to fitness functions perturbed concurrently by additive and multiplicative noise. Both additive noise and multiplicative noise are assumed to take on finitely many values. We explicitly construct a Markov chain that models the evolution of GAs in this noisy environment and analyze it to investigate the algorithms. Our analysis shows that this Markov chain is indecomposable; it has only one positive recurrent communication class. Using this property, we establish a condition that is both necessary and sufficient for GAs to eventually (i.e., as the number of iterations goes to infinity) find a globally optimal solution with probability 1. Similarly, we identify a condition that is both necessary and sufficient for the algorithms to eventually with probability 1 fail to find any globally optimal solution. Our analysis also shows that the chain has a stationary distribution that is also its steady-state distribution. Based on this property and the transition probabilities of the chain, we compute the exact probability that a GA is guaranteed to select a globally optimal solution upon completion of each iteration.     

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.     

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.     

A Markov chain Monte Carlo Algorithm for Approximating Exact Conditional Probabilities

Conditional inference eliminates nuisance parameters by conditioning on their sufficient statistics. For contingency tables conditional inference entails enumerating all tables with the same sufficient statistics as the observed data. For moderately sized tables and/or complex models, the computing time to enumerate these tables is often prohibitive. Monte Carlo approximations offer a viable alternative provided it is possible to obtain samples from the correct conditional distribution. This article presents an MCMC extension of the importance sampling algorithm, using a rounded normal candidate to update randomly chosen cells while leaving the remainder of the table fixed. This local approximation can greatly increase the efficiency of the rounded normal candidate. By choosing the number of cells to be updated at random, a balance is struck between dependency in the Markov chain and accuracy of the candidate.     

Time-inhomogeneous discrete stochastic search methods for optimal Bernoulli parameters

We present two time-inhomogeneous search processes for finding the optimal Bernoulli parameters, where the performance measure cannot be evaluated exactly but must be estimated through Monte Carlo simulation. At each iteration, two neighbouring alternatives are compared and the one that appears to be better is passed on to the next iteration. The first search process uses an increasing sample size of each configuration at each iteration. The second search process uses a sequential sampling procedure with increasing boundaries as the number of iterations increases. At each iteration the acceptance of a new configuration depends on the iterate number, therefore, the search process turns out to be inhomogeneous Markov chain. We show that if the increase occurs slower than a certain rate, these search processes will converge to the optimal set with probability one. © 1998 John Wiley & Sons, Ltd.     

Hybrid schemes for exact conditional inference in discrete exponential families

Exact conditional goodness-of-fit tests for discrete exponential family models can be conducted via Monte Carlo estimation of p values by sampling from the conditional distribution of multiway contingency tables. The two most popular methods for such sampling are Markov chain Monte Carlo (MCMC) and sequential importance sampling (SIS). In this work we consider various ways to hybridize the two schemes and propose one standout strategy as a good general purpose method for conducting inference. The proposed method runs many parallel chains initialized at SIS samples across the fiber. When a Markov basis is unavailable, the proposed scheme uses a lattice basis with intermittent SIS proposals to guarantee irreducibility and asymptotic unbiasedness. The scheme alleviates many of the challenges faced by the MCMC and SIS schemes individually while largely retaining their strengths. It also provides diagnostics that guide and lend credibility to the procedure. Simulations demonstrate the viability of the approach.     

A Marginal Sampler for σ-Stable Poisson–Kingman Mixture Models

We investigate the class of σ-stable Poisson–Kingman random probability measures (RPMs) in the context of Bayesian nonparametric mixture modeling. This is a large class of discrete RPMs, which encompasses most of the popular discrete RPMs used in Bayesian nonparametrics, such as the Dirichlet process, Pitman–Yor process, the normalized inverse Gaussian process, and the normalized generalized Gamma process. We show how certain sampling properties and marginal characterizations of σ-stable Poisson–Kingman RPMs can be usefully exploited for devising a Markov chain Monte Carlo (MCMC) algorithm for performing posterior inference with a Bayesian nonparametric mixture model. Specifically, we introduce a novel and efficient MCMC sampling scheme in an augmented space that has a small number of auxiliary variables per iteration. We apply our sampling scheme to a density estimation and clustering tasks with unidimensional and multidimensional datasets, and compare it against competing MCMC sampling schemes. Supplementary materials for this article are available online.     

Parameter estimation for partially observable systems subject to random failure

In this paper, we present a parameter estimation procedure for a condition‐based maintenance model under partial observations. Systems can be in a healthy or unhealthy operational state, or in a failure state. System deterioration is driven by a continuous time homogeneous Markov chain and the system state is unobservable, except the failure state. Vector information that is stochastically related to the system state is obtained through condition monitoring at equidistant sampling times. Two types of data histories are available — data histories that end with observable failure, and censored data histories that end when the system has been suspended from operation but has not failed. The state and observation processes are modeled in the hidden Markov framework and the model parameters are estimated using the expectation–maximization algorithm. We show that both the pseudolikelihood function and the parameter updates in each iteration of the expectation–maximization algorithm have explicit formulas. A numerical example is developed using real multivariate spectrometric oil data coming from the failing transmission units of 240‐ton heavy hauler trucks used in the Athabasca oil sands of Alberta, Canada. Copyright © 2012 John Wiley & Sons, Ltd.     

Scaled Optimal Path Trust-Region Algorithm

Trust-region algorithms solve a trust-region subproblem at each iteration. Among the methods solving the subproblem, the optimal path algorithm obtains the solution to the subproblem in full-dimensional space by using the eigenvalues and eigenvectors of the system. Although the idea is attractive, the existing optimal path method seems impractical because, in addition to factorization, it requires either the calculation of the full eigensystem of a matrix or repeated factorizations of matrices at each iteration. In this paper, we propose a scaled optimal path trust-region algorithm. The algorithm finds a solution of the subproblem in full-dimensional space by just one Bunch–Parlett factorization for symmetric matrices at each iteration and by using the resulting unit lower triangular factor to scale the variables in the problem. A scaled optimal path can then be formed easily. The algorithm has good convergence properties under commonly used conditions. Computational results for small-scale and large-scale optimization problems are presented which show that the algorithm is robust and effective.     

Interruptible Exact Sampling in the Passive Case

We establish, for various scenarios, whether or not interruptible exact stationary sampling is possible when a finite-state Markov chain can only be viewed passively. In particular, we prove that such sampling is not possible using a single copy of the chain. Such sampling is possible when enough copies of the chain are available, and we provide an algorithm that terminates with probability one.     

Iterative solution of inconsistent systems of linear inequalities

In 1980 S.P. Han proposed a finitely terminating (in exact arithmetic) algorithm for solving an inconsistent system of linear inequalities in a least squares sense. This algorithm uses a singular value decomposition of a submatrix of the problem matrix on each iteration, making it impractical for all but only smaller problems. In this paper we show that a modification of Han's algorithm allows us to introduce an iterative approximation to the singular value decomposition solution. (© 2013 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Optimizing random scan Gibbs samplers

The Gibbs sampler is a popular Markov chain Monte Carlo routine for generating random variates from distributions otherwise difficult to sample. A number of implementations are available for running a Gibbs sampler varying in the order through which the full conditional distributions used by the Gibbs sampler are cycled or visited. A common, and in fact the original, implementation is the random scan strategy, whereby the full conditional distributions are updated in a randomly selected order each iteration. In this paper, we introduce a random scan Gibbs sampler which adaptively updates the selection probabilities or “learns” from all previous random variates generated during the Gibbs sampling. In the process, we outline a number of variations on the random scan Gibbs sampler which allows the practitioner many choices for setting the selection probabilities and prove convergence of the induced (Markov) chain to the stationary distribution of interest. Though we emphasize flexibility in user choice and specification of these random scan algorithms, we present a minimax random scan which determines the selection probabilities through decision theoretic considerations on the precision of estimators of interest. We illustrate and apply the results presented by using the adaptive random scan Gibbs sampler developed to sample from multivariate Gaussian target distributions, to automate samplers for posterior simulation under Dirichlet process mixture models, and to fit mixtures of distributions.     

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.     

Inexact non-interior continuation method for solving large-scale monotone SDCP

For exact Newton method for solving monotone semidefinite complementarity problems (SDCP), one needs to exactly solve a linear system of equations at each iteration. For problems of large size, solving the linear system of equations exactly can be very expensive. In this paper, we propose a new inexact smoothing/continuation algorithm for solution of large-scale monotone SDCP. At each iteration the corresponding linear system of equations is solved only approximately. Under mild assumptions, the algorithm is shown to be both globally and superlinearly convergent.