\section{Introduction}

Conducting the goodness of fit test for discrete exponential family models with sparse contingency tables is a challenging task since we cannot use asymptotic methods, especially when we have large tables.  To do so, one applies sampling procedures, such as Markov Chain Monte Carlo (MCMC) procedures, from the state space \cite{Besag:1989uq,chen2005,Diaconis:1985kx,Guo:1992fk}. 

Even though Besag and Clifford argued that irreducibility of a Markov chain is not essential in \cite{Besag:1989uq}, that view is not conventional since it is very difficult to understand a chain over sparse tables in high dimensions.
To understand the irreducibility of a Markov chain over the set of all possible contingency tables with given constraints, Diaconis and Sturmfels introduced a notion of a {\em Markov basis} \cite{DS}.  Diaconis and Sturmfels showed that a Markov basis, a generating set of a toric ideal associate with the design matrix for a discrete exponential family model, is the set of moves which guarantee to connect all tables for the given model \cite{DS}.  The most significant benefit in computing a Markov basis for a discrete exponential family model is that if we have a Markov basis for the model, we do not have to re-compute it since moves in a Markov basis connect all tables for any observed tables under the model.  However, it is computationally expensive to compute a Markov basis in terms of the dimension of tables.  In fact, \cite{deloera05} showed that the number of elements in a Markov basis for the no-three-way interaction model can be arbitrarily large.   This means that it is not feasible to compute a Markov basis for large tables under the no-three-way interaction model. Therefore, in this research, we focus on the connectivity of a Markov chain for $I \times J \times K$ contingency tables under the no-three-way interaction model. Hara et al.~conducted study of MCMC with a {\em lattice basis} which is much easier to compute than a Markov basis \cite{Hara2011RunningMC}.  However, with a lattice basis, we cannot guarantee that a chain is connected so that estimated distribution of statistics might be biased.  

To try to circumvent the difficulty of computing a
Markov basis which may be large,
\cite{chen-dinwoodie-dobra-huber2005,Bunea:2000fk} 
studied computing a smaller set of moves, such as a lattice basis, by allowing cell counts of the
contingency table to be negative one.  Motivated by their work, in this paper, we focus on a Markov chain using a set of all {\em basic moves} by allowing cell counts to be $-1$.  A basic move of all $2 \times 2\times 2$ minors is a $I \times J \times K$ contingency table with $\{0, \pm 1\}$ entries such that for distinct $i, i' \in \{1, \ldots , I\}$, distinct $j, j' \in \{1, \ldots , J\}$, and distinct $k, k' \in \{1, \ldots , K\}$,  
\begin{eqnarray*}
i & \begin{array}{|c||c|c|}\hline
 & j & j'\\\hline \hline
k & 1 & -1\\\hline
k' & -1 & 1\\\hline
\end{array}\\
i' & \begin{array}{|c||c|c|}\hline
 & j & j'\\\hline \hline
k & -1 & 1\\\hline
k' & 1 & -1\\\hline
\end{array}\\
\end{eqnarray*}
and all zeros otherwise.  One can show that if $I = 2$, $J \geq 2$, and $K \geq 2$, then basic moves of all $2 \times 2 \times 2$ minors connect all tables if we allow cell counts to be negative one. In 2018, Lee showed that if $I = 3$, $J = 3$, and $K \geq 2$, then 
basic moves connect all tables if we allow cell counts $X_{ijk}$ to be $-1$ when we run a Markov chain. For $I > 3$, $J > 3$, and $K \geq 2$, it is still an open problem that basic moves connect all tables if we allow cell counts $X_{ijk}$ to be $-1$ when we run a Markov chain. 

The trade off to this approach is longer
running time of the Markov chains.  
Even using a standard MCMC approach, to sample a table independently from the distribution,
Markov chains can take a long time to 
converge to a stationary distribution in order to satisfy the
independent assumption. 
Since we allow cell counts to drop to $-1$ when we run a Markov chain, it is possible to slow down a mixing time of a Markov chain.  
In order to speed up mixing time we propose to combine a {\em simulated annealing (SA)} approach \cite{doi:10.1126/science.220.4598.671} with a MCMC with basic moves allowing cell counts to be negative one. With SA, we can push a chain more inside of the state space, called a {\em fiber}, of contingency tables with fixed margins so that we can avoid having a fat tails of a distribution of $\chi^2$ statistics.  

In 2018, Kahle, Yoshida, and Garcia-Puente introduced a hybrid scheme of Sequential Importance Sampling (SIS) and MCMC to speed up mixing time and to have a connectivity of a chain \cite{KYG}.  However, from our initial experiments discussed in Section \ref{sec:computation}, SIS keeps sending a Markov chain to near boundary of the state space from another location near the boundary.  Thus, in this paper, we investigate how the hybrid scheme of SIS and MCMC with the set of basic moves of all $2\times 2\times 2$ minors works for $I \times J \times K$ tables under the no-three-way interaction model.
In addition, in order to speed up a mixing time under the no-three-way interaction model, we propose a hybrid procedure of MCMC and Hit and Run (HAR) which Andersen and Diaconis proposed in \cite{AD}. The HAR procedure is a MCMC sampler which sequentially samples points from the state space by drawing a point randomly from a ``line'' in the space.  Andersen and Diaconis suggested in \cite{AD} that we randomly choose a $2 \times 2 \times 2$ sub-table from the observed table and we randomly choose a table from a ``line'' defined by all margins computed from the $2 \times 2 \times 2$ sub-table. 
Computational experiments with simulated and empirical data sets, this hybrid procedure of MCMC and HAR works very well (Section \ref{sec:computation}).  

This paper is organized as follows:  In Section \ref{sec:basic}, we set up basics for our results in this paper.  
Then, in Section \ref{sec:hybrid}, we propose a hybrid scheme between MCMC with basic moves allowing cell counts to be negative one with SA and HAR sampling. Then we apply our hybrid scheme to simulated data and empirical data on U.S. Navy (USN) personnel in Section \ref{sec:computation}.  We end with discussion and open problems.  

\section{Preliminaries}\label{sec:basic}
In this section we set up some basic notation and definitions for our results in this paper.  

Let $A \in \mathbb Z^{m\times d}$ and $b \in \mathbb Z^m$.  We define $b-$fiber of
$A$ is
\[
\mathcal{F}_{A, b}:= \{ u\in \mathbb Z_{\geq 0}^d: A u = b\}.  
\]
For a finite
$\mathcal{M} \subset \ker A \cap \mathbb Z^n$, we define the fiber
graph $\mathcal{G}(\mathcal{F}_{A,b}, \mathcal{M})$ of $A$ as graph having the
$b$-fiber as its vertex set and the set 
\[
\{(u, v): u, v \in \mathcal{F}_{A, b}, \, (u - v) \in \pm \mathcal{M}\} 
\]
as its edge set.

\begin{definition}
We fix a matrix $A \in \mathbb Z^{m\times d}$.  We call $\mathcal{M} \subset \ker A \cap \mathbb Z^n$ as a {\em Markov basis} for a matrix $A$ if the fiber
graph $\mathcal{G}(\mathcal{F}_{A,b}, \mathcal{M})$ of $A$ is connected for any $b \in \mathbb Z^m$.
\end{definition}

In this paper, we focus on assessing interaction between three discrete random variables
$X, \, Y, \, Z$ with finite levels, where $X \in \{1, \ldots , I\}$,
$Y \in \{1, \ldots , J\}$, and $Z \in \{1, \ldots , K\}$ for $I, \,
J, \, K \in \mathbb{N}$.  In order to conduct statistical analysis on such
interactions,  a commonly used model is called a {\em log linear
  model}.  A log linear model is a generalized linear model and it
forms a class of discrete exponential family.  

\begin{definition}
Suppose $A \in \mathbb Z^{m \times d}$ be a non-negative integral matrix such
that 
\[
\sum_{i = 1}^m a_{i1} = \ldots = \sum_{i = 1}^m a_{id}
\]
where $a_{ij}$ is the $(i, j)$th element of the matrix $A$.  
If its likelihood function $L(\theta| \bf{n})$ for an observation $\bf{n}$ can be written as
\[
L(\theta| \bf{n}) = c \cdot \theta^{A\bf{n}}
\]
where $c$ is a normalized constant, then the model is called a {\em
  log linear model} and we call such matrix $A$ as the {\em design matrix} for a model.
\end{definition}

Applying the notation of the log-linear model notation to the
three-way contingency table, we have the parametalization $n_{ijk}
\sim Pois(\mu_{ijk})$, the Poisson distribution with the parameter $\mu_{ijk}$, such that:
\[
\log(\mu_{ijk})=\lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{XZ}+\lambda_{jk}^{YZ}+\lambda_{ijk}^{XYZ}.
\]
If we assume there is no interaction between three categorical
variables, namely, 
\[
\lambda_{ijk}^{XYZ} = 0,
\]
for all $i = 1, \ldots , I$, $j = 1, \ldots , J$, and $k = 1, \ldots ,
K$, then we call this model as 
the {\em no-three-way interaction model}.

It is easy to show that under the no-three-way interaction model, we have 
the $b-$fiber of $A$ with a given observed table ${\bf n}$ is 
\[
\mathcal{F}_{A, b}:= \left\{{\bf u} \Big| \sum_{i = 1}^m u_{ijk} = \sum_{i = 1}^m n_{ijk}, \, \sum_{j = 1}^n
u_{ijk} = \sum_{j = 1}^n n_{ijk}, \, \sum_{k = 1}^l u_{ijk} = \sum_{k
  = 1}^l n_{ijk} \right\} =: \{u \in \mathbb Z^{IJK}_+| Au = b\},
\]
where $n_{ijk}$ is the cell count for the $(i, j, k)$th cell in ${\bf n}$.
\begin{example}
Suppose we have a $2 \times 2 \times 2$ table such that
\[
\begin{array}{|c|c|}\hline
     X_{111}&X_{112}  \\\hline
     X_{121}&X_{122}  \\\hline
\end{array}
\ \
\begin{array}{|c|c|}\hline
     X_{211}&X_{212}  \\\hline
     X_{221}&X_{222}  \\\hline
\end{array} .
\]
The the design matrix under the no-three-way interaction model is 
\[
A = \left(
\begin{array}{cccccccc}
     1 & 1 & 0& 0 &  0 & 0 & 0& 0  \\
     0 & 0 & 1 & 1 & 0& 0 &  0 & 0 \\
     0 & 0 & 0& 0 &  1 & 1 & 0& 0  \\
     0 & 0 & 0& 0 &  0 & 0 & 1 & 1 \\
     1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\
     0 & 1 & 0 & 1 & 0 & 0 & 0 & 0\\
     0 & 0 & 0 & 0 & 1 & 0 & 1 & 0\\
     0 & 0 & 0 & 0 & 0 & 1 & 0 & 1\\
     1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\
     0 & 1 & 0 & 0 & 0& 1 &  0 & 0\\
     0 & 0 & 1 & 0 & 0 & 0 & 1 & 0\\
     0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 \\
\end{array}
\right) .
\]
Then, suppose we have a $2 \times 2 \times 2$ table such that
\[
\begin{array}{|c|c|}\hline
     1&1 \\\hline
     1&1  \\\hline
\end{array}
\ \
\begin{array}{|c|c|}\hline
     1&1 \\\hline
     1&1 \\\hline
\end{array} .
\]
Then the vector of margins is
\[
b = (2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)^T ,
\]
where $x^T$ is the transpose of a vector $x \in \mathbb{R}^m$.
\end{example}


Now we define a ``basic move''.
\begin{definition}
Let $m$ be a $I \times J \times K$ table such that  
\begin{eqnarray*}
i & \begin{array}{|c||c|c|}\hline
 & j & j'\\\hline \hline
k & 1 & -1\\\hline
k' & -1 & 1\\\hline
\end{array}\\
i' & \begin{array}{|c||c|c|}\hline
 & j & j'\\\hline \hline
k & -1 & 1\\\hline
k' & 1 & -1\\\hline
\end{array}\\
\end{eqnarray*}
where $1 \leq i, i' \leq I, \, 1 \leq j, j' \leq
J, \, 1 \leq k, k' \leq K$, $i \not = i'$, $j \not = j'$, and $k \not
= k'$, and other cells are all zero.  We call $m$ a {\em basic move}. 
We denote this table $m$ as $(i, i'; j, j'; k, k')$.
\end{definition}

\begin{example}
Suppose we have a $3 \times 3\times 3$ table such that
\[
\begin{array}{|c|c|c|}\hline
     3t&0 & 0 \\\hline
     0&3t & 0  \\\hline
     0 & 0 & 3t\\\hline
\end{array}
\ \
\begin{array}{|c|c|c|}\hline
     0&3t & 0 \\\hline
     0&0 & 3t  \\\hline
     3t & 0 & 0\\\hline
\end{array} 
\ \
\begin{array}{|c|c|c|}\hline
     0&0 & 3t \\\hline
     3t&0 & 0  \\\hline
     0 & 3t & 0\\\hline
\end{array} ,
\]
where $t = 1, 2, \ldots$.  Then note that a set of basic moves does not connect this table to any other tables with the same margins.  However, it is easy to show that if we allow cell counts to be $-1$, then we can connect this table to other tables in the fiber $\mathcal{F}_{A, b}$.
\end{example}

\begin{remark}
By Lemma 12 in \cite{Sturmfels1995GrobnerBA}, for some positive integer $t$, if we allow $X_{ijk} \geq -t$, then all moves in $M$ connect all
tables in $\mathcal{F}_{A, b}$ for a three way $I\times J\times K$ contingency table
$X_{ijk}$ under the no-three way interaction model.
\end{remark}

\begin{remark}
Suppose ${\bf b}={\bf b}_+ - {\bf b}_-$ is in a Markov basis for a matrix $A$. If we can show that a set of basic moves connect from ${\bf b}_+$ and ${\bf b}_-$ for any move ${\bf b}={\bf b}_+ - {\bf b}_-$ in a Markov basis by allowing cell counts to be negative one, then basic moves connect all tables in $\mathcal{F}_{A, b}$ for any $b \in \mathbb{Z}^m_{\geq 0}$ by allowing cell counts to be $-1$.
\end{remark}

\begin{example}
First we consider a $3 \times 4 \times 6$ table and we consider an
``indispensable move'', which is a necessary move in all Markov bases for a fixed matrix $A$, such that
\begin{eqnarray*}
& \begin{array}{|c|c|c|c|c|c|}\hline
+1 & -1 & 0 & 0 & 0 & 0\\\hline 
0 & +1 & -1 & 0 & 0 & 0\\\hline
0 & 0 & +1 & 0 & 0 & -1\\\hline
-1 & 0 & 0 & 0 & 0 & +1\\\hline
\end{array} \, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & +1 & 0 & -1 & 0 & 0\\\hline 
0 & -1 & 0 & 0 & +1 & 0\\\hline
0 & 0 & 0 & +1 & 0 & -1\\\hline
0 & 0 & 0 & 0 & -1 & +1\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
-1 & 0 & 0 & +1 & 0 & 0\\\hline 
0 & 0 & +1 & 0 & -1 & 0\\\hline
0 & 0 & -1 & -1 & 0 & +2\\\hline
+1 & 0 & 0 & 0 & +1 & -2\\\hline
\end{array}\\
= & 
\begin{array}{|c|c|c|c|c|c|}\hline
1 & 0 & 0 & 0 & 0 & 0\\\hline 
0 & 1 & 0 & 0 & 0 & 0\\\hline
0 & 0 & 1 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 1\\\hline
\end{array} \, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 1 & 0 & 0 & 0 & 0\\\hline 
0 & 0 & 0 & 0 & 1 & 0\\\hline
0 & 0 & 0 & 1 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 1\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 0 & 0 & 1 & 0 & 0\\\hline 
0 & 0 & 1 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 2\\\hline
1 & 0 & 0 & 0 & 1 & 0\\\hline
\end{array}\\
- & 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 1 & 0 & 0 & 0 & 0\\\hline 
0 & 0 & 1 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 1\\\hline
1 & 0 & 0 & 0 & 0 & 0\\\hline
\end{array} \, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 0 & 0 & 1 & 0 & 0\\\hline 
0 & 1 & 0 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 1\\\hline
0 & 0 & 0 & 0 & 1 & 0\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
1 & 0 & 0 & 0 & 0 & 0\\\hline 
0 & 0 & 0 & 0 & 1 & 0\\\hline
0 & 0 & 1 & 1 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 2\\\hline
\end{array}\\
=: & {\bf b^1_{+}} - {\bf b^1_{-}}.
\end{eqnarray*}

If we allow $X_{ijk} \geq -1$ for $1 \leq i \leq I, 1 \leq j \leq J, 1\leq
k \leq K$, then we can show that 
\[
 {\bf b^1_{-}} =  {\bf b^1_{+}} - (2,3;3, 4; 5,6) -(1,3;3,4;1,6) +
 (2,3;2,3;3,5) - (2,3;1,3;1,4) - (1,2;2,3;2,3) - (1,2;1,3;1,2).
\]

In addition we consider another indispensable move such that
\begin{eqnarray*}
& \begin{array}{|c|c|c|c|c|c|}\hline
+1 & -1 & 0 & 0 & 0 & 0\\\hline 
0 & +1 & -1 & 0 & 0 & 0\\\hline
0 & 0 & +1 & -1 & 0 & 0\\\hline
-1 & 0 & 0 & +1 & 0 & 0\\\hline
\end{array} \, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
-1 & 0 & 0 & 0 & 0 & +1\\\hline 
0 & 0 & +1 & 0 & -1 & 0\\\hline
0 & 0 & -1 & 0 & 0 & +1\\\hline
+1 & 0 & 0 & 0 & +1 & -2\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & +1 & 0 & 0 & 0 & -1\\\hline 
0 & -1 & 0 & 0 & +1 & 0\\\hline
0 & 0 & 0 & +1 & 0 & -1\\\hline
0 & 0 & 0 & -1 & -1 & +2\\\hline
\end{array}\\
= & 
\begin{array}{|c|c|c|c|c|c|}\hline
1 & 0 & 0 & 0 & 0 & 0\\\hline 
0 & 1 & 0 & 0 & 0 & 0\\\hline
0 & 0 & 1 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 1 & 0 & 0\\\hline
\end{array} \, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 0 & 0 & 0 & 0 & 1\\\hline 
0 & 0 & 1 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 1\\\hline
1 & 0 & 0 & 0 & 1 & 0\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 1 & 0 & 0 & 0 & 0\\\hline 
0 & 0 & 0 & 0 & 1 & 0\\\hline
0 & 0 & 0 & 1 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 2\\\hline
\end{array}\\
- & 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 1 & 0 & 0 & 0 & 0\\\hline 
0 & 0 & 1 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 1 & 0 & 0\\\hline
1 & 0 & 0 & 0 & 0 & 0\\\hline
\end{array} \, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
1 & 0 & 0 & 0 & 0 & 0\\\hline 
0 & 0 & 0 & 0 & 1 & 0\\\hline
0 & 0 & 1 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 2\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 0 & 0 & 0 & 0 & 1\\\hline 
0 & 1 & 0 & 0 & 0 & 0\\\hline
0 & 0 & 0 & 0 & 0 & 1\\\hline
0 & 0 & 0 & 1 & 1 & 0\\\hline
\end{array}\\
=: & {\bf b^2_{+}} - {\bf b^2_{-}}.
\end{eqnarray*}

If we allow $X_{ijk} \geq -1$ for $1 \leq i \leq I, 1 \leq j \leq J, 1\leq
k \leq K$, then we can show that 
\[
 {\bf b^2_{-}} =  {\bf b^2_{+}} + (2,3;2, 4; 5,6) +(2,3;3,4;4,6) +
 (2,3;1,2;2,6) - (1,2;3,4;1,4) - (1,2;1,3;1,3) + (1,2;1,2;2,3).
\]

\end{example}

\begin{remark}
One notices that by ways we use basic moves to connect from ${\bf
  b_+^1}$ to ${\bf b_-^1}$ and ${\bf b_+^2}$ to ${\bf b_-^2}$ we
cannot connect by allowing just one cell to be $-1$.  We have to allow at
least three cells to be $-1$.  So it is interesting to find out
the upper bounds $K$ such that all moves in $M$ connect all
tables in $F$ if  we allow at most $K$ cells $X_{ijk} \geq -1$.
\end{remark}

\begin{remark}
Dinwoodie et al.~discussed connecting tables via allowing cell counts to be negative one and they discussed a condition for a set of moves so that it connects tables in the fiber by allowing cell counts to be negative one (Proposition 0.2.1 \cite{Ian2008}).  However, for the no-three-way interaction model, the {\em ideal} generated by binomials with exponents for the positive part and negative part of all possible basic moves for $3 \times 3 \times 4$ is not radical.  Therefore, we cannot use Proposition 0.2.1 for the connectivity of tables with allowing cell counts to be $-1$. However, using a known Markov basis for $3 \times 3 \times K$, for $K \geq 2$, by Aoki and Takemura in \cite{aoki-takemura-2003anz}, we can show that basic moves of all $2 \times 2 \times 2$ minors can connect tables for $3 \times 3 \times 4$ case.
\end{remark}

In general we do not know that a set of basic moves of $2 \times 2 \times 2$ minors connect all tables in $\mathcal{F}_{A, b}$ under the no-three-way interaction model for any $b \in \mathbb{R}^{m}_{\geq 0}$ if we allow cell counts to be $-1$.  


\section{Hybrid Scheme of MCMC with SA and HAR}\label{sec:hybrid}

It is well-known that MCMC with a Markov basis for a discrete exponential family model does not mix rapidly \cite{tobias}.
Thus, in this section, we propose hybrid schemes of MCMC with other sampling methodologies including SA, SIS and HAR.  

Recall that if we do not allow cell counts to be $-1$, then a set of basic moves cannot connect all tables in $\mathcal{F}_{A, b}$.  Thurs, throughout this paper, we assume that we allow cell counts to be $-1$ to make sure a chain is connected. 

The Hit-and Run (HAR) algorithm is a Markov Chain Monte Carlo method to sample a point from a closed bounded set from $\mathbb{R}^d$ \cite{Zabinsky2013}.  The HAR algorithm to sample a point from a set $S \subset \mathbb{R}^d$ is shown in Algorithm \ref{alg:HAR}.


\begin{algorithm}[H]
\centering
\fbox{\parbox{0.91\linewidth}{
\scalebox{.82}{
\begin{minipage}{.9\linewidth}
\KwData{Initial point $x_0$ $\in$ $S$ and maximum iteration $I\ge1$
}
\KwResult{A random point $x$ $\in$ $S$.}

 $k \gets 0$;

 \For{$k = 0$ \KwTo $I-1$}{
 Generate random direction $D^k$ uniformly distributed over surface of a unit hypersphere centered around $x_k$\
 
 Generate a random point $x_{k+1}$ from a line $L_k:={y \in S : y = x_k + \lambda, \lambda \in D^k}$}


 Return $x:=x_I$\;
\end{minipage}%
    }
    }}
    \caption{Hit-and-Run algorithm \cite{Zabinsky2013}}\label{alg:HAR}
\end{algorithm}


Andersen and Diaconis defined the HAR algorithm for sampling three way contingency tables (Section 2.3 in \cite{AD}) under the no-three-way interaction model from its fiber. Since a classical HAR method samples points from a line over $\mathbb{R}^d$ and since the HAR algorithm defined by Andersen and Diaconis samples contingency tables from the fiber which is a lattice point inside of a polytope, we distinguish between the two by naming the HAR algorithm defined in \cite{AD} {\em discrete Hit and Run} (DHR).   Algorithm \ref{alg:DHR4} describes the outline of their DHR algorithm for three way contingency tables.


\newgeometry{top=20mm}
\begin{algorithm}[H]
\centering
\fbox{\parbox{0.9\linewidth}{
\scalebox{.82}{
\begin{minipage}{.9\linewidth}
\KwData{An $I \times J \times K$ contingency table, ${\bf n_0}$ where $n_{ijk}$ is the cell count for the (i,j,k)th cell
}
\KwResult{Sample $\bf z$ from the fiber $\mathcal{F}_{A, b}$, where $b$ is the vector of margins and $A$ is the design matrix for ${\bf n}_0$ without Metropolis weighting.}
\begin{enumerate}
   

 \item Randomly choose $(i,i') \subset I$, $(j,j') \subset J$, and $(k,k') \subset K$;

 \item Let
        \[
        S((i, i'), (j, j'), (k, k')) := \left\{y \in \mathbb Z^8\Big|
        \begin{array}{ccc}
        y_{ijk}+ y_{ijk'}&=& n_{ijk}+n_{ijk'}\\
        y_{ij'k}+ y_{ij'k'}&=& n_{ij'k}+n_{ij'k'}\\
        y_{ijk}+ y_{ij'k}&=& n_{ijk}+n_{ij'k}\\
        y_{ijk'}+ y_{ij'k'}&=& n_{ijk'}+ n_{ij'k'}\\
        y_{i'jk}+ y_{i'jk'}&=& n_{i'jk}+ n_{i'jk'}\\
        y_{i'j'k}+ y_{i'j'k'}&=& n_{i'j'k}+ n_{i'j'k'}\\
        y_{i'jk}+ y_{i'j'k}&=& n_{i'jk}+ n_{i'j'k}\\
        y_{i'jk'}+ y_{i'j'k'}&=& n_{i'jk'}+ n_{i'j'k'}\\
        y_{ijk} + y_{i'jk} &=& n_{ijk} + n_{i'jk}\\
        y_{ij'k} + y_{i'j'k} &=& n_{ij'k} + n_{i'j'k}\\
        y_{ijk'} + y_{i'jk'} &=& n_{ijk'} + n_{i'jk'}\\
        y_{ij'k'} + y_{i'j'k'} &=& n_{ij'k'} + n_{i'j'k'}\\
        \end{array}
        \right\}
        \]

 \item Randomly sample $u \in \mathbb Z_+$ uniformly from $[L, U]$ where
        \[
        L \gets \arg \min \{y_{ijk} \in \mathbb Z^8| y \in S((i, i'), (j, j'), (k, k'))\}
        \]
        \[
        U \gets \arg \max \{y_{ijk} \in \mathbb Z^8| y \in S((i, i'), (j, j'), (k, k'))\}
        \]
 \item Sample $u \in \mathbb Z_+$ from the hyper-geometric distribution with the probability 
\[
p(x) = \frac{\binom{U}{x}\binom{U}{U+L-x}}{\binom{2U}{L+U}}.
\]
 \item Let ${\bf z}$ be an $I\times J \times K$ contingency table with the $(i'', j'', k'')$th coordinate $z_{i'', j'', k''}$ for $i'' \in \{1, \ldots, I\}$, $j'' \in \{1, \ldots, J\}$, and $k'' \in \{1, \ldots, K\}$, such that 
    \[
    z_{i''j''k''} = \begin{cases}
    u & \mbox{if } i'' = i, \, j'' = j, \, k'' = k\\
    n_{ijk}+n_{ij'k}- u & \mbox{if } i'' = i, \, j'' = j', \, k'' = k\\
    n_{ijk}+n_{ijk'}-u & \mbox{if } i'' = i, \, j'' = j, \, k'' = k'\\
    n_{ijk}+n_{i'jk}-u & \mbox{if } i'' = i', \, j'' = j, \, k'' = k\\
    n_{ijk'}+n_{ij'k'}-z_{ijk'} & \mbox{if } i'' = i, \, j'' = j', \, k'' = k'\\
    n_{ij'k}+n_{i'j'k}-z_{ij'k} & \mbox{if } i'' = i', \, j'' = j', \, k'' = k\\
    n_{ijk'}+n_{i'jk'}-z_{ijk'} & \mbox{if } i'' = i', \, j'' = j, \, k'' = k'\\
    n_{ij'k'}+n_{i'j'k'}-z_{ij'k'} & \mbox{if } i'' = i', \, j'' = j', \, k'' = k'\\
    n_{i''j''k''} & \mbox{if } i'' \not \in \{i, i'\},  j'' \not \in \{j, j'\}, \mbox{ and }k'' \not \in \{k, k'\}. 
    \end{cases}
    \]
 \item Return ${\bf z}$\;
\end{enumerate}
\end{minipage}%
    }
    }}
    \hspace{8mm}
    \caption{Discrete Hit and Run (DHR) with Hyper-Geometric Distribution \cite{AD}}\label{alg:DHR4}
\end{algorithm}
\begin{adjustwidth}{0pt}{5pt}
Instead of using the hyper-geometric distribution, one can also use a uniform distribution at Step 4 in Algorithm \ref{alg:DHR4}.  

Now we outline three hybrid schemes: two for MCMC and DHR and one for MCMC and SIS to sample a table from $\mathcal{F}_{A, b}$.  Algorithms~\ref{alg:DHR} and~\ref{alg:SIS} use pre-defined count thresholds to call DHR or SIS, respectively. Algorithm~\ref{alg:DHR_DIS} uses pre-fixed intervals to call DHR. 
\end{adjustwidth}
   
\restoregeometry

\newgeometry{top=5mm}
\begin{algorithm}[H]
\centering
\fbox{\parbox{0.9\linewidth}{
\scalebox{.82}{
\begin{minipage}{.9\linewidth}
\KwData{\\
The observed table ${\bf N}_0$\\
$W>1$ --- desired number of iterations\\
$l \ge 1$ --- current iteration\\
${\bf n}_0$ --- sample point from fiber $\mathcal{F}_{A, b}$ (${\bf n}_0 = {\bf N}_0$ for $l = 1$)\\
$\beta > 1$ --- threshold to execute DHR\\
$ 0 \le c\le \beta$ --- current count from prior iterations\\
$\gamma>0$ --- threshold to execute DHR from ${\bf N}_0$
}
\KwResult{Sample $\bf z$ from the fiber $\mathcal{F}_{A, b}$, where $b$ is the vector of margins and $A$ is the design matrix for ${\bf n}_0$.}

$T \gets 1-\frac{l}{W}$\;
\eIf{$c < \beta$}{
Apply basic move to ${\bf n}_0$ to form ${\bf z}$\;
\eIf{$n_{ijk} \geq -1$ for all $n_{ijk} \in {\bf z}$}
{
\If{$n_{ijk} \ge 0$ for all $n_{ijk} \in {\bf z}$}{
Compute
        \[
        r \gets \frac{\prod_{{n_{ijk} \in {\bf z} }}  n_{ijk}!}{\prod_{n'_{ijk} \in {\bf n_0} } n'_{ijk}!}.
      \]
Set
        \[ 
        \bf n_1 \gets \begin{cases}
        \bf z & \text{with probability } \min\{r,1\}\\
        \bf n_0 & \text{Otherwise.}
        \end{cases} 
        \]
Set
        \[ 
        c \gets \begin{cases}
        0 & \bf {n_1}=\bf {z}\\
        $c + 1$ & \text{Otherwise.}
        \end{cases} 
        \]
}
}{$c \gets c + 1$\;
Define S where $(i,j,k) \in S$ such that $n_{ijk} \neq n'_{ijk}$\

Compute
            \[
            q \gets \sum_S n_{i,j,k}-\sum_S n'_{i,j,k}
            \]
Set
        \[ 
        \bf n_1 \gets \begin{cases}
        \bf z & \text{with probability } \min\{e^{\frac{q}{T}},1\}\\
        \bf n_0 & \text{Otherwise.}
        \end{cases} 
        \]
}

}
{$k\gets0$\\
Define ${\bf T_0}$ as I x J x K table such that $t_{ijk} = -1$  $\forall$ $t_{ijk}$ $\in $ ${\bf T_0}$

\While{$\exists$ $t_{ijk}=-1$}{
\eIf{$k<\gamma$}{Obtain ${\bf z}$ using Algorithm 2\

\eIf{$n_{ijk}\ge0$ $\forall$ $n_{ijk}$ $\in$ ${\bf z}$}{${\bf n_1} \gets {\bf z}$}{${\bf T_0} = {\bf z}$} }{Obtain ${\bf z}$ using Algorithm 2 with ${\bf N_0}$ as input\

\eIf{$n_{ijk}\ge0$ $\forall$ $n_{ijk}$ $\in$ ${\bf z}$\ 
}{${\bf n_1} \gets {\bf z}$}
{${\bf T_0} = {\bf z}$}
}
}$c\gets 0$ }
Return ${\bf n_1}$,c\;
\end{minipage}%
    }
    }}
    \hspace{4mm}
    \caption{Sampling with MCMC-DHR and Simulated Annealing}\label{alg:DHR}
\end{algorithm}
\restoregeometry




\newgeometry{top=25mm}
\begin{algorithm}[H]
\centering
\fbox{\parbox{0.9\linewidth}{
\scalebox{.82}{
\begin{minipage}{.9\linewidth}
\KwData{\\
The observed table ${\bf N}_0$\\
$W>1$ --- desired number of iterations\\
$l \ge 1$ --- current iteration\\
${\bf n}_0$ --- sample point from fiber $\mathcal{F}_{A, b}$ (${\bf n}_0 = {\bf N}_0$ for $l = 1$)\\
$\delta$ --- integer value indicating interval when to call DHR
}
\KwResult{Sample $\bf z$ from the fiber $\mathcal{F}_{A, b}$, where $b$ is the vector of margins and $A$ is the design matrix for ${\bf n}_0$.}

$T \gets 1-\frac{l}{W}$\;
\eIf{$(l\mod \delta) \ne 0$}{
Apply basic move to ${\bf n}_0$ to form ${\bf z}$\;
\eIf{$n_{ijk} \geq -1$ for all $n_{ijk} \in {\bf z}$}
{
\If{$n_{ijk} \ge 0$ for all $n_{ijk} \in {\bf z}$}{
Compute
        \[
        r \gets \frac{\prod_{{n_{ijk} \in {\bf z} }}  n_{ijk}!}{\prod_{n'_{ijk} \in {\bf n_0} } n'_{ijk}!}.
      \]
Set
        \[ 
        \bf n_1 \gets \begin{cases}
        \bf z & \text{with probability } \min\{r,1\}\\
        \bf n_0 & \text{Otherwise.}
        \end{cases} 
        \]
}
}{
Define $S$ where $(i,j,k) \in S$ such that $n_{ijk} \neq n'_{ijk}$\

Compute
            \[
            q \gets \sum_S n_{i,j,k}-\sum_S n'_{i,j,k}
            \]
Set
        \[ 
        \bf n_1 \gets \begin{cases}
        \bf z & \text{with probability } \min\{e^{\frac{q}{T}},1\}\\
        \bf n_0 & \text{Otherwise.}
        \end{cases} 
        \]
}

}
{\If{$n_{ijk} \ge 0$ $\forall$ \; $n_{ijk} \in {\bf n_0}$}
Define ${\bf T_0}$ as I x J x K table such that $t_{ijk} = -1$  $\forall$ $t_{ijk}$ $\in $ ${\bf T_0}$

\While{$\exists$ $t_{ijk}=-1$}{
Obtain ${\bf z}$ using Algorithm 2\

\eIf{$n_{ijk}\ge0$ $\forall$ $n_{ijk}$ $\in$ ${\bf z}$}{${\bf n_1} \gets {\bf z}$}{${\bf T_0} = {\bf z}$} 
}
}}
Return ${\bf n_1}$\;
\end{minipage}%
    }
    }}
   
    \caption{Sampling with MCMC-DHR and Simulated Annealing using pre-fixed DHR intervals}\label{alg:DHR_DIS}
\end{algorithm}

\restoregeometry



\begin{algorithm}[H]
\centering
\fbox{\parbox{.91\linewidth}{
\scalebox{.82}{
\begin{minipage}{.9\linewidth}
\KwData{\\
The observed table ${\bf N}_0$\\
$A$ --- $m \times n$ dimensional design matrix associated with ${\bf N}_0$\\
$b$ --- $m$ dimensional vector of row, column, and layer margins of ${\bf N}_0$\\
$W>1$ --- desired number of iterations\\
$l \ge 1$ --- current iteration\\
${\bf n}_0$ --- sample point from fiber $\mathcal{F}_{A, b}$ (${\bf n}_0 = {\bf N}_0$ for $l = 1$)\\
$\alpha > 1$ --- threshold to execute SIS\\
$ 0 \le c\le \alpha$ --- current count from prior iterations\\}
\KwResult{Sampled point $\bf z$ from the fiber $\mathcal{F}_{A, b}$, where $b$ is the vector of margins and $A$ is the design matrix for ${\bf n}_0$.}

$T \gets 1-\frac{l}{W}$\;
\eIf{$c < \beta$}{
Apply basic move to ${\bf n}_0$ to form ${\bf z}$\;
\eIf{$n_{ijk} \geq -1$ for all $n_{ijk} \in {\bf z}$}
{
\If{$n_{ijk} \ge 0$ for all $n_{ijk} \in {\bf z}$}{
Compute
        \[
        r \gets \frac{\prod_{{n_{ijk} \in {\bf z} }}  n_{ijk}!}{\prod_{n'_{ijk} \in {\bf n_0} } n'_{ijk}!}.
      \]
Set
        \[ 
        \bf n_1 \gets \begin{cases}
        \bf z & \text{with probability } \min\{r,1\}\\
        \bf n_0 & \text{Otherwise.}
        \end{cases} 
        \]
Set
        \[ 
        c \gets \begin{cases}
        0 & \bf {n_1}=\bf {z}\\
        $c + 1$ & \text{Otherwise.}
        \end{cases} 
        \]
}
}{$c \gets c + 1$\;
Define S where $(i,j,k) \in S$ such that $n_{ijk} \neq n'_{ijk}$\

Compute
            \[
            q \gets \sum_S n_{i,j,k}-\sum_S n'_{i,j,k}
            \]
Set
        \[ 
        \bf n_1 \gets \begin{cases}
        \bf z & \text{with probability } \min\{e^{\frac{q}{T}},1\}\\
        \bf n_0 & \text{Otherwise.}
        \end{cases} 
        \]
}

}
{

$D\gets FALSE$

\hspace{1.52mm}$c\gets0$

\While{$D=$ \textrm{FALSE}}{
Obtain ${\bf z}$ using SIS\\
\If{$n_{ijk}\ge0$ $\forall$ $n_{ijk}$ $\in$ ${\bf z}$}{${\bf n_1} \gets {\bf z}$

\hspace{1.1mm}$D\gets TRUE$}
}
}
Return ${\bf n_1}$,c\
\end{minipage}%
    }
    }}
    \hspace{4mm}
    \caption{Sampling with MCMC-SIS and Simulated Annealing}\label{alg:SIS}
\end{algorithm}

Each hybrid algorithm allows cell, $n_{ijk}$, to be $-1$ in order to connect tables throughout $\mathcal{F}_{A, b}$. However, any statistic taken from a sampled point outside the fiber, $\mathcal{F}_{A, b}$, is discarded.  These algorithms also call simulated annealing to increase the probability of sampling from inside the $\mathcal{F}_{A, b}$ as simulation iterations increase.

In Section~\ref{sec:computation} we apply each hybrid algorithm to simulated and real-world data for computational experiments.

\section{Computational Experiments}\label{sec:computation}

To demonstrate the effectiveness of hybrid Algorithms~\ref{alg:DHR},~\ref{alg:DHR_DIS}, and~\ref{alg:SIS} versus the basic MCMC approach, which is a MCMC algorithm with a set of basic moves (we define this MCMC scheme as the {\em basic MCMC} to discriminate from other MCMC schemes from here), we conduct computational experiments.  In each experiment we apply the hybrid algorithms in two ways: 1) allowing movement outside the fiber (i.e., allowing $n_{ijk}\ge-1$) and 2) only running a chain inside of the fiber itself.  In this paper, we focus on the first case.  The results of the second case are shown in the supplement.  

We use three data sets (two simulated and one empirical) for the experiments.  We describe these data sets and computational results below. 

\subsection{First Computational Experiment}\label{sec:simp}

For the initial experiment, we use simulated data on a $4 \times 4 \times 5$ table, $\bf N_0$, where its cell counts, $n_{ijk}$, are randomly chosen from $[0,7]$.  We describe the $4 \times 4 \times 5$ table below.


\begin{eqnarray*}
\begin{array}{|c|c|c|c|c|c|}\hline
0 & 3 & 5 & 6\\\hline 
0 & 3 & 0 & 1\\\hline
7 & 0 & 4 & 0\\\hline
1 & 0 & 6 & 0\\\hline
\end{array} \, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
5 & 7 & 0 & 2\\\hline 
0 & 6 & 2 & 1\\\hline
2 & 5 & 2 & 6\\\hline
0 & 6 & 4 & 3\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
5 & 6 & 5 & 1\\\hline 
5 & 6 & 1 & 4\\\hline
0 & 5 & 2 & 6\\\hline
5 & 6 & 4 & 1\\\hline
\end{array}\\
\end{eqnarray*}
\begin{eqnarray*}
\begin{array}{|c|c|c|c|c|c|}\hline
6 & 2 & 3 & 3\\\hline 
5 & 3 & 6 & 3\\\hline
5 & 7 & 5 & 7\\\hline
4 & 4 & 1 & 3\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
3 & 0 & 4 & 4\\\hline 
0 & 3 & 0 & 5\\\hline
0 & 5 & 1 & 6\\\hline
3 & 6 & 4 & 5\\\hline
\end{array}\\
\end{eqnarray*}

We conduct this first computational experiment with four methods: the MCMC with a set of basic moves,  Algorithm~\ref{alg:DHR}, Algorithm~\ref{alg:DHR_DIS}, and Algorithm~\ref{alg:SIS}. We iterate each chain $1.5$ million iterations.

We evaluate each result using the un-normalized log likelihood (UNLL) as shown by Kahle, Yoshida, and Garcia-Puente \cite{KYG}.  This is a modification of Fisher's exact test defined on a table with $r$ cells:

\begin{align*}
    f(\mathbf{t|(A, b)}) &= \frac{\frac{1}{n_1!...n_r!}}{\sum_{\mathcal{X}\;\in\mathcal{F}_{A,b}}{\frac{1}{x_1!...X_r!}}} \; \text{for }\mathbf{n} \in \mathcal{F}_{A, b}.
\end{align*}

Because we cannot compute its normalizing constant, we instead use the un-normalized likelihood function:

\begin{align*}
    \tilde{f}(\mathbf{n}|\mathbf{(A, b)}) = \frac{1}{n_1!...n_r!}\; \text{for }\mathbf{n} \in \mathcal{F}_{A, b}.
\end{align*}

Then we take logarithm on the unnormalized likelihood function, which we call the {\em UNLL}.  While this is tractable for many of our purposes, {\tt R} cannot handle the UNLL with very large cell counts. 
Thus, we also use Pearson $\chi^2$ statistics calculated from tables sampled from the fiber.  Generally, we observe that we arrive at similar conclusions using either metric.  For the UNLL, larger values suggest sampled tables are closer to the expected table under the no-three-way interaction model while the same may be inferred for small $\chi^2$ statistics.

Using the simulated data defined above, we compare the two metrics using Algorithm~\ref{alg:DHR} with a $\beta=1,000$.  We also include the result from the MCMC scheme with the set of basic moves, the basic MCMC, from the same data.

Figure~\ref{fig:MCMC_Pear} shows the progression of UNLL values (left) and $\chi^2$ statistics (right) for basic MCMC and MCMC-DHR using Algorithm~\ref{alg:DHR} on the simulated data.   
\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{SIM_UNLL_1000.pdf}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{SIM_DHR_1000.pdf}
  
 \end{minipage}
 \caption{\centering UNLL values (left) and $\chi^2$ statistics (right) for simulated data using basic MCMC (red) and MCMC-DHR (blue).}
 \label{fig:MCMC_Pear}
\end{figure}

Results suggest that either metric derived from tables sampled from the fiber progressively converge asymptotically using Algorithm~\ref{alg:DHR}.
For the rest of this section focuses on UNLL vice $\chi^2$ statistics.

\subsubsection{Experimenting with Simulated Data}

First, we apply Algorithm~\ref{alg:DHR} to the simulated data above by assigning different values to $\beta$  since $\beta$ is a tuning parameter where $n_{ijk} \ge -1$.  From our computational experiments, if we decrease (increase) $\beta$ then Algorithm~\ref{alg:DHR} turns on the DHR more (less) often.  Figure~\ref{fig:BASEMCMC} shows the estimated distribution of $\chi^2$ statistics via the basic MCMC plotted against the asymptotic distribution with $ (I-1)(J-1)(K-1)$ degrees of freedom which is the degrees of freedom used for the no-three-way-interaction model.  We use this plot for comparison against the hybrid methods.

\begin{figure}[H]
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{BASEMCMC.pdf}
  \caption{Basic MCMC on simulated data.}
  \label{fig:BASEMCMC}
\end{figure}
 













   

With Algorithm~\ref{alg:DHR} with this $\beta$, we observe often simulated annealing in a chain because we wish to keep a chain in the interior of $\mathcal{F}_{A, b}$ when we allow cell counts to be $-1$.
Results are shown in Figures~\ref{fig:MCMC_DHR_4000},~\ref{fig:MCMC_DHR_2500},~\ref{fig:MCMC_DHR_1000},~\ref{fig:MCMC_DHR_500}, and~\ref{fig:MCMC_DHR_250} depending on the tuning parameter $\beta$. 



\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{MCM-DHR_NEG_4000.pdf}
  \caption{UNLL with $\beta=4,000$.}
  \label{fig:MCMC_DHR_4000}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{MCMC-DHR_NEG_2500.pdf}
  \caption{UNLL with $\beta=2,500$.}
  \label{fig:MCMC_DHR_2500}
 \end{minipage}
\end{figure}

\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{MCMD-DHR_NEG_1000.pdf}
  \caption{UNLL with $\beta=1,000$.}
  \label{fig:MCMC_DHR_1000}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{MCMC-DHR_NEG_500.pdf}
  \caption{UNLL with $\beta=500$.}
  \label{fig:MCMC_DHR_500}
 \end{minipage}
\end{figure}

\begin{figure}[H]
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{MCMC_DHR_NEG_250.pdf}
  \caption{UNLL with $\beta=250$.}
  \label{fig:MCMC_DHR_250}
\end{figure}

Regardless of the threshold used, the results show each chain converges to approximately the same values.  This is a consequence of simulated annealing.  As iterations increase, we observe DHR less often because ``cooling'' process forces a chain inside of $\mathcal{F}_{A, b}$.  For results shown in Figure~\ref{fig:MCMC_DHR_4000}, we never observe DHR with Algorithm~\ref{alg:DHR}.  At $\beta = 4,000$, Algorithm~\ref{alg:DHR} exclusively uses simulated annealing to move inside of $\mathcal{F}_{A, b}$.

For lower thresholds, we observe more fluctuation in the UNLL and $\chi^2$ statistics since with lower thresholds DHR occurs more often.  

Figures~\ref{fig:DIST_DHR_4000},~\ref{fig:DIST_DHR_2500A},~\ref{fig:DIST_DHR_1000},~\ref{fig:DIST_DHR_500}, and~\ref{fig:DIST_DHR_250} show estimated distributions using MCMC-DHR with varying $\beta$. Table~\ref{tab:DHR_NEG} shows the number of DHR called for each $\beta$. 

\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_NEG_DHR4000.pdf}
  \caption{\centering MCMC-DHR with $\beta=4,000$.}
  \label{fig:DIST_DHR_4000}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.8cm,height=6.9cm]{DIST_NEG_DHR2500.pdf}
  \caption{\centering MCMC-DHR with $\beta=2,500$.}
  \label{fig:DIST_DHR_2500A}
 \end{minipage}
\end{figure}

\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_NEG_DHR1000.pdf}
  \caption{\centering MCMC-DHR with $\beta=1,000$.}
  \label{fig:DIST_DHR_1000}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{DIST_NEG_DHR500.pdf}
  \caption{\centering MCMC-DHR with $\beta=500$.}
  \label{fig:DIST_DHR_500}
 \end{minipage}
\end{figure}

\begin{figure}[H]
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_NEG_DHR250.pdf}
  \caption{\centering MCMC-DHR with $\beta=250$.}
  \label{fig:DIST_DHR_250}
\end{figure}


\begin{table}[H]
\caption{Number of DHR excursions with varying thresholds.}
\begin{center}
\begin{tabular}{|c|c|}
    \hline
    $\beta$&DHR Iterations\\\hline
    4,000&0\\ \hline
    2,500&0 \\ \hline
    1,000&1\\\hline
    500&90 \\ \hline
    250&597\\\hline
\end{tabular}
\end{center}
\label{tab:DHR_NEG}
\end{table}

An estimated distribution of $\chi^2$ statistics or UNLL with Algorithm~\ref{alg:DHR} depends on the tuning parameter, $\beta$.  From our initial experiments, we recommend that we should use small $\beta$ for DHR for larger cell counts in Algorithm~\ref{alg:DHR}.  



\subsubsection{Experimenting with DHR for Pre-fixed Intervals}

The tuning parameter, $\delta$, in Algorithm~\ref{alg:DHR_DIS} is a user-defined, pre-fixed interval. In order to see how different intervals affect running a Markov chain we also vary $\delta$.   

Figure~\ref{fig:DHR_PRE} shows how UNLL changes over the iteration for the simulated data set with allowing cell counts to be $-1$. The pre-fixed interval is set intentionally low so that we have DHR in each chain.  For this experiment we have DHR in every 10 iterations.  

\begin{figure}[H]
  \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DHR_PRE_10.pdf}
 \caption{\centering MCMC-DHR with $\delta=10$ $n_{ijk} \ge -1$.}
 \label{fig:DHR_PRE}
\end{figure}

Figure~\ref{fig:DIST_DHR_PRE} shows the estimated distribution of $\chi^2$ statistics using Algorithm~\ref{alg:DHR_DIS} allowing cell counts to be $-1$.   With Algorithm~\ref{alg:DHR_DIS} we have an estimated distribution which is close to the asymptotic distribution of $\chi^2$ statistics under the no-three-way interaction model.

\begin{figure}[H]
  \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_PRE_10.pdf}
 \caption{\centering MCMC-DHR with $\delta=10$ for $n_{ijk} \ge -1$.}
 \label{fig:DIST_DHR_PRE}
\end{figure}

\subsubsection{Experimenting with Random Starts}

In this section, we also consider random starts for multiple chains via Algorithm~\ref{alg:DHR} on the simulated data.   Random starts allow for re-running the algorithm multiple times using the final sampled table of the prior run as a starting point.  Ideally, this method provides more opportunity for sampling the interior of $\mathcal{F}_{A, b}$ without bias.  In this experiment, we use a single threshold, $\beta=750$, for each of seven repetitions of the algorithm.  Figure~\ref{fig:RESTART_DHR_750} shows how UNLL changes over the iteration.

\begin{figure}[H]
    \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DHR_UNLL.pdf}
  \caption{\centering MCMC-DHR with random starts and $\beta=750$.}
  \label{fig:RESTART_DHR_750}

\end{figure}

The method significantly increases the number of points inside  $\mathcal{F}_{A,b}$ and each restart (which can be identified by the large negative value spikes) converges to approximately the same values.
From Figure~\ref{fig:RESTART_DHR_750}, Algorithm~\ref{alg:DHR} samples points from broader locations in $\mathcal{F}_{A, b}$.

Figure~\ref{fig:RESTART_DIST_750} shows the estimated distribution of $\chi^2$ statistics using different restarts. An estimated distribution of $\chi^2$ statistics differs noticeably from the asymptotic distribution.  We believe this happens due to the combination of chains being outside of $\mathcal{F}_{A, b}$ by allowing cell counts to be $-1$ and the effect of DHR.


\begin{figure}[H]
    \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_REST.pdf}
  \caption{\centering MCMC-DHR with random starts and $\beta=750$.}
  \label{fig:RESTART_DIST_750}

\end{figure}

\subsubsection{Experimenting with MCMC and SIS}

Similar to Algorithm~\ref{alg:DHR} where $\beta$ is a user-defined threshold parameter, we have a threshold tuning parameter, $\alpha$ in Algorithm~\ref{alg:SIS} to determine when to execute the SIS portion of the algorithm.  As $\alpha$ increases, SIS is called less.  Figures~\ref{fig:MCMC_SIS_1000},~\ref{fig:MCMC_SIS_500},~\ref{fig:MCMC_SIS_250},~\ref{fig:MCMC_SIS_100}, and~\ref{fig:MCMC_SIS_50} show how UNLL changes for different values of $\alpha$.

\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{MCMC_SIS_1000.pdf}
  \caption{UNLL for $\alpha = 1,000$.}
  \label{fig:MCMC_SIS_1000}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{MCMC_SIS_500.pdf}
  \caption{UNLL for $\alpha = 500$.}
  \label{fig:MCMC_SIS_500}
 \end{minipage}
\end{figure}

\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{MCMC_SIS_250.pdf}
  \caption{UNLL for $\alpha = 250$.}
  \label{fig:MCMC_SIS_250}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{MCMC_SIS_100.pdf}
  \caption{UNLL for $\alpha = 100$.}
  \label{fig:MCMC_SIS_100}
 \end{minipage}
\end{figure}

\begin{figure}[H]
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{MCMC_SIS_50.pdf}
  \caption{UNLL for $\alpha = 50$.}
  \label{fig:MCMC_SIS_50}
\end{figure}

As previously discussed, SIS is called more with smaller values of $\alpha$ and it causes a chain to move further away from the expected table.  As iterations increase, simulated annealing will have a greater impact on movement in $\mathcal{F}_{A, b}$.  However, this is affected by $\alpha$: the more often SIS is called, the more often we sample from $\mathcal{F}_{A, b}$ further away from the expected table.  

Figures~\ref{fig:DIST_SIS_1000},~\ref{fig:DIST_SIS_500},~\ref{fig:DIST_SIS_250},~\ref{fig:DIST_SIS_100}, and~\ref{fig:DIST_SIS_50} show estimated distributions of $\chi^2$ statistics using Algorithm~\ref{alg:SIS}.  Like MCMC-DHR in Algorithm~\ref{alg:DHR}, the lower the tuning parameter the more skewed the distribution becomes.  While this is also true for MCMC-DHR, it is much less pronounced with varying values of $\beta$ in Algorithm~\ref{alg:DHR}. 

\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_SIS_1000.pdf}
  \caption{\centering MCMC-SIS with $\alpha = 1,000$.}
  \label{fig:DIST_SIS_1000}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_SIS_500.pdf}
  \caption{\centering MCMC-SIS with $\alpha = 500$.}
  \label{fig:DIST_SIS_500}
 \end{minipage}
\end{figure}

\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_SIS_250.pdf}
  \caption{\centering MCMC-SIS with $\alpha = 250$.}
  \label{fig:DIST_SIS_250}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{DIST_SIS_100.pdf}
  \caption{\centering MCMC-SIS with $\alpha = 100$.}
  \label{fig:DIST_SIS_100}
 \end{minipage}
\end{figure}

\begin{figure}[H]
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_SIS_50.pdf}
  \caption{\centering MCMC-SIS with $\alpha = 50$.}
  \label{fig:DIST_SIS_50}
\end{figure}

Table~\ref{tab:SIS_RESb} shows the number of times SIS is called with different values of $\alpha$.  As $\alpha$ decreases, the number of times SIS is called increases. We note that when $\alpha$ is small enough, relatively small changes to the tuning parameter results in significant changes to how many times SIS is used (e.g., $\alpha=50$ to $\alpha=100$ in the simulated data).  However, how small $\alpha$ needs to be to see significant change is dependent on the size of the contingency table and the cell counts.

\begin{table}[H]
\caption{Number of SIS excursions with varying thresholds.}
\begin{center}
\begin{tabular}{|c|c|}
    \hline
    $\alpha$&SIS Iterations\\\hline
    1,000&6\\ \hline
    500&82 \\ \hline
    250&612\\\hline
    100&5,284\\\hline
    50&18,573\\\hline
\end{tabular}
\end{center}
\label{tab:SIS_RESb}
\end{table}

Experimenting with random restarts in conjunction with MCMC-SIS was conducted similarly to MCMC-DHR with random restarts above.  We use five restarts with the tuning parameter is set to $\alpha=100$. Figures~\ref{fig:SIS_RESTAN} and~\ref{fig:SIS_RESTAN_DIST} show random restart results incorporated with MCMC-SIS.

\begin{figure}[H]
 \centering
  \includegraphics[width=6.8cm,height=6.6cm]{SIS_SP_RESTAN.pdf}
  \caption{\centering MCMC-SIS with random starts.}
  \label{fig:SIS_RESTAN}
 \end{figure}

\begin{figure}[H]
 \centering
  \includegraphics[width=6.8cm,height=6.6cm]{SIS_DIST_RESTAN.pdf}
  \caption{\centering MCMC-SIS with random starts.}
  \label{fig:SIS_RESTAN_DIST}
 \end{figure}

We observe little difference between a single iteration or using restarts for SIS.  The distribution estimate is approximately identical.

In order to compare we use the hybrid scheme of SIS and MCMC with a set of basic moves described in Algorithm~\ref{alg:SIS} introduced by Kahle et al.~\cite{KYG}.  This means that because we use random starting tables sampled by SIS, we can sample via multiple chains from $\mathcal{F}_{A, b}$ so we hope to solve the problem on connectivity using a set of basic moves.  This means that we do not allow cell counts to be $-1$ in this experiments.  
Figures~\ref{fig:SIS_15_REST} show results from MCMC-SIS with $\alpha=15$ for a single iteration (left) and with random starts (right). 

\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_SIS_15.pdf}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_SIS_15_REST.pdf}
  \end{minipage}
  

 \caption{\centering MCMC-SIS with cell counts $n_{ijk} \ge0$ and with $\alpha = 15$ for a single chain (left) and multiple chains with random starts (right).}\label{fig:SIS_15_REST}
\end{figure}

\begin{comment}
\begin{figure}[H]
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_SIS_15_ex.pdf}
  \end{minipage}\hfill
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_SIS_15_REST_ex.pdf}
 \end{minipage}
 \caption{\centering Estimated distribution of $\chi^2$ statistics captured using MCMC-SIS with cell counts $n_{ijk} \ge0$ and with $\alpha = 15$ for a single chain (left) and multiple chains with random starts (right).}\label{fig:SIS_15_REST_ex}
\end{figure}
\end{comment}

Results show that an estimated distribution of $\chi^2$ statistics is close to the asymptotic distribution of $\chi^2$ statistics.

\subsection{Experimenting with Sparse Data}

Because it is challenging to run a Markov chain with sparse contingency tables, we also examine how Algorithm~\ref{alg:DHR},~\ref{alg:DHR_DIS}, and~\ref{alg:SIS} perform on sparse data. 

First, we use Algorithms~\ref{alg:DHR} in two ways: 1) single iteration and 2) with restarts. The data is illustrated below.

\begin{eqnarray*}
\begin{array}{|c|c|c|c|c|c|}\hline
6 & 2 & 3 & 3\\\hline 
5 & 3 & 6 & 3\\\hline
5 & 7 & 5 & 7\\\hline
4 & 4 & 1 & 3\\\hline
\end{array}\, \, 
\begin{array}{|c|c|c|c|c|c|}\hline
4 & 0 & 0 & 0\\\hline 
0 & 4 & 0 & 0\\\hline
0 & 0 & 0 & 4\\\hline
0 & 0 & 4 & 0\\\hline
\end{array}\\
\end{eqnarray*}

In this experiment, we set $I=4$, $J=4$, and $K=2$. We intentionally use small cell counts in the second dimension of the $K$ direction to simulate significant difference between the $K$ categories.  Contingency tables of smaller $I,\, J,\, K$ require smaller $\beta$ threshold values in order to consistently call the DHR portion of Algorithm~\ref{alg:DHR}. For this experiment, we use a threshold of 50.

Figures~\ref{fig:UNLL_DHR_SP_50},~\ref{fig:UNLL_DHR_SP}, and~\ref{fig:UNLL_MCMC_SP} show how UNLL changes for each iteration of MCMC-DHR with thresholds of 50 and 500, respectively, and with basic MCMC.  

\begin{figure}[H]
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{DHR_SP_50.pdf}
  \caption{\centering MCMC-DHR with $\beta=50$.}
  \label{fig:UNLL_DHR_SP_50}
 \end{minipage}\hfill
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.9cm,height=6.6cm]{DHR_SP.pdf}
  \caption{\centering MCMC-DHR with $\beta=500$.}
  \label{fig:UNLL_DHR_SP}
  \end{minipage}
\end{figure}

\begin{figure}[H]
 \centering
  \includegraphics[width=6.8cm,height=6.6cm]{MCMC_SP.pdf}
  \caption{\centering UNLL for MCMC.}
  \label{fig:UNLL_MCMC_SP}
 \end{figure}

There is little difference between the three methods with respect to UNLL.  We can see differences in the larger UNLL associated with the MCMC-DHR iterations at the beginning of each run.  Despite the differences, each iteration converges to approximately the same values.

Figures~\ref{fig:DHR_SP_50} and ~\ref{fig:DHR_SP_500} show estimates distributions of $\chi^2$ statistics using different values of $\beta$.  Figure~\ref{fig:MCMC_SPa} shows the estimated distribution using the basic MCMC algorithm applied to the sparse data.

\begin{figure}[H]
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{DHR_SP_DIST.pdf}
  \caption{\centering MCMC-DHR with $\beta=50$.}
  \label{fig:DHR_SP_50}
 \end{minipage}\hfill
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.9cm,height=6.6cm]{DHR_DIST_HSP.pdf}
  \caption{\centering MCMC-DHR with $\beta=500$.}
  \label{fig:DHR_SP_500}
  \end{minipage}
\end{figure}

\begin{figure}[H]
 \centering
  \includegraphics[width=6.8cm,height=6.6cm]{MCMC_SP_DIST.pdf}
  \caption{\centering Basic MCMC on sparse data.}
  \label{fig:MCMC_SPa}
 \end{figure}

With random restarts, we use Algorithm~\ref{alg:DHR} to sample from the interior of $\mathcal{F}_{A, b}$.  Figures~\ref{fig:MCMC_SPb} and~\ref{fig:MCMC_SP} show how UNLL changes over each iteration and an estimated distribution of $\chi^2$ statistics, respectively.


\begin{figure}[H]
 \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DHR_SP_REST.pdf}
  \caption{\centering MCMC-DHR with random starts.}
  \label{fig:MCMC_SPb}
 \end{figure}

\begin{figure}[H]
 \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DIST_DHR_REST_SP.pdf}
  \caption{\centering MCMC-DHR with random starts.}
  \label{fig:MCMC_SP}
 \end{figure}
 
 Figure~\ref{fig:PRE_SP_DIST} shows an estimated distribution of $\chi^2$ statistics with Algorithm~\ref{alg:DHR_DIS}. In this case, the DHR is applied every fiftieth iteration.
 
 \begin{figure}[H]
    \centering
  \includegraphics[width=6.8cm,height=6.6cm]{DHR_PRE_SP.pdf}
  \caption{\centering MCMC-DHR with $\delta=50$.}
  \label{fig:PRE_SP_DIST}

\end{figure}

Lastly, we use Algorithm~\ref{alg:SIS} on the sparse data as a single iteration and with random restarts.  We use in two separate iterations where  $\alpha = 25$ and $\alpha = 250$.  For sparse data, we recommend to use smaller $\alpha$ to call SIS in Algorithm~\ref{alg:SIS}.
With $\alpha=250$, only simulated annealing affects sampling, similar to MCMC-DHR with larger values of $\beta$.  Figures~\ref{fig:UNLL_SIS_SP_25} and~\ref{fig:UNLL_SIS_SP_250} illustrate how UNLL changes over each iteration. 

\begin{figure}[H]
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{MCMC_SIS_25sp.pdf}
  \caption{\centering MCMC-SIS with $\alpha=25$.}
  \label{fig:UNLL_SIS_SP_25}
 \end{minipage}\hfill
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.9cm,height=6.6cm]{MCMC_SIS_250sp.pdf}
  \caption{\centering MCMC-SIS with $\alpha=250$.}
  \label{fig:UNLL_SIS_SP_250}
  \end{minipage}
\end{figure}

We observe that with $\alpha=25$, SIS is called throughout the duration of the iteration.  When $\alpha=250$, we have similar results shown in Figures~\ref{fig:UNLL_DHR_SP_50} and~\ref{fig:UNLL_DHR_SP}.  Notably, with the MCMC-DHR algorithm with  $\beta=50$ UNLL fluctuates less compared with MCMC-SIS with a value of $\alpha=25$.  While in either case, DHR and SIS, respectively, are called more often due to their tuning parameters, MCMC-DHR results suggest that DHR sends a chain to tables far from the expected table.

Figures~\ref{fig:DIST_SIS_SP_25} and~\ref{fig:DIST_SIS_SP_250} show estimated distributions with MCMC-SIS on sparse data and 
iterations with $\alpha=25$ result in a long right tail.

\begin{figure}[H]
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{DIST_SIS_25sp.pdf}
  \caption{\centering MCMC-SIS with $\alpha=25$.}
  \label{fig:DIST_SIS_SP_25}
 \end{minipage}\hfill
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.9cm,height=6.6cm]{DIST_SIS_250sp.pdf}
  \caption{\centering MCMC-SIS with $\alpha=250$.}
  \label{fig:DIST_SIS_SP_250}
  \end{minipage}
\end{figure}

We apply Algorithm~\ref{alg:SIS} to use random restarts where the starting table is the final sampled table of a prior iteration.  We use a single $\alpha$ and conduct the experiment in five replications where each replication iterates 1.5 million times.  Figures~\ref{fig:SIS_SP_RESTA} and~\ref{fig:SIS_SP_REST} show how UNLL and $\chi^2$ statistics change over each iteration.   

\begin{figure}[H]
 \centering
  \includegraphics[width=6.8cm,height=6.6cm]{SIS_SP_REST.pdf}
  \caption{\centering MCMC-SIS with random starts.}
  \label{fig:SIS_SP_RESTA}
 \end{figure}

\begin{figure}[H]
 \centering
  \includegraphics[width=6.8cm,height=6.6cm]{SIS_DIST_SP_REST.pdf}
  \caption{\centering MCMC-SIS with random starts.}
  \label{fig:SIS_SP_REST}
 \end{figure}

Sparse simulated data better represents data we find in the real-world.  Experimental results suggest that with MCMC-DHR using Algorithm~\ref{alg:DHR}, regardless of $\beta$ value,  provides distribution estimates similar to those of the basic MCMC.  This is not the case with MCMC-DHR in Algorithm~\ref{alg:DHR}.
Like experiments with Algorithm~\ref{alg:DHR} on sparse data, similar inference can be made regarding the use of MCMC-SIS in Algorithm~\ref{alg:SIS}. However, if we use the hybrid of MCMC and SIS with random starts, we might be able to guarantee the connectivity of a Markov chain while whereas the basic MCMC might not.  Generally, all methods performed similarly with the exception of interval-based MCMC-DHR (i.e., Algorithm~\ref{alg:DHR_DIS}.   

\subsection{USN Officer and Enlisted Population by Race, Rank, and Gender}


We now examine the performance of the algorithms using an empirical data set from \cite{NavyData}.  The data set consists of the USN officer and enlisted population as of January 20, 2021.  The data is organized into a three way contingency table stratified by race, rank, and gender.  The rank category consists of 19 levels.  The race category consists of six levels and gender is binary.  This results in a $19 \times 6 \times 2$ contingency table.  The values are displayed in Tables~\ref{tab:maloff},~\ref{tab:malenl},~\ref{tab:femoff}, and~\ref{tab:femenl}.  Note that because of the magnitude of many of the cell counts, we cannot use UNLL.

\begin{table}[H]
\caption{Male Naval Officer Population by Race and Rank as of January 20,2021}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
    \hline
    Race&Adm.&O-6&O-5&O-4&O-3&O-2&O-1&W-4&W-3&W-2\\\hline
    Nat. Am.&1&9&29&92&154&64&66&7&14&15\\ \hline
    Asian&2&96&225&417&832&347&349&12&51&34 \\ \hline
    Af. Am.&6&167&357&551&1,006&329&396&81&148&142\\\hline
    Pac. Isl.&2&107&263&280&413&127&118&14&38&50\\\hline
    Multi-Race&2&36&123&242&853&311&342&13&19&14\\\cline{1-11}
    White&192&2,452&4,752&6,517&11,635&4,038&4,038&222&413&371\\
    \hline
\end{tabular}
\end{center}
\label{tab:maloff}
\end{table}

\begin{table}[H]
\caption{Male Enlisted Population by Race and Rank as of January 20,2021}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
    \hline
    Race&E-9&E-8&E-7&E-6&E-5&E-4&E-3&E-2&E-1\\\hline
    Nat. Am&44&210&777&1,863&1,142&379&254&110&56\\ \hline
    Asian&129&466&1,257&2,720&3,405&2,229&2,051&693&452 \\ \hline
    Af. Am.&409&994&3,151&7,514&9,963&6,533&5,968&2,820&2,036\\\hline
    Pac. Isl.&93&308&798&1,387&2,593&3,087&1,922&333&30\\\hline
    Multi-Race&48&162&942&4,679&5,004&2,189&1,279&476&348\\\hline
    White&1,809&4,242&11,591&26,435&34,716&25,716&20,871&9,140&6,502\\
    \hline
\end{tabular}
\end{center}
\label{tab:malenl}
\end{table}

\begin{table}[H]
\caption{Female Naval Officer Population by Race and Rank as of January 20,2021}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
    \hline
    Race&Adm.&O-6&O-5&O-4&O-3&O-2&O-1&W-4&W-3&W-2\\\hline
    Nat. Am.&0&2&11&12&45&27&32&0&0&0\\ \hline
    Asian&1&28&56&136&323&129&126&0&3&4 \\ \hline
    Af. Am.&0&50&103&226&466&143&141&13&29&38\\\hline
    Pac. Isl.&1&14&39&81&215&55&45&0&8&6\\\hline
    Multi-Race&0&5&33&87&327&130&144&1&3&1\\\hline
    White&13&294&677&1,447&2,955&1,089&1,117&12&21&33\\
    \hline
\end{tabular}
\end{center}
\label{tab:femoff}
\end{table}

\begin{table}[H]
\caption{Female Enlisted Population by Race and Rank as of January 20,2021}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
    \hline
    Race&E-9&E-8&E-7&E-6&E-5&E-4&E-3&E-2&E-1\\\hline
    Nat. Am&6&24&114&294&314&173&124&53&23\\ \hline
    Asian&14&71&214&545&863&652&677&198&129 \\ \hline
    Af. Am.&88&915&3,151&2,509&4,453&3,039&2,983&1,050&894\\\hline
    Pac. Isl.&13&50&157&297&779&920&714&132&5\\\hline
    Multi-Race&5&24&181&975&1,536&757&544&165&122\\\hline
    White&120&324&1,267&3,536&7,555&6,669&5,914&2,308&1,748\\
    \hline
\end{tabular}
\end{center}
\label{tab:femenl}
\end{table}


Of note, the table consists of a wide range of cell counts.  In general, cell counts decrease as rank increases which is expected as there are fewer billets for higher ranks thus reflecting the USN's ``up or out'' promotion structure.  Further, we note that there are considerably fewer females than males across all ranks and races.  For this reason, several sub-categories were combined to ensure that the chain runs smoothly. 

Figure~\ref{fig:MCMC_DHR_NAVY} shows results after calling Algorithm~\ref{alg:DHR} on the data with varying $\beta$ and allowing movement outside the fiber.  

\begin{figure}[H]
  \centering
  \includegraphics[width=10cm,height=8cm]{NAV_BETA.pdf}
   \caption{Results using Algorithm~\ref{alg:DHR} (MCMC-DHR) and basic MCMC on USN data with varying thresholds}
   \label{fig:MCMC_DHR_NAVY}
    \end{figure}

Note that the basic MCMC provides much smaller Pearson $\chi^2$ statistics than any variation of the MCMC-DHR model. 
These results suggest that when DHR is called more often, a chain moves widely across $\mathcal{F}_{A, b}$ yielding a wide range of $\chi^2$ statistics.  As $\beta$ increases, it stays interior of $\mathcal{F}_{A, b}$.  This is partly a result of simulated annealing since, as the iteration increases, the ``cooling'' temperature decreases, resulting in the threshold count being reached fewer times.  Table~\ref{tab:DHR_RES} illustrates the number of DHR called in Algorithm~\ref{alg:DHR}.

\begin{table}[H]
\caption{Number of DHR called with varying thresholds.}
\begin{center}
\begin{tabular}{|c|c|}
    \hline
    $\beta$&DHR Iterations\\\hline
    500&157\\ \hline
    1,200&1 \\ \hline
    2,000&0\\\hline
\end{tabular}
\end{center}
\label{tab:DHR_RES}
\end{table}

Note that in the MCMC-DHR method using a threshold of 2,000 DHR is not called in the algorithm.  As with the simulated data, this is equivalent to exclusively using simulated annealing in a method.  Figures~\ref{fig:DHR_NAV_2000} and~\ref{fig:DHR_NAV_1200} show an estimated distributions of $\chi^2$ statistics for threshold of 2,000 and 1,200, respectively.

\begin{figure}[H]
  \begin{minipage}[c]{.5\linewidth}
  \centering
  \includegraphics[width=6.9cm,height=6.6cm]{NAV_DHR_DIST2K.pdf}
  \caption{\centering MCMC-DHR with $\beta=2,000$.}
  \label{fig:DHR_NAV_2000}
 \end{minipage}\hfill
 \begin{minipage}[c]{.5\linewidth}
\centering
  \includegraphics[width=6.9cm,height=6.6cm]{NAV_DHR_DIST12.pdf}
  \caption{\centering MCMC-DHR with $\beta=1,200$.}
  \label{fig:DHR_NAV_1200}
  \end{minipage}
\end{figure}

We also use Algorithm~\ref{alg:DHR_DIS} on the USN data.  From the previous example with Algorithm~\ref{alg:DHR}, we expect differences in the result depending on the chosen interval.  Figure~\ref{fig:MCMC_DHR_NAVYa} shows results with three different intervals and the basic MCMC method.

\begin{figure}[H]
  \centering
  \includegraphics[width=10cm,height=8cm]{NAV_PRE_DET.pdf}
   \caption{Results using MCMC-DHR (and basic MCMC) on USN data with pre-fixed intervals.}
   \label{fig:MCMC_DHR_NAVYa}
    \end{figure}

The more often DHR is called (i.e., in smaller intervals), the more iterations to converge.  For an interval of 10, the algorithm does not have a chance to consistently call the basic MCMC or simulated annealing portions of Algorithm~\ref{alg:DHR_DIS} resulting in large $\chi^2$ statistics.

Finally, Algorithm~\ref{alg:SIS} is applied on the USN personnel data.  Table~\ref{tab:SIS_RES} shows the number of times SIS is called for three values $\alpha$. Figure~\ref{fig:MCMC_SIS_NAVYa} shows the results.  Compared with MCMC-DHR results shown in Figure~\ref{fig:MCMC_DHR_NAVY}, $\chi^2$ statistics in Figure~\ref{fig:MCMC_SIS_NAVYa} suggest that when SIS is called, it takes longer to move back towards the expected table.  By comparison, even though MCMC-DHR sampling in Algorithm~\ref{alg:DHR} results in large $\chi^2$ statistics, the chain moves towards the expected table more readily than the MCMC-SIS algorithm.



\begin{table}[H]
\caption{Number of SIS excursions with varying thresholds.}
\begin{center}
\begin{tabular}{|c|c|}
    \hline
    $\alpha$&SIS Iterations\\\hline
    750&131\\ \hline
    1,000&46 \\ \hline
    2,500&0\\\hline
\end{tabular}
\end{center}
\label{tab:SIS_RES}
\end{table}

\begin{figure}[H]
  \centering
  \includegraphics[width=10cm,height=8cm]{NAV_ALPHA.pdf}
   \caption{MCMC-SIS results on USN data with different values of $\alpha$.}
   \label{fig:MCMC_SIS_NAVYa}
    \end{figure}

Figure~\ref{fig:MCMC_DIST_NAVYa} shows an estimated distribution of $\chi^2$ statistics using Algorithm~\ref{alg:SIS} which is similar to results in Figure~\ref{fig:DHR_NAV_2000}. This is because in neither DHR or SIS is called in these iterations.  For smaller  $\alpha$ (e.g., $\alpha=1,200$), range of $\chi^2$ statistics becomes very large.

    \begin{figure}[H]
  \centering
  \includegraphics[width=10cm,height=8cm]{DIST_SIS_2500nav.pdf}
   \caption{MCMC-SIS with $\alpha=2,500$.}
   \label{fig:MCMC_DIST_NAVYa}
    \end{figure}

\section{Conclusion}

In this paper we conduct MCMC with basic moves by allowing cell counts to be negative one and we also introduce two algorithms which merge basic MCMC methods, simulated annealing, and DHR and one that replaces DHR with SIS. 
Experimental results imply that the basic MCMC algorithm provides an estimated distribution of $\chi^2$ statistics closer to the asymptotic distribution, $\chi^2$ distribution with the degree of freedom $(I-1)(J-1)(K-1)$, than any hybrid model in most cases.  While MCMC-DHR and MCMC-SIS methods provide for broader sampling in the interior of the fiber there is difficulty in controlling the movement when either DHR or SIS is called.  This is especially problematic in high dimension tables with large cell counts.  However, hybrid models using Algorithms~\ref{alg:DHR} and~\ref{alg:SIS} perform suitably for sparse tables.

In addition, it is still an open problem to prove that $I \times J \times K$ tables for $I > 3, \, J > 3, \, K>3$ under the no-three-way interaction model moves of $2\times 2\times 2$ minors connect all tables in the fiber if we allow cell counts to be negative one.  Specifically 
\begin{conjecture}
Consider a three way $I \times J\times J$ contingency table
$X_{ijk}$ for $I > 3, \, J > 3, \, K>3$ under the no-three-way interaction model.
If we allow $X_{ijk} \geq -1$, then all moves of $2\times 2\times 2$ minors connect all tables in $\mathcal{F}_{A, b}$, where $A$ is a design matrix for the model and $b$ is a vector of margins computed from the observed table.
\end{conjecture}


