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 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.     

Markov basis and Gröbner basis of Segre–Veronese configuration for testing independence in group-wise selections

We consider testing independence in group-wise selections with some restrictions on combinations of choices. We present models for frequency data of selections for which it is easy to perform conditional tests by Markov chain Monte Carlo (MCMC) methods. When the restrictions on the combinations can be described in terms of a Segre–Veronese configuration, an explicit form of a Gröbner basis consisting of binomials of degree two is readily available for performing a Markov chain. We illustrate our setting with the National Center Test for university entrance examinations in Japan. We also apply our method to testing independence hypotheses involving genotypes at more than one locus or haplotypes of alleles on the same chromosome.     

Estimating the number of zero-one multi-way tables via sequential importance sampling

In 2005, Chen et al. introduced a sequential importance sampling (SIS) procedure to analyze zero-one two-way tables with given fixed marginal sums (row and column sums) via the conditional Poisson (CP) distribution. They showed that compared with Monte Carlo Markov chain (MCMC)-based approaches, their importance sampling method is more efficient in terms of running time and also provides an easy and accurate estimate of the total number of contingency tables with fixed marginal sums. In this paper, we extend their result to zero-one multi-way ( $d$ -way, $d \ge 2$ ) contingency tables under the no $d$ -way interaction model, i.e., with fixed $d-1$ marginal sums. Also, we show by simulations that the SIS procedure with CP distribution to estimate the number of zero-one three-way tables under the no three-way interaction model given marginal sums works very well even with some rejections. We also applied our method to Samson’s monks data set.     

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.     

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.     

Markov chain mixing time on cycles

Mixing time quantifies the convergence speed of a Markov chain to the stationary distribution. It is an important quantity related to the performance of MCMC sampling. It is known that the mixing time of a reversible chain can be significantly improved by lifting, resulting in an irreversible chain, while changing the topology of the chain. We supplement this result by showing that if the connectivity graph of a Markov chain is a cycle, then there is an Ω(n²) lower bound for the mixing time. This is the same order of magnitude that is known for reversible chains on the cycle.     

Hybrid Monte Carlo on Hilbert spaces

The Hybrid Monte Carlo (HMC) algorithm provides a framework for sampling from complex, high-dimensional target distributions. In contrast with standard Markov chain Monte Carlo (MCMC) algorithms, it generates nonlocal, nonsymmetric moves in the state space, alleviating random walk type behaviour for the simulated trajectories. However, similarly to algorithms based on random walk or Langevin proposals, the number of steps required to explore the target distribution typically grows with the dimension of the state space. We define a generalized HMC algorithm which overcomes this problem for target measures arising as finite-dimensional approximations of measures π which have density with respect to a Gaussian measure on an infinite-dimensional Hilbert space. The key idea is to construct an MCMC method which is well defined on the Hilbert space itself.We successively address the following issues in the infinite-dimensional setting of a Hilbert space: (i) construction of a probability measure Π in an enlarged phase space having the target π as a marginal, together with a Hamiltonian flow that preserves Π; (ii) development of a suitable geometric numerical integrator for the Hamiltonian flow; and (iii) derivation of an accept/reject rule to ensure preservation of Π when using the above numerical integrator instead of the actual Hamiltonian flow. Experiments are reported that compare the new algorithm with standard HMC and with a version of the Langevin MCMC method defined on a Hilbert space.     

Another look at independence of hitting place and time for the simple random walk

Results of Samuels and Wendel for the simple random walk with drift on the integers which assert independence of interarrival times at the sets {a?r, a+r}, r=1, 2…, k and the arrival position in the set {a?k,a+k}, where a is the starting point, are reobtained by treating the walk as a Markov chain, and considering related chains conditional on absorption at a specified barrier.     

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.     

Markov Bases of Binary Graph Models

This paper is concerned with the topological invariant of a graph given by the maximum degree of a Markov basis element for the corresponding graph model for binary contingency tables. We describe a degree four Markov basis for the model when the underlying graph is a cycle and generalize this result to the complete bipartite graph K_2,n. We also give a combinatorial classification of degree two and three Markov basis moves as well as a Buchberger-free algorithm to compute moves of arbitrary given degree. Finally, we compute the algebraic degree of the model when the underlying graph is a forest.AMS Subject Classification: 05C99, 13P10, 62Q05.     

Markov bases and subbases for bounded contingency tables

In this paper we study the computation of Markov bases for contingency tables whose cell entries have an upper bound. It is known that in this case one has to compute universal Gröbner bases, and this is often infeasible also in small- and medium-sized problems. Here we focus on bounded two-way contingency tables under independence model. We show that when these bounds on cells are positive the set of basic moves of all 2 × 2 minors connects all tables with given margins. We also give some results about bounded incomplete table and we conclude with an open problem on the necessary and sufficient condition on the set of structural zeros so that the set of basic moves of all 2 × 2 minors connects all incomplete contingency tables with given positive margins.     

Sparse Variational Analysis of Linear Mixed Models for Large Data Sets

It is increasingly common to be faced with longitudinal or multi-level data sets that have large numbers of predictors and/or a large sample size. Current methods of fitting and inference for mixed effects models tend to perform poorly in such settings. When there are many variables, it is appealing to allow uncertainty in subset selection and to obtain a sparse characterization of the data. Bayesian methods are available to address these goals using Markov chain Monte Carlo (MCMC), but MCMC is very computationally expensive and can be infeasible in large p and/or large n problems. As a fast approximate Bayes solution, we recommend a novel approximation to the posterior relying on variational methods. Variational methods are used to approximate the posterior of the parameters in a decomposition of the variance components, with priors chosen to obtain a sparse solution that allows selection of random effects. The method is evaluated through a simulation study, and applied to an epidemiological application.     

Multicanonical MCMC for sampling rare events: an illustrative review

Multicanonical MCMC (Multicanonical Markov Chain Monte Carlo; Multicanonical Monte Carlo) is discussed as a method of rare event sampling. Starting from a review of the generic framework of importance sampling, multicanonical MCMC is introduced, followed by applications in random matrices, random graphs, and chaotic dynamical systems. Replica exchange MCMC (also known as parallel tempering or Metropolis-coupled MCMC) is also explained as an alternative to multicanonical MCMC. In the last section, multicanonical MCMC is applied to data surrogation; a successful implementation in surrogating time series is shown. In the appendix, calculation of averages and normalizing constant in an exponential family, phase coexistence, simulated tempering, parallelization, and multivariate extensions are discussed.     

Bayesian Model Choice of Grouped t-Copula

One of the most popular copulas for modeling dependence structures is t-copula. Recently the grouped t-copula was generalized to allow each group to have one member only, so that a priori grouping is not required and the dependence modeling is more flexible. This paper describes a Markov chain Monte Carlo (MCMC) method under the Bayesian inference framework for estimating and choosing t-copula models. Using historical data of foreign exchange (FX) rates as a case study, we found that Bayesian model choice criteria overwhelmingly favor the generalized t-copula. In addition, all the criteria also agree on the second most likely model and these inferences are all consistent with classical likelihood ratio tests. Finally, we demonstrate the impact of model choice on the conditional Value-at-Risk for portfolios of six major FX rates.     

基于MCMC的随机库存系统中提前期需求分布的计算

在需求和提前期均是随机的库存系统中,提前期需求的分布是由提前期分布与需求分布复合而成的,这个复合分布的计算通常是困难的。本文采用基于Gibbs抽样的马尔科夫链蒙特卡洛(MCMC,Markov chain Monte Carlo)方法,抽取条件分布样本作为提前期需求分布的样本,通过样本来计算提前期需求分布密度、服务水平和损失函数。这种方法避免了直接求解复杂积分计算上的困难,也克服了近似分布拟合偏差过大的问题,有效地解决了随机需求与随机提前期的复杂库存系统中提前期需求确定问题。理论与数值分析结果表明:与现有方法相比较,基于MCMC的方法具有计算简便、拟合精度高、通用性好等特点。     

Adaptive Markov Chain Monte Carlo for Bayesian Variable Selection

We describe adaptive Markov chain Monte Carlo (MCMC) methods for sampling posterior distributions arising from Bayesian variable selection problems. Point-mass mixture priors are commonly used in Bayesian variable selection problems in regression. However, for generalized linear and nonlinear models where the conditional densities cannot be obtained directly, the resulting mixture posterior may be difficult to sample using standard MCMC methods due to multimodality. We introduce an adaptive MCMC scheme that automatically tunes the parameters of a family of mixture proposal distributions during simulation. The resulting chain adapts to sample efficiently from multimodal target distributions. For variable selection problems point-mass components are included in the mixture, and the associated weights adapt to approximate marginal posterior variable inclusion probabilities, while the remaining components approximate the posterior over nonzero values. The resulting sampler transitions efficiently between models, performing parameter estimation and variable selection simultaneously. Ergodicity and convergence are guaranteed by limiting the adaptation based on recent theoretical results. The algorithm is demonstrated on a logistic regression model, a sparse kernel regression, and a random field model from statistical biophysics; in each case the adaptive algorithm dramatically outperforms traditional MH algorithms. Supplementary materials for this article are available online.     

Generalized Nonparametric Mixed Effects Models

Generalized linear mixed effects models (GLMM) provide useful tools for correlated and/or over-dispersed non-Gaussian data. This article considers generalized nonparametric mixed effects models (GNMM), which relax the rigid linear assumption on the conditional predictor in a GLMM. We use smoothing splines to model fixed effects. The random effects are general and may also contain stochastic processes corresponding to smoothing splines. We show how to construct smoothing spline ANOVA (SS ANOVA) decompositions for the predictor function. Components in a SS ANOVA decomposition have nice interpretations as main effects and interactions. Experimental design considerations help determine which components are fixed or random. We estimate all parameters and spline functions using stochastic approximation with Markov chain Monte Carlo (MCMC). As iteration increases we increase the MCMC sample size and decrease the step-size of the parameter update. This approach guarantees convergence of the estimates to the expected fixed points. We evaluate our methods through a simulation study.     

A Markov chain Monte Carlo Expectation Maximization Algorithm for Statistical Analysis of DNA Sequence Evolution with Neighbor-Dependent Substitution Rates

The evolution of DNA sequences can be described by discrete state continuous time Markov processes on a phylogenetic tree. We consider neighbor-dependent evolutionary models where the instantaneous rate of substitution at a site depends on the states of the neighboring sites. Neighbor-dependent substitution models are analytically intractable and must be analyzed using either approximate or simulation-based methods. We describe statistical inference of neighbor-dependent models using a Markov chain Monte Carlo expectation maximization (MCMC-EM) algorithm. In the MCMC-EM algorithm, the high-dimensional integrals required in the EM algorithm are estimated using MCMC sampling. The MCMC sampler requires simulation of sample paths from a continuous time Markov process, conditional on the beginning and ending states and the paths of the neighboring sites. An exact path sampling algorithm is developed for this purpose.