arXiv:1607.07600v1 [] 26 Jul 2016

An introduction to computational algebraic statistics
Satoshi Aoki
November 15, 2017
Abstract
In this paper, we introduce the fundamental notion of a Markov basis,
which is one of the first connections between commutative algebra and statistics. The notion of a Markov basis is first introduced by Diaconis and Sturmfels ([9]) for conditional testing problems on contingency tables by Markov
chain Monte Carlo methods. In this method, we make use of a connected
Markov chain over the given conditional sample space to estimate the p values
numerically for various conditional tests. A Markov basis plays an importance
role in this arguments, because it guarantees the connectivity of the chain,
which is needed for unbiasedness of the estimate, for arbitrary conditional
sample space. As another important point, a Markov basis is characterized as
generators of the well-specified toric ideals of polynomial rings. This connection between commutative algebra and statistics is the main result of [9]. After this first paper, a Markov basis is studied intensively by many researchers
both in commutative algebra and statistics, which yields an attractive field
called computational algebraic statistics. In this paper, we give a review of
the Markov chain Monte Carlo methods for contingency tables and Markov
bases, with some fundamental examples. We also give some computational
examples by algebraic software Macaulay2 ([11]) and statistical software R.
Readers can also find theoretical details of the problems considered in this
paper and various results on the structure and examples of Markov bases in
[3].

1

Conditional tests for contingency tables

A contingency table is a cross-classified table of frequencies. For example, suppose
40 students in some class took examinations of two subjects, Algebra and Statistics.
Suppose that both scores are classified to one of the categories, {Excellent, Good,
Fair}, and are summarized in Table 1. This is a typical example of two-way contingency tables. Since this table has 3 rows and 3 columns, this is called a 3 × 3
contingency table. The two subjects, Algebra and Statistics, are called factors of
the table, and the outcomes (i.e., scores) of each factor, {Excellent, Good, Fair},
are called levels of each factor. The cells of the I × J contingency table is the IJ
1

Table 1: Scores of Algebra and Statistics for 40 students (imaginary data)
Alg\Stat Excellent
Good
Fair
Total
Excellent
11
5
2
18
Good
4
9
1
14
Fair
2
3
3
8
Total
17
17
6
40

possible combinations of outcomes. Three-way, four-way or higher dimensional contingency tables are defined similarly. For example, adding to the data of Table 1,
if the scores of another subject (Geometry, for example) are also given, we Q
have a
three-way contingency table. An I1 × · · · × Im (m-way)contingency table has m
i=1 Ii
cells, where Ii is the number of levels for the ith factor, i = 1, . . . , m. In statistical
data analysis, the development of methods for analyzing contingency tables began
in the 1960s. We refer to [4] for standard textbook in this field.
We begin with simple I × J cases, and will consider generalizations to m-way
cases afterward. In statistical inference, we consider underlying random variables
and statistical models for observed data such as Table 1, and treat the observed
data as one realization of the random variables. In the case of Table 1, it is natural
to deal with the two-dimensional discrete random variables
(V1 , W1 ), (V2 , W2 ), . . . , (Vn , Wn ),

(1)

where n is the sample size, (n = 40 for Table 1) and (Vk , Wk ) is the couple of scores
obtained by the kth student. The random couples (Vk , Wk ) for k = 1, . . . , n are
drawn independently from the same distribution
P (Vk = i, Wk = j) = θij , i ∈ [I], j ∈ [J], k ∈ [n].
Here we use a notation [r] = {1, 2, . . . , r} for r ∈ Z≥0 , where Z≥0 is the set of
nonnegative integers. Note that we use appropriate coding such as 1: Excellent, 2:
Good, 3: Fair. The probability θ = (θij ) satisfies the condition
I X
J
X

θij = 1,

i=1 j=1

and is called a parameter. The parameter space
(
∆IJ−1 =

(θ11 , . . . , θIJ ) ∈ RIJ
≥0 :

I X
J
X
i=1 j=1

is called an IJ − 1 dimensional probability simplex.
2

θij = 1

)

To consider the data in the form of a contingency table, we also summarize the
underlying random variable (1) to the form of the contingency tables as
Xij =

n
X

1(Vk = i, Wk = j),

k=1

for i ∈ [I], j ∈ [J], where 1(·) is the indicator function. By this aggregation from
the raw scores to the contingency table, we neglect the order of observations in (1),
that is considered to have no information for estimating the parameter θ. Then the
data x = (xij ) ∈ ZIJ
≥0 is treated as a realization of X = (Xij ). The distribution of
X is a multinomial distribution given by
n!
p(x) = P (X = x) = I J
YY

xij !

I Y
J
Y
i=1 j=1

x
θijij ,

I X
J
X

xij = n.

(2)

i=1 j=1

i=1 j=1

We see that the multinomial distribution (2) is derived from the joint probability
function for n individuals under the assumption that each outcome is obtained
independently.
By summarizing the data in the form of contingency tables for fixed sample
size n, the degree of freedom of the observed frequency x becomes IJ − 1. which
coincides the degree of freedom of the parameter θ ∈ ∆IJ−1 . Here, we use “degree
of freedom” as the number of elements that are free to vary, that is a well-used
terminology in statistical fields. We can see the probability simplex ∆IJ−1 as an
example of statistical models, called a saturated model. Statistical model is called
saturated if the degree of freedom of the parameter equals to the degree of freedom
of data.
The saturated model is also characterized as the statistical model having the
parameter with the largest degree of freedom. In this sense, the saturated model
is the most complex statistical model. In other words, the saturated model is the
statistical model that fits the observed data perfectly, i.e., fits the data without error.
In fact, the parameter θ in the saturated model ∆IJ−1 is estimated from the data as
xij
θ̂ij =
, i ∈ [I], j ∈ [J],
(3)
n
that is also called an empirical probability of data. Because we assume that the
data x is obtained from some probability function such as multinomial distribution
(2) with some randomness, we want to consider more simple statistical model, i.e.,
a subset of the saturated model, M ⊂ ∆IJ−1 .
In the two-way contingency tables, a natural, representative statistical model is
an independence model.
Definition 1.1. The independence model for I × J contingency tables is the set
Mindp = {θ ∈ ∆IJ−1 : θij = θi+ θ+j , ∀i, ∀j},
3

(4)

where
θi+ =

J
X
j=1

θij , θ+j =

I
X

θij for i ∈ [I], j ∈ [J].

i=1

Remarks 1.2. Here we consider that only the sample size n is fixed. However,
several different situations can be considered for I × J contingency tables. The
situation that we consider here is called a multinomial sampling scheme. For other
sampling schemes such as Poisson, binomial and so on, see Chapter 2 of [4] or
Chapter 4 of [14]. Accordingly, the corresponding independence model Mindp is
called in different way for other sampling schemes. For example, it is called a
common proportions model for (product of) binomial sampling scheme where the
row sums are fixed, and main effect model for Poisson sampling scheme where no
marginal is fixed. Though there are also a little differences between the descriptions
of these models, we can treat these models almost in the same way by considering the
conditional probability function, which we consider afterward. Therefore we restrict
our arguments to the multinomial sampling scheme in this paper.
There are several equivalent descriptions for the independence model Mindp . The
most common parametric description in statistical textbooks is
Mindp = {θ ∈ ∆IJ−1 : θij = αi βj for some (αi ), (βj )} .

(5)

For other equivalent parametric descriptions or implicit descriptions, see Section 1
of [17], for example.
The meaning of Mindp in Table 1 is as follows. If Mindp is true, there are no
relations between the scores of two subjects. Then we can imagine that the scores
of two subjects follow the marginal probability functions for each score respectively,
and are independent, and the discrepancy we observed in Table 1 is obtained “by
chance”. However, it is natural to imagine some structure between the two scores
such as “there is a tendency that the students having better scores in Algebra are
likely to have better scores in Statistics”, because these subjects are in the same
mathematical category. In fact, we see relatively large frequencies 11 and 9 in the
diagonals of Table 1, which seem to indicate a positive correlation. Therefore one of
natural questions for Table 1 is “Is there some tendency between the two scores that
breaks independence?”. To answer this question, we evaluate the fitting of Mindp
by hypothetical testing.
The hypothetical testing problem that we consider in this paper is as follows.
H0 : θ ∈ Mindp v.s. H1 : θ ∈ ∆IJ−1 \ Mindp .
Here we call H0 a null hypothesis and H1 an alternative hypothesis. The terms null
model and alternative model are also used. The hypothetical testing in the above
form, i.e., a null model is a subset of a saturated model, M ⊂ ∆IJ−1 , and the
alternative model is the complementary set of M into the saturated model, is called
a goodness-of-fit test of model M. The testing procedures are composed of steps
such as choosing a test statistics, choosing a significance level, and calculating the
p value. We see these steps in order.
4

Choosing a test statistic. First we have to choose a test statistic to use. In
general, the term statistic means a function of the random variable X = (Xij ). For
example, (Xi+ ) and (X+j ) given by
Xi+ =

J
X

Xij , X+j =

j=1

I
X

Xij for i ∈ [I], j ∈ [J]

i=1

are examples of statistics called the row sums and the column sums, respectively.
Other examples of statistics are the row mean X̄i+ = Xi+ /J and the column mean
X̄+j /I for i ∈ [I], j ∈ [J]. To perform the hypothetical testing, we first select an appropriate statistic, called a test statistic, to measure the discrepancy of the observed
data from the null model. One of the common test statistic for the goodness-of-fit
test is a Pearson goodness-of-fit χ2 given by
χ2 (X) =

I X
J
X
(Xij − m̂ij )2
i=1 j=1

m̂ij

,

where m̂ij is the fitted value of Xij under H0 , i.e., an estimator of E(Xij ) = mij =
nθij , given by
xi+ x+j
.
(6)
m̂ij = nθ̂ij =
n
Here we use the maximum likelihood estimate of the parameter under the null model,
θ̂ = (θ̂ij ), given by
xi+ x+j
θ̂ij =
,
(7)
n2
that is obtained by maximizing the log-likelihood
Const +

I X
J
X

xij log θij

i=1 j=1

under the constraint θ ∈ Mindp . The meaning of this estimate is also clear in a
parametric description (5) since the maximum likelihood estimates of (αi ), (βj ) are
given by
xi+
x+j
α̂i =
, β̂j =
,
n
n
respectively. The fitted value for Table 1 under Mindp is given in Table 2.
There are various test statistics other than the Pearson goodness-of-fit χ2 that
can be used in our problem. Another representative is the (twice log) likelihood
ratio given by
I X
J
X
Xij
,
(8)
2
Xij log
m̂
ij
i=1 j=1
where m̂ij is given by (6). In general, test statistic should be selected by considering
their power, i.e., the probability that the null hypothesis is rejected if the alternative
5

Table 2: The fitted value under Mindp for Table 1
Alg\Stat Excellent
Good
Fair
Total
Excellent
7.65
7.65
2.70
18
Good
5.95
5.95
2.10
14
Fair
3.40
3.40
1.20
8
Total
17
17
6
40

hypothesis is true. See textbooks such as [15] for the theory of the hypothetical
testing, the optimality of the test statistics, examples and the guidelines for choosing
test statistics for various problems.
Choosing a significance level. Once we choose a test statistic to use, as the
Pearson goodness-of-fit χ2 for example, the hypothetical testing procedure is written
by
χ2 (xo ) ≥ cα ⇒ Reject H0 ,
where xo is the observed data, and cα is the critical point at the significance level α
satisfying
P (χ2 (X) ≥ cα | H0 ) ≤ α.
(9)
The probability of the left hand side of (9) is called a type I error. Equivalently, we
define the p-value by
p = P (χ2(X) ≥ χ2 (xo ) | H0 ),
(10)
then the testing procedure is written by
p ≤ α ⇒ Reject H0 .
The meaning of the p-value for the data xo is the conditional probability that “more
or equally discrepant results are obtained than the observed data if the null hypothesis is true”. Therefore, if p-value is significantly small, we conclude that null
hypothesis is unrealistic, because it is doubtful that such an extreme result xo is
obtained. This is the idea of the statistical hypothetical testing. In this process, the
significance level α plays a threshold to decide the p-value is “significantly small”
to reject the null hypothesis. In statistical and scientific literature, it is common to
choose α = 0.05 or α = 0.01. Readers can find various topics on p-value in [19].
Calculating the p-value. Once we choose a test statistic and a significance level,
all we have to do is to calculate the p-value given in (10) for observed data xo . The
observed value of the Pearson goodness-of-fit χ2 for Table 1 is
3 X
3
X
(xoij − m̂ij )2
(11 − 7.65)2
(3 − 1.20)2
=
+···+
= 8.6687,
χ (x ) =
m̂ij
7.65
1.20
i=1 j=1
2

o

6

therefore the p-value for our xo is
p = P (χ2 (X) ≥ 8.6687 | H0 ).
This probability is evaluated based on the probability function of the test statistic
χ2 (X) under H0 , which we call a null distribution hereafter. Unfortunately, the null
distribution depends on the unknown parameter θ ∈ Mindp and the p-values cannot
be calculated in most cases in principle. One naive idea to evaluate the p-values for
such cases is to calculate its supremum in Mindp and perform the test as the form
sup P (χ2 (X) ≥ χ2 (xo ) | H0 ) ≤ α ⇒ Reject H0 .

(11)

θ∈Mindp

However, this idea is hard to implement in general, i.e., it is usually difficult to
evaluate the left-hand side of (11) or to seek tests that are powerful under (11).
Then, what should we do? We consider the following three strategies for calculating
p-values in this paper.
(a) Using the asymptotic distribution of the test statistic.
(b) Exact calculation based on the conditional distribution.
(c) Estimate the p-value by the Monte Carlo method.
The aim of this paper is to introduce strategy (c). We will consider each strategy
in order.
(a) Using the asymptotic distribution of the test statistic. In applications,
it is common to rely on various asymptotic theories for the test statistics. As for
the Pearson goodness-of-fit test χ2 test, the following result is known.
Theorem 1.3. Under the null model Mindp , the Pearson goodness-of-fit χ2 (X)
asymptotically follows the χ2 distribution with (I − 1)(J − 1) degree of freedom,
i.e.,
lim P (χ2 (X) ≥ u) = P (V ≥ u) for u > 0,
n→∞

χ2(I−1)(J−1) ,

where V ∼
degree of freedom.

i.e., V is distributed to the χ2 distribution with (I − 1)(J − 1)

This theorem is shown as a consequence of the central limit theorem. In addition,
the same asymptotic distribution is given when we consider the conditional limit,
i.e., consider n → ∞ under the condition that Xi+ /n → ai and X+j /n → bj for
i ∈ [I], j ∈ [J] for some fixed 0 < ai , bj < 1. See [6] or [18] for detail. Anyway, these
asymptotic properties are the reason why we call this test as Pearson goodness-of-fit
“χ2 test”. Similarly to the Pearson goodness-of-fit χ2 , there are several test statistics
that have the χ2 distribution as the asymptotic distribution. An important example
is the likelihood ratio test statistic, which is given in (8) for our setting. Moreover,
7

0.00

0.05

0.10

0.15

several asymptotic good properties of likelihood ratio test statistics are known. See
[15] for details. Note also that our methods, Markov chain Monte Carlo methods,
can be applicable for arbitrary type of test statistics, though we only consider the
Pearson goodness-of-fit χ2 in this paper.
Following Theorem 1.3, it is easy to evaluate the asymptotic p-value of the
Pearson goodness-of-fit χ2 test. For our data, the observed value of test statistic,
χ2 (xo ) = 8.6687, is less than the upper 5 percent point of the χ2 distribution with
4 degrees of freedom, χ24,0.05 = 9.488. Therefore, for the significance level α = 0.05,
we cannot reject the null hypothesis H0 , i.e., we cannot say that “the fitting of the
model Mindp to Table 1 is poor”. Equivalently, the asymptotic p-value is calculated
as the upper probability of χ24 , which is 0.0699 and is greater than α = 0.05. Figure
1 presents the probability density function of the χ24 distribution. The above results

0

5

10

15

20

25

30

Figure 1: χ2 distribution with degree of freedom 4. The vertical solid line indicates
the observed value χ2 (xo ) = 8.6687, and the dotted line indicates the critical point
for the significance level α = 0.05, χ24,0.05 = 9.488.
can be obtained numerically by the following codes of the statistical software R.
> x <- matrix(c(11,5,2,4,9,1,2,3,3), byrow=T, ncol=3, nrow=3)
> x
[,1] [,2] [,3]
[1,]
11
5
2
[2,]
4
9
1
[3,]
2
3
3
> chisq.test(x)
Pearson’s Chi-squared test
8

data: x
X-squared = 8.6687, df = 4, p-value = 0.06994
> pchisq(8.6687,4, lower.tail=F)
[1] 0.06993543
> qchisq(0.05,4,lower.tail=F)
[1] 9.487729

# critical point

As we see above, using asymptotic null distribution is easy way to evaluate
p-values, and one of the most common approaches in applications. One of the disadvantages of strategy (a) is that there might not be a good fit with the asymptotic
distribution. In fact, because sample size is only n = 40 for Table 1, it is doubtful that we can apply the asymptotic result of n → ∞. Besides, it is well known
that there is cases that the fitting of the asymptotic distributions are poor for data
with relatively large sample sizes. One such case is sparse data case, another one is
unbalanced case. See [12] for these topics.
(b) Exact calculation based on the conditional distribution. If we want
to avoid asymptotic approaches as strategy (a), an alternative choice is to calculate
p-values exactly. For the cases that the null distribution of the test statistics depend
on the unknown parameters, we can formulate the exact methods based on the
conditional probability functions for fixed minimal sufficient statistics under the null
model Mindp . The key notion here is the minimal sufficient statistics.
Definition 1.4. Let X be a discrete random variable with the probability function
p(x) with the parameter θ. The statistic T(X), i.e., a vector or a scalar function of
X, is called sufficient for θ if the conditional probability function of X for a given
T,
p(x | t) = P (X = x | T(X) = t),
(12)
does not depend on θ. The sufficient statistic T(X) is minimal if there is no other
sufficient statistics that is a function of T(X).
The meaning of the minimal sufficient statistic is explained as follows. If we
know the value of T, then knowing X provides no further information about the
parameter θ. Therefore for the parameter estimation or hypothetical testing, it is
sufficient to consider the methods based on the minimal sufficient statistic. The
minimal sufficient statistics for our two-way problem is as follows.
• Under the saturated model θ ∈ ∆IJ−1 , a minimal sufficient statistic is the
contingency table X. Adding the additional information such as the scores
of the kth student, (Vk , Wk ), in (1) gives us no additional information on the
estimation of θ. Indeed, under the saturated model, the maximum likelihood
estimate of the parameter is the empirical probability (3), that is a function
of the minimal sufficient statistic.
9

• Under the independence model Mindp , a minimal sufficient statistic is the row
sums {Xi+ , i ∈ [I]} and the column sums {X+j , j ∈ [J]}, as we see below.
Indeed, we have already seen that the maximum likelihood estimate of the
parameter under the independence model is (7), that is a function of the row
sums and column sums. Note that X itself is also the sufficient statistic under
the independence model, but is not minimal.
To see that a given statistic T(X) is sufficient for a parameter θ, a useful way is
to rely on the following theorem.
Theorem 1.5. T(X) is a sufficient statistic for θ if and only if the probability
function of X is factored as
p(x; θ) = h(x)g(T (x); θ),

(13)

where g(·) is a function that depends on the parameter θ and h(·) is a function that
does not.
For the case of discrete probability function, this theorem, called a factorization
theorem, is easily (i.e., without measure theories) proved from the definition of the
sufficient statistic. Generally, to obtain such a factorization is easier than to compute explicitly the conditional distribution (12). For example, under the parametric
description θij = αi βj , the probability function of the multinomial distribution (2)
is written as
!
!
Y x
Y x
n!
αi i+
p(x; θ) = Q Q
βj +j
xij !
i
j

and we see that T (X) = ({Xi+ }, {X+j }) is a sufficient statistic for the parameter
θ ∈ Mindp .
Here, for later generalization, we introduce a configuration matrix A and express
a minimal sufficient statistic by A as follows. Let the number of the cells of the
contingency table X be p and treat X as a p-dimensional column vector. Let T (X)
be a d-dimensional sufficient statistic for the parameter θ ∈ M. For example of
the independence model Mindp for I × J contingency tables, we have ν = IJ,
X = (X11 , X12 , . . . , XIJ )′ and
T(X) = (X1+ , . . . , XI+ , X+1 , . . . , X+J )′
and d = I + J. Then we see that T (X) is written as
T (X) = AX

(14)

for d × ν integer matrix A. For the 3 × 3 contingency tables, A is written as follows:


1 1 1 0 0 0 0 0 0
 0 0 0 1 1 1 0 0 0 


 0 0 0 0 0 0 1 1 1 

(15)
A=
 1 0 0 1 0 0 1 0 0 .


 0 1 0 0 1 0 0 1 0 
0 0 1 0 0 1 0 0 1
10

Following the sufficiency of T (X) = AX, the conditional probability function for
given T = t does not depend on the parameter. For the case of the independence
model Mindp for two-way contingency tables, it is


h(x) = P (X = x | AX = t, H0 ) =




Y
i




xi+ !

n!

Y

Y
j



x+j !


xij !

,

x ∈ Ft ,

(16)

i,j

where
Ft = {x ∈ Zν≥0 : Ax = t}
is the conditional sample space, which is called a t-fiber in the arguments of Markov
bases. The conditional probability function h(x) is called a hypergeometric distribution. Using this conditional probability, the conditional p-value can be defined
by
X
g(x)h(x)
(17)
p = EH0 (g(X) | AX = Axo ) =
x∈FAxo

for the observed table xo , where g(x) is the test function

1, χ2 (x) ≥ χ2 (xo ),
g(x) =
0, otherwise.

Now calculate the conditional p-value exactly for Table 1. For the observed table xo ,
i.e., Table 1, we consider the independence model Mindp . The configuration matrix
A for Mindp is given in (15). The t-fiber including xo , i.e., Axo -fiber, is the set of all
contingency tables that have the same value of the row sums and the column sums
to xo ,


x
x
x
18


11
12
13




x
x
x
14
21
22
23
9
FAxo = x ∈ Z≥0 :
.
x31 x32 x33 8 





17 17 6 40

There are 2366 elements in this FAxo . For each 2366 elements in FAxo , the conditional probability is given by
(18!14!8!) (17!17!6!) Y 1
h(x) =
, x ∈ FAxo .
40!
xij !
i,j
Then we have the exact conditional p-value
X
g(x)h(x) = 0.07035480,
p=
x∈FAxo

where the test function is
g(x) =



1, χ2 (x) ≥ 8.6687,
0, otherwise.
11

As a result, we cannot reject H0 at significance level 0.05, which is the same result
to strategy (a).
Example 1.6. The following toy example should help the reader in understanding
the method. Let consider the 2 × 3 contingency table with the row sums and the
column sums given as follows.
x11 x12 x13 3
x21 x22 x23 2
2
2
1 5
There are 5 elements in the fiber as

2 1 0
2 0 1
,
,
F(3,2,2,2,1) =
0 1 1
0 2 0
= {x1 , x2 , x3 , x4 , x5 }.

1 2 0
,
1 0 1

1 1 1
,
1 1 0

0 2 1
2 0 0



1.2 1.2 0.6
. Then the Pearson goodness-of0.8 0.8 0.4
fit χ2 for each element is calculated as
The fitted value under the Mindp is

(χ2 (x1 ), χ2 (x2 ), χ2 (x3 ), χ2 (x4 ), χ2 (x5 )) = (2.917, 5, 2.917, 0.833, 5).
The conditional probabilities
h(x) =

2Y 1
3!2!2!2! Y 1
=
5!
xij !
5 i,j xij !
i,j

for each element are calculated as
(h(x1 ), h(x2 ), h(x3 ), h(x4 ), h(x5 )) = (0.2, 0.1, 0.2, 0.4, 0.1).
Therefore the conditional p-value for x4 is 1.0, that for x1 or x3 is 0.6, and that for
x2 or x5 is 0.2.
Remark 1.7. We briefly mention the generalization of the above method to general
problems and models. First important point is the existence of the minimal sufficient
statistics in the form of (14). It is known that, for the exponential family, well-known
family of the distribution, minimal sufficient statistics exist, and for a special case
of the exponential family, called the toric model, minimal sufficient statistics of the
form (14) exist. The toric model is relatively new concept arising in the field of
the computational algebraic statistics and is defined from the configuration matrix
A = (aij ) ∈ Zd×ν
≥0 as follows. For the jth column vector aj = (a1j , . . . , adj ) of A,
j ∈ [ν], define the monomial
θ

aj

=

d
Y

a

θi ij , j ∈ [ν].

i=1

12

Then the toric model of A is the image of the orthant Rd>0 under the map
f : Rd → Rν , θ 7→ Pν

1

aj
j=1 θ

(θa1 , . . . , θaν ).

See Chapter 1.2 of [17] for detail. The toric model specified by the configuration
matrix A ∈ Zd×ν
≥0 is also written by
MA = {θ = (θi ) ∈ ∆ν−1 : log θ ∈ rowspan(A)},
where rowspan(A) = image(A′ ) is the linear space spanned by the rows of A, and
log θ = (log θ1 , . . . , log θν )′ , where ′ is a transpose. In statistical fields, this is
called a log-linear model. In fact, for example of the independence model Mindp of
2 × 3 tables, that is a log-linear model, the parametric description θij = αi βj can be
written as

 


1 0 1 0 0 
log θ11
α
1
 log θ12   1 0 0 1 0 


 

 log θ13   1 0 0 0 1   α2 


 

 log θ21  =  0 1 1 0 0   β1  .

 


 log θ22   0 1 0 1 0  β2
β3
0 1 0 0 1
log θ23

The conditional probability function, i.e., the generalization of the hypergeometric
distribution h(x) in (16) is as follows. For the model specified by the configuration
matrix A, the conditional probability function for given sufficient statistic Axo is
1
−1
P (X = x | AX = Axo ) = CAx
o Y

xi !

,

i∈[ν]

where
CAxo =

X

y∈FAxo

1
Y

yi !

(18)

i∈[ν]

is a normalizing constant. Based on this conditional probability function, we can
calculate the conditional p-values by (17).
Finally, we note an optimality of the method briefly. The conditional procedure
mentioned above is justified if we consider the hypothetical testing to the class of
similar tests and the minimal sufficient statistics is complete. For the class of the
exponential family, it is known that the minimal sufficient statistic is complete. See
Chapter 4.3 of [15] for detail.
(c) Estimate the p-value by the Monte Carlo method. The two strategies to
evaluate p-values we have considered, asymptotic evaluation and exact computation,
have both advantages and disadvantages. The asymptotic evaluations relying on the
13

asymptotic χ2 distribution are easy to carry out, especially by various packages in
softwares such as R. However, poor fitting to the asymptotic distribution can not be
ignorable for sparse or unbalanced data even with relatively large sample sizes. The
exact calculation of the conditional p-values is the best method if it is possible to
carry out. In fact, various exact methods and algorithms are considered for problems
of various types of the contingency tables, statistical models and test statistics. See
the survey paper [2] for this field. However, for large size samples, the cardinality of
the fiber |FAxo | can exceed billions, making exact computations difficult to be carried
out. In fact, it is known that the cardinality of a fiber increases exponentially in
the sample size n. (An approximation for the cardinality of a fiber is given by [10].)
For these cases, the Monte Carlo methods can be effective.
The Monte Carlo methods estimate the p-values as follows. To compute the
conditional p-value (17), generate samples
P x1 , . . . , xN from the null distribution h(x).
Then the p-value is estimated as p̂ = N
i=1 g(xi )/N, that is an unbiased estimate
of the p-value. We can set N according to the performance of our computer. As
an advantage of the Monte Carlo method, we can also estimate the accuracy, i.e.,
variance ofpthe estimate. For example, a conventional 95% confidence interval of
p, p̂ ± 1.96 p̂(1 − p̂)/N, is frequently used. The problem here is how to generate
samples from the null distribution. We consider Markov chain Monte Carlo methods,
often abbreviated as the MCMC methods, in this paper.
Following MCMC methods setup, we construct an ergodic Markov chain on the
fiber F = FAxo whose stationary distribution is prescribed, given by (16). Let the
elements of F be numbered as
F = {x1 , . . . , xs }.
We write the null distribution on F as
π = (π1 , . . . , πs ) = (h(x1 ), . . . , h(xs )).
Here, by standard notation, we treat π as a row vector. We write the transition
probability matrix of the Markov chain {Zt , t ∈ Z≥0 } over F as Q = (qij ), i.e., we
define
qij = P (Zt+1 = xj | Zt = xi ).
Then a probability distribution θ ∈ ∆s−1 is called a stationary distribution if it
satisfies θ = θQ. The stationary distribution uniquely exists if the Markov chain
is irreducible, (i.e., connected in this case) and aperiodic. Therefore for the connected and aperiodic Markov chain, starting from an arbitrary state Z0 = xi , the
distribution of Zt for large t is close to its stationary distribution. If we can construct a connected and aperiodic Markov chain with the stationary distribution π,
by running the Markov chain and discarding a large number t of initial steps (called
burn-in steps), we can treat Zt+1 , Zt+2 , . . . to be samples from the null distribution
π and use them to estimate p-values. Then the problem becomes how to construct
a connected and aperiodic Markov chain with the stationary distribution as the null
14

distribution π over F . Among these conditions, the conditions for the stationary
distribution can be solved easily. Once we construct an arbitrary connected chain
over F , we can modify its stationary distribution to the given null distribution π as
follows.
Theorem 1.8 (Metropolis-Hastings algorithm). Let π be a probability distribution
on F . Let R = (rij ) be the transition probability matrix of a connected, aperiodic and
symmetric Markov chain over F . Then the transition probability matrix Q = (qij )
defined by


πj
qij = rij min 1,
, i 6= j
πi
X
qij
qii = 1 −
j6=i

satisfies π = πQ.

This theorem is a special case of [13]. Though the symmetry assumption (rij =
rji ) can be removed easily, we only consider symmetric R for simplicity. The proof
of this theorem is easy and is omitted. See [13] or Chapter 4.1 of [14], for example.
Instead, we consider the algorithm for data of small size.
Example 1.9. Consider the small example in Example 1.6. As we have seen, the
fiber is
F = {x1 , x2 , x3 , x4 , x5 }
and the null distribution is
π = (π1 , . . . , π5 ) = (h(x1 ), . . . , h(x5 )) = (0.2, 0.1, 0.2, 0.4, 0.1).
Using the Markov basis we consider in the next section, we can construct a connected, aperiodic and symmetric Markov chain with the transition probability matrix


1/2 1/6 1/6 1/6 0
 1/6 2/3 0 1/6 0 



(19)
R=
 1/6 0 1/2 1/6 1/6  .
 1/6 1/6 1/6 1/3 1/6 
0
0 1/6 1/6 2/3
Following Theorem 1.8, we modify the Markov
bility matrix

7/12 1/12 1/6
 1/6 2/3
0

1/6
0
7/12
Q=

 1/12 1/24 1/12
0
0
1/6

chain to have the transition proba-


1/6
0
1/6
0 

1/6 1/12 
 .
3/4 1/24 
1/6 2/3

We can check that the eigenvector from the left of Q with the eigenvalue 1 is π. We
can also check that each row vector of QT for large T converges to π.
15

An important advantage of the Markov chain Monte Carlo method is that it does
not require the explicit evaluation of the normalizing constant of the null distribution. As is shown in Theorem 1.8, we only need to know π up to a multiplicative
constant, because the normalizing constant, (18) in the general form, canceled in the
ratio πj /πi . With Theorem 1.8, the remaining problem is to construct an arbitrary
connected and aperiodic Markov chain over F , that is solved by the Gröbner basis
theory.

2

Markov bases and ideals

As stated in the previous section, the main task for estimating p-values thanks
to MCMC methods is to construct a connected and aperiodic Markov chain over
F = FAxo with stationary distribution given by (16). Here, A ∈ Zd×ν is a given
configuration matrix, xo ∈ Zν≥0 is the observed contingency table and FAxo , a Axo fiber, is the set of all contingency tables with the same value of the minimal sufficient
statistics to xo ,
FAxo = {x ∈ Zν≥0 : Ax = Axo }.
We write the integer kernel of A as
KerZ (A) = Ker(A) ∩ Zν = {z ∈ Zν : Az = 0}.
An element of KerZ (A) is called a move. Note that x − y ∈ KerZ (A) if and only
if x, y ∈ FAx . Then for a given subset B ⊂ KerZ (A) and t ∈ Zd≥0 , we can define
undirected graph Gt,B = (V, E) by
V = Ft ,

E = {(x, y) : x − y ∈ B or y − x ∈ B}.

Definition 2.1 (A Markov basis). B ⊂ KerZ (A) is a Markov basis for A if Gt,B is
connected for arbitrary t ∈ Zd≥0 .
Once we obtain a Markov basis B for A, we can construct a connected Markov
chain over FAxo easily as follows. For each state x ∈ FAxo , randomly choose a move
z ∈ B and a sign ε ∈ {−1, 1} and consider x + εz. If x + εz ∈ FAxo , then x + εz is
the next state, otherwise stay at x. Then we have the connected Markov chain over
Axo . We see these arguments in an example.
Example 2.2. Again we consider a small data of Example 1.6, where the fiber is
redisplayed below.


2 1 0
2 0 1
1 2 0
1 1 1
0 2 1
F(3,2,2,2,1) =
,
,
,
,
0 1 1
0 2 0
1 0 1
1 1 0
2 0 0
= {x1 , x2 , x3 , x4 , x5 }.

16

The integer kernel for the configuration matrix

1 1 1 0
 0 0 0 1

A=
 1 0 0 1
 0 1 0 0
0 0 1 0

0
1
0
1
0

includes moves such as
z1 =

0
1
0
0
1








(20)

1 −1 0
1 0 −1
0
1 −1
, z2 =
, z3 =
,
−1
1 0
−1 0
1
0 −1
1

2 −1 −1
, ....
−2
1
1

From these, we consider some sets of moves. If we consider B1 = {z1 }, corresponding undirected graph G(3,2,2,2,1),B1 is given in Figure 2(a), which is not connected.
Therefore B1 is not a Markov basis. If we consider B2 = {z1 , z2 }, corresponding
undirected graph G(3,2,2,2,1),B2 is given in Figure 2(b), which is connected. However,
B2 is also not a Markov basis, because there exists t ∈ Z5≥0 where Gt,B2 is not connected. An example of such t is t = (1, 1, 0, 1, 1), with the corresponding t-fiber is
a two-element set


0 1 0
0 0 1
,
.
(21)
F(1,1,0,1,1) =
0 0 1
0 1 0
The above example shows that a Markov basis includes z3 to connect the two elements above. In fact, B = {z1 , z2 , z3 } is a Markov basis for this A, with the corresponding undirected graph G(3,2,2,2,1),B3 in Figure 2(c). The transition probability
x♥
1
x♥
2

❅
❅

x♥
4

x♥
1

x♥
3
x♥
5

(a)

✁✁ ❅❅

x♥
2 ✁
✁
✁
x♥
4

x♥
1

x♥
3
x♥
5

(b)

x♥
2 ✁

✁✁ ❅❅

x♥
3
✁ ✑✑
✑
✁✑
x♥
x♥
4
5
(c)

Figure 2: Undirected graphs for B1 , B2 , B3 for t = (3, 2, 2, 2, 1).
,
matrix (19) in Example 1.9 corresponds to a Markov chain constructed from B3 as
“in each step, choose 3 elements in B3 and its sign {−1, 1} with equal probabilities”.
At first sight, we may feel the cases such as (21) are trivial and may imagine that
“if we only consider the cases with t ∈ Zd>0 , i.e., cases with strictly positive minimal
sufficient statistics (that may be realistic situations in the actual data analysis), it is
easy to connect the fiber Ft ”. However, it is not so. We will see an example where
complicated moves are needed even for the fiber with positive t.
17

The connection between the Markov basis and a toric ideal of a polynomial ring
by [9] is as follows. Let k[u] = k[u1 , u2 , . . . , uν ] denote the ring of polynomials in ν
variables over a field k. Let a contingency table x ∈ Zν≥0 be mapped to the monomial
ux ∈ k[u], and a move, i.e., an element of the integer kernel z = z+ − z− ∈ KerZ (A),
+
−
be mapped to the binomial uz − uz ∈ k[u]. For the case of the independence
model for the 3 × 3 contingency tables, examples of these correspondences are as
follows.
11
5
2
5 2 4 9
2 3 3
4
9
1 ⇐⇒ u11
11 u12 u13 u21 u22 u23 u31 u32 u33
2
3
3
2 −1 −1
−3
1
2 ⇐⇒ u211 u22 u223 u31 − u12 u13 u321 u33
1
0 −1
The binomial ideal in k[u] generated by the set of binomials corresponding to the
set of moves for A,
Dn +
oE
z
z−
+
−
IA =
u − u : z − z ∈ KerZ (A) ,

is the the toric ideal of configuration A.

Theorem 2.3 (Theorem 3.1 of [9]). B = {z1 , . . . , zL } ⊂ KerZ (A) is a Markov basis
+
−
for A if and only if {uzi − uzi , i = 1, . . . , L} generates IA .
A proof of Theorem 2.3 is given in the original paper [9]. We can also find more
detailed proof in Chapter 4 of [14]. In these proofs, the sufficiency and the necessity
are shown by induction on some integer. In the proof of sufficiency, this integer
represents the number of steps of the chain, and the argument is straightforward.
On the other hand, in the proof of necessity, this integer represents the number of
terms in the expansion that we want to show in the proof, and is not necessarily
equal to the number of steps of the chain. Theorem 2.3 shows a non-trivial result
on this point.
To calculate a Markov basis for a given configuration matrix A, we can use the
elimination theory. For this purpose, we also prepare variables v = {v1 , . . . , vd } for
the minimal sufficient statistic t and consider the polynomial ring k[v] = k[v1 , . . . , vd ].
The relation t = Ax can be expressed by the homomorphism
ψA : k[u] → k[v]
a
a
a
uj 7→ v1 1j v2 2j · · · vd dj .
Then the toric ideal IA is also expressed as IA = Ker(ψA ). We now have the
following.
Corollary 2.4 (Theorem 3.2 of [9]). Let IA∗ be the ideal of k[u, v] given by
IA∗ = h−ψA (uj ) + uj , j = 1, . . . , νi ⊂ k[u, v].
Then we have IA = IA∗ ∩ k[u].
18

Corollary 2.4 suggests that we can obtain a generator of IA as its Gröbner basis
for an appropriate term order called an elimination order. For an ideal J ∈ k[u] and
a term order ≺, a set of polynomials {g1 , . . . , gs }, g1 , . . . , gs ∈ J, is called a Gröbner
basis of I with respect to a term order ≺, if {in≺ (g1 ), . . . , in≺ (gs )} generates an initial
ideal of J defined by in≺ (J) = h{in≺ (f ) : 0 6= f ∈ J}i. Here we write in≺ (f ) as an
initial term of f with respect to a term order ≺. For more theories and results on
Gröbner bases, see textbooks such as [7]. The elimination theory is one of the useful
applications of Gröbner bases and is used for our problem as follows. For the reduced
Gröbner basis G∗ of IA∗ for any term order satisfying {v1 , . . . , vd } ≻ {u1, . . . , uν },
G∗ ∩ k[u] is a reduced Gröbner basis of IA . Because the Gröbner basis is a generator
of IA , we can obtain a Markov basis for A as the reduced Gröbner basis in this way.
The computations of Gröbner bases can be carried out by various algebraic
softwares such as Macaulay2 ([11]), SINGULAR ([8]), CoCoA ([5]), Risa/Asir ([16])
and 4ti2 ([1]). Here, we show some computations by Macaulay2, because we can
also rapidly use it online at the website 1 . We start with a simple example.
Example 2.5. In Example 2.2, we give a Markov basis for the independence model
for 2 × 3 contingency tables without any proof or calculations. Here we check that
the set of 3 moves


1 −1 0
1 0 −1
0
1 −1
, z2 =
, z3 =
z1 =
−1
1 0
−1 0
1
0 −1
1
constitute a Markov basis for A given in (20). In other words, we check that the
corresponding toric ideal IA is generated by 3 binomials
{u11 u22 − u12 u21 , u11 u23 − u13 u21 , u12 u23 − u13 u22 }.

(22)

Following Corollary 2.4, we prepare the variable v = (v1 , . . . , v5 ) for the row sums
and column sums of x as
x11 x12 x13 v1
x21 x22 x23 v2
v3 v4 v5
and consider the homomorphism
u11 7→ v1 v3 , u12 7→ v1 v4 , u13 7→ v1 v5 ,
u21 7→ v2 v3 , u22 7→ v2 v4 , u23 7→ v2 v5 .
Then under the elimination order v ≻ u, compute the reduced Gröbner basis of the
toric ideal
IA∗ = h−v1 v3 + u11 , −v1 v4 + u12 , . . . , −v2 v5 + u23 i .
These calculations are done by Macaulay2 as follows.
1

Macaulay2 online: http://habanero.math.cornell.edu.3690

19

i1 : R=QQ[v1,v2,v3,v4,v5,u11,u12,u13,u21,u22,u23,MonomialOrder=>{5,6}]
o1 = R
o1 : PolynomialRing
i2 : I=ideal(-v1*v3+u11,-v1*v4+u12,-v1*v5+u13,-v2*v3+u21,-v2*v4+u22,-v2*v5+u23)
o2 = ideal (- v1*v3 + u11, - v1*v4 + u12, - v1*v5 + u13, - v2*v3 + u21, - v2*v4
-------------------------------------------------------------------------+ u22, - v2*v5 + u23)
o2 : Ideal of R
i3 : G=gb(I); g=gens(G)
o4 = | u13u22-u12u23 u13u21-u11u23 u12u21-u11u22 v4u23-v5u22 v4u13-v5u12
-------------------------------------------------------------------------v3u23-v5u21 v3u22-v4u21 v3u13-v5u11 v3u12-v4u11 v1u23-v2u13 v1u22-v2u12
-------------------------------------------------------------------------v1u21-v2u11 v2v5-u23 v1v5-u13 v2v4-u22 v1v4-u12 v2v3-u21 v1v3-u11 |
1
o4 : Matrix R

18
<--- R

i5 : selectInSubring(1,g)
o5 = | u13u22-u12u23 u13u21-u11u23 u12u21-u11u22 |
1
o5 : Matrix R

3
<--- R

The output o4 shows the reduced Gröbner basis of IA∗ under the elimination
(reverse lexicographic) order v ≻ u, and the output o5 shows the reduced Gröbner
basis of IA , which we can use as a Markov basis. We have now checked a Markov
basis (22).
From the Markov basis (22), we may imagine that the set of moves corresponding
to the binomials
{uij ui′ j ′ − uij ′ ui′ j , 1 ≤ i < i′ ≤ I, 1 ≤ j < j ′ ≤ J}
forms a Markov basis for the independence model of the I × J contingency tables,
which is actually true. This fact is given and proved as Theorem 2.1 of [3], for
example.
Now we are ready to estimate p-value for our original problem of 3×3 contingency
table in Table 1. The Markov basis for this problem is formed by 9 moves of the
above type. Using this Markov basis, we calculate the conditional p-values for Table
1 by the Markov chain Monte Carlo method. For each step of the chain, we choose
an element of the Markov basis randomly, and modify the transition probability by
20

0.00

0.05

Density

0.10

0.15

Theorem 1.8. We start the chain at the observed table xo of Table 1, discard initial
50000 steps as the burn-in steps, and have 100000 samples of the Pearson goodnessof-fit χ2 . Figure 3 is a histogram of the sampled Pearson goodness-of-fit χ2 with
the asymptotic χ24 distribution. In these 100000 samples, 6681 samples are larger

0

5

10

15

20

25

30

Pearson chi−squares

Figure 3: A histogram of sampled Pearson χ2 goodness-of-fit for Table 1 generated
by a Markov chain Monte Carlo method. The dotted curve is the corresponding
asymptotic χ24 distribution.
than or equal to the observed value χ2 (xo ) = 8.6687, then we have the estimate
p̂ = 0.06681. Therefore we cannot reject H0 at significance level 0.05, which is the
same result to the other strategies (a) and (b). Though the difference from the exact
value p = 0.07035480 from the simulated value is slightly larger than the asymptotic
estimate (p̂ = 0.0699), we may increase the accuracy of the estimates by increasing
the sample sizes. To compare the three strategies for Table 1, we compute the
upper percentiles of 90%, 95%, 99%, 99.9% for (a) asymptotic χ24 distribution, (b)
exact conditional distribution and (c) Monte Carlo simulated distribution in Table
3.
Finally, we give an example for which the structure of the Markov basis is complicated. The model we consider is a no three-factor interaction model for three-way
contingency tables. The parametric description of the no three-factor interaction
model is given by
Mn3 = {θ ∈ ∆ : θijk = αij βik γjk for some (αij ), (βik ), (γjk )}.
This is one of the most important statistical models in the statistical data analysis
of three-way contingency tables. The minimal sufficient statistics for Mn3 is the
21

Table 3: The upper percentiles for three strategies of Pearson
Table 1.
90% 95%
(a) Asymptotic χ24 distribution
7.779 9.488
(b) Exact null distribution
7.766 9.353
(c) Monte Carlo simulated distribution 7.684 9.287

goodness-of-fit χ2 for
99%
13.28
12.78
12.73

99.9%
18.47
17.99
18.58

two-dimensional marginals
{xij+ }, {xi+k }, {x+jk },
where we define
xij+ =

K
X

xijk , xi+k =

J
X

xijk , x+jk =

j=1

k=1

I
X

xijk .

i=1

We only consider 3 × 3 × 3 case (i.e., I = J = K = 3) here. Then the configuration
matrix A is 27 × 27 matrix written as follows.
 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
0








A=







0

0

1

1

1

0

0

0

0

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

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
1
0

0
0
1

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
1
0

0
0
1

0
0
0

0
0
0

0

0

0

0

0

0

0

0

1

0

0

0

0

0
..
.
0
0
0
..
.
0
0
0
..
.
0

0

0

0

0

0

0

0

0

0

0

0

0

0

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
0
0

1
0
0

1
0
0

1
0
0

0
0
0

0
0
0

0
0
0

0
0
0

0
1
0

0
0
1

1
0
0

0
0
0

0
0
0

1
0
0

0
0
0

0
0
0

1
0
0

0

0

0

1

0

0

0

0

0

0

0

0

1
















For this model we see that the “simplest moves”, i.e., the moves with the minimum
degree, correspond to the binomials of degree 4 such as
u111 u122 u212 u221 − u112 u121 u211 u222 ,

(23)

which is called a basic move. There are 9 such moves for the case of 3 × 3 × 3 tables.
Unfortunately, however, the set of these 9 moves does not become a Markov basis.
To see this consider the following example.
Example 2.6. Consider the 3 × 3 × 3 contingency tables with the fixed twodimensional marginals
(xij+ ) = (xi+k ) = (x+jk ) = (2, 1, 1, 1, 2, 1, 1, 1, 2)′.
22

(24)

We write a 3 × 3 × 3 table as follows.
x111 x112 x113
x121 x122 x123
x131 x132 x133

x211 x212 x213
x221 x222 x223
x231 x232 x233

x311 x312 x313
x321 x322 x323
x331 x332 x333

Then the fixed marginals (24) are displayed as
x111 x112 x113
x121 x122 x123
x131 x132 x133
2 1 1

2 x211 x212 x213
1 x221 x222 x223
1 x231 x232 x233
4 1 2 1

1 x311 x312 x313
2 x321 x322 x323
1 x331 x332 x333
4 1 1 2

1
1
2
4

2
1
1
4

1
2
1
4

1 4
1 4
2 4
4 12

where the rightmost table shows the marginal {x+jk }. There are 18 elements in this
fiber as follows.
200

1: 010
001
200

4: 010
001
110

7: 100
001
110

10 : 0 0 1
100
101

13 : 1 0 0
010
101

16 : 0 1 0
100

010
110
001

001
001
110

2: 010

200

001
110
010

010
001
101

5: 001

100
020
001

001
001
110

8: 100

100
020
001

001
100
011

11 : 0 0 1

010
011
100

100
010
002

14 : 0 1 0

010
101
010

100
010
002

17 : 1 0 0

001
200
010
110
001
110
100
101
100
011
100

010
101
010

001
010
101

3: 010

200

010
110
001

001
010
101

6: 001

100
011
010

001
010
101

9: 100

001
110
010

100
010
002

12 : 1 0 0

100
011
010

010
100
002

15 : 0 1 0

100
020
001

100
001
011

18 : 1 0 0

001
200
010
110
001
101
010
101
100
011
100

010
011
100

001
100
011

001
020
100

010
100
002

001
020
100

100
001
011

100
020
001

010
001
101

010
110
001

100
001
011

100
011
010

100
010
002

Now consider connecting these elements by the set of 9 basic moves such as (23).
The undirected graph we obtain is Figure 4. Because this is not connected, the set
of the basic moves is not a Markov basis. This example shows that we need moves
such as
u111 u122 u133 u213 u221 u232 − u113 u121 u132 u211 u222 u233
(25)
to constitute a Markov basis.
Now calculate a Markov basis by Macaulay2 for this example. Using a, b, c for
the sufficient statistics instead of v, the following is the commands to calculate a
reduced Gröbner basis for this problem.

23

4
3

5
6
1

9

7

2

8

15

17

10

16

18

12

14

13

11

Figure 4: Undirected graph obtained from the set of the basic moves.

R = QQ[a11,a12,a13,a21,a22,a23,a31,a32,a33,
b11,b12,b13,b21,b22,b23,b31,b32,b33,
c11,c12,c13,c21,c22,c23,c31,c32,c33,
x111,x112,x113,x121,x122,x123,x131,x132,x133,
x211,x212,x213,x221,x222,x223,x231,x232,x233,
x311,x312,x313,x321,x322,x323,x331,x332,x333,
MonomialOrder=>{27,27}]
I = ideal(x111-a11*b11*c11,x112-a11*b12*c12,x113-a11*b13*c13,
x121-a12*b11*c21,x122-a12*b12*c22,x123-a12*b13*c23,
x131-a13*b11*c31,x132-a13*b12*c32,x133-a13*b13*c33,
x211-a21*b21*c11,x212-a21*b22*c12,x213-a21*b23*c13,
x221-a22*b21*c21,x222-a22*b22*c22,x223-a22*b23*c23,
x231-a23*b21*c31,x232-a23*b22*c32,x233-a23*b23*c33,
x311-a31*b31*c11,x312-a31*b32*c12,x313-a31*b33*c13,
x321-a32*b31*c21,x322-a32*b32*c22,x323-a32*b33*c23,
x331-a33*b31*c31,x332-a33*b32*c32,x333-a33*b33*c33)
G = gb(I); g = gens(G)
selectInSubring(1,g)

Unfortunately, this calculation may be hard to carry out for average PC. In fact,
I could not finish the above calculation within one hour by my slow laptop (with
2.80 GHz CPU, 8.00 GB RAM, running on vmware). Instead, check the calculation
for 2 × 3 × 3 cases. With the similar input commands, we have the output instantly
in this case. From the output, we see that there are 1417 elements in the reduced
Gröbner basis of IA∗ , and 15 elements in the reduced Gröbner basis of IA as follows.
i10 : selectInSubring(1,g)
o10 = | x122x133x223x232-x123x132x222x233 x112x133x213x232-x113x132x212x233
------------------------------------------------------------------------x121x133x223x231-x123x131x221x233 x121x132x222x231-x122x131x221x232
------------------------------------------------------------------------x111x133x213x231-x113x131x211x233 x111x132x212x231-x112x131x211x232

24

------------------------------------------------------------------------x112x123x213x222-x113x122x212x223 x111x123x213x221-x113x121x211x223
------------------------------------------------------------------------x111x122x212x221-x112x121x211x222
------------------------------------------------------------------------x112x121x133x211x223x232-x111x123x132x212x221x233
------------------------------------------------------------------------x111x122x133x213x221x232-x113x121x132x211x222x233
------------------------------------------------------------------------x111x122x133x212x223x231-x112x123x131x211x222x233
------------------------------------------------------------------------x113x121x132x212x223x231-x112x123x131x213x221x232
------------------------------------------------------------------------x112x121x133x213x222x231-x113x122x131x212x221x233
------------------------------------------------------------------------x111x123x132x213x222x231-x113x122x131x211x223x232 |
1
15
o10 : Matrix R <--- R

We see that the set of the basic moves and the degree 6 moves such as (25)
actually constitutes a Markov basis for 2 × 3 × 3 cases.
The calculation for 3 × 3 × 3 cases can be carried out by faster software such as
4ti2 ([1]), that can be also used in Macaulay 2 as follows.
loadPackage "FourTiTwo"
A = matrix "1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1;
1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1;
1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0;
0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0;
0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0;
25

0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0;
0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0;
0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0;
0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0;
0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0;
0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1"
R = QQ[x111,x112,x113,x121,x122,x123,x131,x132,x133,
x211,x212,x213,x221,x222,x223,x231,x232,x233,
x311,x312,x313,x321,x322,x323,x331,x332,x333]
I = toricMarkov(A,R)
This calculation is finished within 1 second by my laptop. From the output, we see
that 27 basic moves such as (23) and 54 moves of degree 6 such as (25) constitute
a minimal Markov basis 2 . Using this minimal Markov basis, we can construct a
connected Markov chain for this fiber. The corresponding undirected graph is Figure
5.
4
3

5
6
1

9

7

2

8

15

17

10

16

18

12

14

13

11

Figure 5: Undirected graph obtained from a minimal Markov basis.

Interestingly, for the problems of the larger sizes, the structure of the Markov
basis becomes more complicated. For example, for the no three-factor interaction
model of 3 × 3 × 4 tables, the set of degree 4, 6, 8 moves becomes a Markov basis,
and for 3 × 3 × 5 tables, the set of degree 4, 6, 8, 10 moves becomes a Markov basis.
These results are summarized in Chapter 9 of [3].
2

The 4ti2 command toricMarkov gives a minimal Markov basis as the output. We can also
obtain a Gröbner basis by the command toricGroebner.

26

Acknowledgments
The author thanks the referee for his constructive suggestions that led to significant
improvement of this paper.

References
[1] 4ti2 team. 4ti2 – A software package for algebraic, geometric and combinatorial
problems on linear spaces. Available at www.4ti2.de.
[2] A. Agresti (1992). A survey of exact inference for contingency tables. Statistical
Science, 7, 131–177.
[3] S. Aoki, H. Hara and A. Takemura (2012). Markov bases in algebraic statistics.
Springer Series in Statistics.
[4] A. Agresti (2013). Categorical data analysis, 3rd ed. Wiley Series in Probability
and Statistics, Wiley.
[5] CoCoA Team. A system for doing computations in commutative algebra,
http://cocoa.dima.unige.it
[6] J. Cornfield (1956). A statistical problem arising from retrospective studies, in
Proceedings of 3rd Berkeley Symposium on Mathematical Statistics and Probability. University of California Press, Vol. 4, 135–148.
[7] D. Cox, J. Little and D. O’Shea. (2007). Ideals, Varieties, and Algorithms, 3rd
edn, New York, Springer.
[8] W. Decker, G. -M. Greuel, G. Pfister and H. Schönemann. Singular 3-1-2, A computer algebra system for polynomial computations,
http://www.singular.uni-kl.de
[9] P. Diaconis and B. Sturmfels (1998). Algebraic algorithms for sampling from
conditional distributions. Annals of Statistics, 26, 363–397.
[10] M. H. Gail and N. Mantel (1977). Counting the number of r × c contingency
tables with fixed margins. J. Am. Stat. Assoc., 72, 859–862.
[11] D. R. Grayson and M. E. Stillman. Macaulay2, a software system for research
in algebraic geometry, http://www.math.uiuc.edu/Macaulay2/
[12] S. J. Haberman (1988). A warning on the use of chi-squared statistics with
frequency tables with small expected cell counts. J. Am. Stat. Assoc. 83, 555–
560.

27

[13] W. K. Hastings (1970). Monte Carlo sampling methods using Markov chains
and their applications. Biometrika, 57, 97–109.
[14] T. Hibi (ed.) (2013). Grob̈ner bases, Statistics and software systems, Springer.
[15] E. L. Lehmann and J. P. Romano (2005). Testing statistical hypotheses, 3rd ed.
Springer Texts in Statistics, Springer.
[16] N.
Noro,
N.
Takayama,
H.
Nakayama,
and
K.
Ohara.
Risa/Asir:
A
computer
http://www.math.kobe-u.ac.jp/Asir/asir.html

K.
Kishiyama
algebra
system,

[17] L. Pachter and B. Sturmfels (eds.), (2005). Algebraic statistics for computational
biology, Cambridge University Press, New York.
[18] R. L. Plackett (1981). The analysis of categorical data, 2nd ed. Griffin, London.
[19] A. Vickers (2009). What is a p-value anyway? 34 Stories to help you actually
understand statistics. Addison-Wesley.

28

