Counting Contingency Tables via Multistage Markov Chain Monte Carlo

Abstract

This article describes a multistage Markov chain Monte Carlo (MSMCMC) procedure for estimating the count, , where denotes the set of all a × b contingency tables with specified row and column sums and m total entries. On each stage s = 1, …, r, Hastings–Metropolis (HM) sampling generates states with equilibrium distribution , with inverse-temperature schedule β1 = 0 < β2 < ⋅⋅⋅ < β r < β r + 1 = ∞; nonnegative penalty function H, with global minima H(x) = 0 only for ; and superset , which facilitates sampling by relaxing column sums while maintaining row sums. Two kernels are employed for nominating states. For small β s , one admits moderate-to-large penalty changes. For large β s , the other allows only small changes. Neither kernel admits local minima, thus avoiding an impediment to convergence. Preliminary sampling determines when to switch from one kernel to the other to minimize relative error. Cycling through stages in the order r to 1, rather than 1 to r, speeds convergence for large β s . A comparison of estimates for examples, whose dimensions range from 15 to 24, with exact counts based on Barvinok’s algorithm and estimates based on sequential importance sampling (SIS) favors these alternatives. However, the comparison strongly favors MSMCMC for an example with 64 dimensions. It estimated count with 3.3354 × 10− 3 relative error, whereas the exact counting method was unable to produce a result in more than 182 CPU (computational) days of execution, and SIS would have required at least 42 times as much CPU time to generate an estimate with the same relative error. This latter comparison confirms the known limitations of exact counting methods and of SIS for larger-dimensional problems and suggests that MSMCMC may be a suitable alternative. Proofs not given in the article appear in the Appendix in the online supplemental materials.