Discrete Statistical Models with Rational Maximum Likelihood Estimator

TL;DR

This paper characterizes discrete statistical models with rational maximum likelihood estimators using real algebraic geometry, providing an algorithm to construct such models, including Bayesian networks and graphical models.

Contribution

It offers a complete characterization of models with rational MLEs and introduces an algorithm for constructing these models, connecting algebraic geometry with statistical modeling.

Findings

Identified all models with rational MLEs using algebraic geometry.

Developed an algorithm to construct models with rational MLEs.

Demonstrated the approach on Bayesian networks and graphical models.

Abstract

A discrete statistical model is a subset of a probability simplex. Its maximum likelihood estimator (MLE) is a retraction from that simplex onto the model. We characterize all models for which this retraction is a rational function. This is a contribution via real algebraic geometry which rests on results due to Huh and Kapranov on Horn uniformization. We present an algorithm for constructing models with rational MLE, and we demonstrate it on a range of instances. Our focus lies on models familiar to statisticians, like Bayesian networks, decomposable graphical models, and staged trees.

Equations

$$\begin{array}{c}{\hat{p}}_0=\frac{(u_0+u_1)(u_0+u_2)}{|u{|}^2},{\hat{p}}_1=\frac{(u_0+u_1)(u_1+u_3)}{|u{|}^2},\\ {\hat{p}}_2=\frac{(u_2+u_3)(u_0+u_2)}{|u{|}^2},{\hat{p}}_3=\frac{(u_2+u_3)(u_1+u_3)}{|u{|}^2}.\end{array}$$

$${\phi }_{(H,\lambda )}:{\mathbb{R}}_{>0}^{n+1}\to {\mathbb{R}}_{>0}^{n+1}.$$

$$\Phi ={\phi }_{(H,\lambda )}={\varphi }_{(\Delta ,\mathbf{m})}\circ H.$$

$${\Delta }_1\to {\Delta }_2,(s_0,s_1)\mapsto (s_0^2,s_0{s}_1,s_1)where{s}_0,s_1>0\; and{s}_0+s_1=1.$$

$${\hat{s}}_0=\frac{2u_0+u_1}{2u_0+2u_1+u_2}\text{and}{\hat{s}}_1=\frac{u_1+u_2}{2u_0+2u_1+u_2}.$$

$$\Phi(u_{0},u_{1},u_{2})\,\,=\,\,(\hat{s}_{0}^{2},\hat{s}_{0}\hat{s}_{1},\hat{s}_{1})\,\,=\,\,\,\small\biggl{(}\frac{(2u_{0}+u_{1})^{2}}{(2u_{0}{+}2u_{1}{+}u_{2})^{2}}\,,\,\frac{(2u_{0}{+}u_{1})(u_{1}{+}u_{2})}{(2u_{0}+2u_{1}+u_{2})^{2}}\,,\,\frac{u_{1}+u_{2}}{2u_{0}{+}2u_{1}{+}u_{2}}\biggr{)}.$$

$$H=\left(\begin{array}{ccc}2& 1& 0\\ 0& 1& 1\\ -2& -2& -1\end{array}\right)\mathrm{and}\lambda =(1,1,-1).$$

$$\Delta =x_3^2-x_1^2-x_1{x}_2+x_2{x}_3=(x_3-x_1)(x_1+x_2+x_3).$$

$$\phi_{(\Delta,{\bf m})}\,\,:\,\,(x_{1},x_{2},x_{3})\,\mapsto\,\biggl{(}\,\frac{x_{1}^{2}}{x_{3}^{2}}\,,\,\frac{x_{1}x_{2}}{x_{3}^{2}}\,\,,-\frac{x_{2}}{x_{3}}\biggr{)}.$$

$$(Hu)^{h_{j}}\,\,:=\,\,\prod_{i=1}^{m}\bigl{(}h_{i0}u_{0}+h_{i1}u_{1}+\cdots+h_{in}u_{n}\bigr{)}^{h_{ij}}\qquad{\rm for}\,\,j=0,1,\ldots,n.$$

$${\lambda }_0(Hu{)}^{h_0}+{\lambda }_1(Hu{)}^{h_1}+\cdots +{\lambda }_n(Hu{)}^{h_n}=1.$$

$$\mathbb{R}^{n+1}\,\to\,\mathbb{R}^{n+1},\,\,u\,\mapsto\,\bigl{(}\lambda_{0}(Hu)^{h_{0}},\,\lambda_{1}(Hu)^{h_{1}},\,\ldots,\,\lambda_{n}(Hu)^{h_{n}}\bigr{)}.$$

$$H^{\prime }=\left(\begin{array}{ccc}0& 2& 2\\ 2& 1& 0\\ 0& -1& -1\\ -2& -2& -1\end{array}\right)\mathrm{and}{\lambda }^{\prime }=\left(1,-\frac{1}{4},\frac{1}{4}\right)$$

$${\lambda }_j^{\prime }=\frac{{\lambda }_j{\mu }^{\mu \cdot h_{kj}}}{(1+\mu {)}^{(1+\mu )h_{kj}}}\text{for all j=0,…,n.}$$

$$H=\left(\begin{array}{cccc}-1& -1& -2& -2\\ 1& 0& 3& 2\\ 1& 3& 0& 2\\ -1& -2& -1& -2\end{array}\right).$$

$$\lambda\,=\,(\lambda_{0},\lambda_{1},\lambda_{2},\lambda_{3})\,=\,\small\biggl{(}\frac{2}{3}\,,\,-\frac{4}{27}\,,\,-\frac{4}{27}\,,\,\frac{1}{27}\biggr{)}.$$

$${\phi }_{(H,\lambda )}:{\mathbb{R}}_{>0}^4\to {\mathbb{R}}_{>0}^4,\left(\begin{array}{c}{u}_0\\ u_1\\ u_2\\ u_3\end{array}\right)\mapsto \left(\begin{array}{c}\frac{2(u_0+3u_2+2u_3)(u_0+3u_1+2u_3)}{3(u_0+u_1+2u_2+2u_3)(u_0+2u_1+u_2+2u_3)}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \frac{4(u_0+3u_1+2u_3{)}^3}{27(u_0+u_1+2u_2+2u_3)(u_0+2u_1+u_2+2u_3{)}^2}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \frac{4(u_0+3u_2+2u_3{)}^3}{27(u_0+u_1+2u_2+2u_3{)}^2(u_0+2u_1+u_2+2u_3)}\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \frac{(u_0+3u_2+2u_3{)}^2(u_0+3u_1+2u_3{)}^2}{27(u_0+u_1+2u_2+2u_3{)}^2(u_0+2u_1+u_2+2u_3{)}^2}\end{array}\right).$$

$$9p_1{p}_2-8p_0{p}_3=p_0^2-12p_3=0.$$

$$\gamma_{A}\,\,:\,(\mathbb{R}^{*})^{r}\to\mathbb{\mathbb{R}}\mathbb{P}^{m-1}\,,\,\,\,(t_{1},\ldots,t_{r})\,\mapsto\,\biggl{(}\,\prod_{i=1}^{r}t_{i}^{a_{i1}}:\,\prod_{i=1}^{r}t_{i}^{a_{i2}}:\,\,\cdots\,\,:\,\prod_{i=1}^{r}t_{i}^{a_{im}}\biggr{)}.$$

$$\bigl{\{}\,x\in(\mathbb{R}\mathbb{P}^{m-1})^{\vee}\,\mid\,\gamma_{A}^{-1}(H_{x}\cap Y_{A})\,\text{\rm{ is singular}}\,\bigr{\}}.$$

$$Y_{A,\sigma}^{*}\,\,=\,\,\bigl{\{}\,x\in Y_{A}^{*}\,:\,\sigma_{i}\cdot x_{i}>0\,\,\hbox{for}\,\,i=1,2,\ldots,m\,\bigr{\}}\quad\subset\,\,\,\mathbb{R}\mathbb{P}^{m-1}.$$

$$A=\left(\begin{array}{cccc}3& 2& 1& 0\\ 0& 1& 2& 3\end{array}\right).$$

$$Y_{A}\,=\,\overline{\bigl{\{}(t_{1}^{3}:t_{1}^{2}t_{2}:t_{1}t_{2}^{2}:t_{2}^{3})\in\mathbb{R}\mathbb{P}^{3}\,:\,t_{1},t_{2}\in{\mathbb{R}}^{*}\bigr{\}}}.$$

$${\Delta }_A=\underline{27x_1^2{x}_4^2}-18x_1{x}_2{x}_3{x}_4+4x_1{x}_3^3+4x_2^3{x}_4-x_2^2{x}_3^2.$$

$$\begin{array}{c}(t+1{)}^2(t+3),(t+5{)}^2(t-1),(t-1{)}^2(t+3),(t+5{)}^2(t-8),\\ (t-3{)}^2(t+1),(t-1{)}^2(t-3),\underline{(t-2{)}^2(t+3)},(t+1{)}^2(t-3).\end{array}$$

$$\mathbf{m}\,=\,x^{\mathrm{max}_{k}(h_{k}^{-})}\quad\text{ and }\quad\Delta\,=\,\mathbf{m}\cdot\biggl{(}1-\sum_{k=0}^{n}\lambda_{k}x^{h_{k}}\biggr{)}.$$

$$\phi_{(\Delta,{\bf m})}\,:\,(\mathbb{R}^{*})^{m}\rightarrow\mathbb{R}^{n+1},\,\,x\,\mapsto\,\bigl{(}\lambda_{0}x^{h_{0}},\,\lambda_{1}x^{h_{1}},\,\ldots,\,\lambda_{n}x^{h_{n}}\bigr{)}.$$

$$\phi_{(\Delta,{\bf m})}\,:\,(x_{1},x_{2},x_{3},x_{4})\,\mapsto\,\biggl{(}\frac{2}{3}\frac{x_{2}x_{3}}{x_{1}x_{4}},-\frac{4}{27}\frac{x_{3}^{3}}{x_{1}x_{4}^{2}},-\frac{4}{27}\frac{x_{2}^{3}}{x_{1}^{2}x_{4}},\frac{1}{27}\frac{x_{2}^{2}x_{3}^{2}}{x_{1}^{2}x_{4}^{2}}\biggr{)}.$$

$$u_{I,C}=u_{i+k+}=\sum _{j=1}^{r_2}\sum _{l=1}^{r_4}{u}_{ijkl}.$$

$${\hat{p}}_I=\frac{{\prod }_C{u}_{I,C}}{{\prod }_S{u}_{I,S}},$$

Full text

Maximum likelihood estimator,

graphical models,

algebraic statistics,

likelihood geometry,

real algebraic geometry,

discrete statistical models,

arXiv:1903.06110 \startlocaldefs

\endlocaldefs

,

and

t5Partial affiliation with Max-Planck-Institut für Mathematik in den Naturwissenschaften, Leipzig t2Corresponding author t3Partial affiliation with University of California, Berkeley

MSC 2010 Subject Classification: 62F10, 13P25, 14P10, 14M25.

1 Introduction

A discrete statistical model   is a subset $\mathcal{M}$ of the open probability simplex $\Delta_{n}$. Each point $p$ in $\Delta_{n}$ is a probability distribution on the finite state space $\{0,1,\ldots,n\}$, i.e. $p=(p_{0},p_{1},\ldots,p_{n}),$ where the $p_{i}$ are positive real numbers that satisfy $p_{0}+p_{1}+\cdots+p_{n}=1$. The model $\mathcal{M}$ is the set of all distributions $p\in\Delta_{n}$ that are relevant for an application.

In data analysis we are given an empirical distribution $u=(u_{0},u_{1},\ldots,u_{n})$. This is the point in $\Delta_{n}$ whose $i$th coordinate $u_{i}$ is the fraction of samples in state $i$. The maximum likelihood estimator (MLE) of $\mathcal{M}$ is a function $\Phi\colon\Delta_{n}\rightarrow\mathcal{M}$ that takes the empirical distribution $u$ to a distribution $\hat{p}=(\hat{p}_{0},\hat{p}_{1},\ldots,\hat{p}_{n})$ that best explains the given observations. Here “best” is understood in the sense of likelihood inference, so that $\hat{p}=\Phi(u)$ is the point in $\mathcal{M}$ that maximizes the log-likelihood function $p\mapsto\sum_{i=0}^{n}u_{i}\cdot{\rm log}(p_{i})$. For any vector $u$ in $\mathbb{R}^{n+1}_{>0}$, we set $\Phi(u):=\Phi(u/|u|)$ where $|u|=u_{0}+\cdots+u_{n}$.

Likelihood inference is consistent. This means that $\Phi(u)=u$ for $u\in\mathcal{M}$. This follows from the fact that the log-likelihood function is strictly concave on $\Delta_{n}$ and its unique maximizer is $p=u$. Hence, the MLE $\Phi$ is a retraction from the simplex onto the model.

This point is fundamental for two fields at the crossroads of mathematics and data science. Information Geometry [2] views the MLE as the nearest point map of a Riemannian metric on $\Delta_{n}$, given by the Kullback-Leibler divergence of probability distributions. Algebraic Statistics [6, 18] is concerned with models $\mathcal{M}$ whose MLE $\Phi$ is an algebraic function of $u$. This happens when the constraints that define $\mathcal{M}$ are given in terms of polynomials in $p$. In this article we address a question that is fundamental for both fields:

For which models $\mathcal{M}$ is the MLE $\,\Phi$ a rational function in the empirical distribution $u$?

The most basic example where the MLE is rational is the independence model for two binary random variables $(n=3)$. Here $\mathcal{M}$ is a surface in the tetrahedron $\Delta_{3}$. That surface is a familiar picture that serves as a point of entry for both Information Geometry and Algebraic Statistics. Points in $\mathcal{M}$ are positive rank one $2\times 2$ matrices $\begin{bmatrix}p_{0}&p_{1}\\ p_{2}&p_{3}\end{bmatrix}$ whose entries sum to one. The data takes the form of a nonnegative integer $2\times 2$ matrix $u$ of counts of observed frequencies. Hence $\,|u|=u_{0}{+}u_{1}{+}u_{2}{+}u_{3}$ is the sample size, and $u/|u|$ is the empirical distribution. The MLE $\hat{p}=\Phi(u)$ is evaluated by multiplying the row and column sums of $u$:

[TABLE]

These four expressions are rational, homogeneous of degree zero, and their sum is equal to $1$. See [11, Example 2] for a discussion of these formulas from our present perspective.

The surface $\mathcal{M}$ belongs to the class of graphical models [15]. Fix an undirected graph $G$ whose nodes represent random variables with finitely many states. The undirected graphical model $\mathcal{M}_{G}$ is a subset of $\Delta_{n}$, where $n{+}1$ is the number of states in the joint distribution. The graphical model $\mathcal{M}_{G}$ is decomposable if and only if the graph $G$ is chordal. Each coordinate $\hat{p}_{i}$ of its MLE is an alternating product of linear forms given by maximal cliques and minimal separators of $G$. A similar formula exists for directed graphical models, which are also known as Bayesian networks.

In both cases, the coordinates of the MLE are not only rational functions, but even alternating products of linear forms in $u=(u_{0},u_{1},\ldots,u_{n})$. This is no coincidence. Huh [11] proved that if $\Phi$ is a rational function then each of its coordinates is an alternating product of linear forms, with numerator and denominator of the same degree. Huh further showed that this alternating product must take a very specific shape. That shape was discovered by Kapranov [13], who named it the Horn uniformization. The results by Kapranov and Huh are valid for arbitrary complex algebraic varieties. They make no reference to a context where the coordinates are real, positive, and add up to $1$.

The present paper makes the leap from complex varieties back to statistical models. Building on the remarkable constructions by Kapranov and Huh, we here work in the setting of real algebraic geometry that is required for statistical applications. Our main result (Theorem 1) characterizes all models $\mathcal{M}$ in $\Delta_{n}$ whose MLE is a rational function. It is stated in Section 2 and all its ingredients are presented in a self-contained manner.

In Section 3 we examine models with rational MLE that are familiar to statisticians, such as decomposable graphical models and Bayesian networks. Our focus lies on staged tree models, a far-reaching generalization of discrete Bayesian networks, described in the book by Collazo, Görgen and Smith [4]. We explain how our main result applies to these models. The proof of Theorem 1 is presented in Section 4. This is the technical heart of our paper, building on the likelihood geometry of [12, §3]. We also discuss the connection to toric geometry and geometric modeling developed by Clarke and Cox [3]. In Section 5 we present our algorithm for constructing models with rational MLE, and we discuss its implementation and some experiments. The input is an integer matrix representing a toric variety, and the output is a list of models derived from that matrix. Our results suggest that only a very small fraction of Huh’s varieties in [11] are statistical models.

2 How to be Rational

Let $\mathcal{M}$ be a discrete statistical model in the open simplex $\Delta_{n}$ that has a well-defined maximum likelihood estimator $\Phi:\Delta_{n}\rightarrow\mathcal{M}$. We also write $\Phi:\mathbb{R}^{n+1}_{>0}\rightarrow\mathcal{M}$ for the induced map $u\mapsto\Phi(u/|u|)$ on positive vectors. If the $n+1$ coordinates of $\Phi$ are rational functions in $u$, then we say that $\mathcal{M}$ has rational MLE. The following is our main result.

Theorem 1**.**

The following are equivalent for the statistical model $\mathcal{M}$ with MLE $\Phi$:

• (1)

The model $\mathcal{M}$ has rational MLE.

• (2)

There exists a Horn pair $(H,\lambda)$ such that $\mathcal{M}$ is the image of the Horn map

[TABLE]

• (3)

There exists a discriminantal triple $(A,\Delta,{\bf m})$ such that $\mathcal{M}$ is the image under the monomial map $\phi_{(\Delta,{\bf m})}$ of precisely one orthant (9) of the dual toric variety $Y_{A}^{*}$.

The MLE of the model satisfies the following relation on the open orthant $\mathbb{R}^{n+1}_{>0}$:

[TABLE]

This theorem matters for statistics because it reveals when a model has an MLE of the simplest possible closed form. Property (2) says that the polynomials appearing in the numerators and denominators of the rational formulas must factor into linear forms with positive coefficients. Property (3) offers a recipe, based on toric geometry, for explicitly constructing such models. The advance over [11] is that Theorem 1 deals with positive real numbers. It hence furnishes the definitive solution in the case of applied interest.

The goal of this section is to define all the terms seen in parts (2) and (3) of Theorem 1.

Example 2**.**

We first discuss Theorem 1 for a simple experiment: Flip a biased coin. If it shows heads, flip it again. This is the model with $n=2$ given by the tree diagram

s_{0}$$s_{1}$$s_{0}$$s_{1}$$p_{0}$$p_{1}$$p_{2}.

The model $\mathcal{M}$ is a curve in the probability triangle $\Delta_{2}$. The tree shows its parametrization

[TABLE]

The implicit representation of the curve $\mathcal{M}$ is the equation $p_{0}p_{2}-(p_{0}+p_{1})p_{1}=0$. Let $(u_{0},u_{1},u_{2})$ be the counts from repeated experiments. A total of $2u_{0}+2u_{1}+u_{2}$ coin tosses were made. We estimate the parameters as the empirical frequency of heads resp. tails:

[TABLE]

The MLE is the retraction from the triangle $\Delta_{2}$ to the curve $\mathcal{M}$ given by the formula

[TABLE]

Hence $\mathcal{M}$ has rational MLE. We see that the Horn pair from part (2) in Theorem 1 has

[TABLE]

We next exhibit the discriminantal triple $(A,\Delta,\bf m)$ in part (3) of Theorem 1. The matrix $A=\begin{pmatrix}1&1&1\end{pmatrix}$ gives a basis of the left kernel of $H$. The second entry is the polynomial

[TABLE]

The third entry marks the leading term ${\bf m}=x_{3}^{2}$. These data define the monomial map

[TABLE]

The toric variety of the matrix $A$ is the point $Y_{A}=\{(1:1:1)\}$ in $\mathbb{P}^{2}$. Our polynomial $\Delta$ vanishes on the line $Y_{A}^{*}=\{x_{1}+x_{2}+x_{3}=0\}$ that is dual to $Y_{A}$. The relevant orthant is the open line segment $Y^{*}_{A,\sigma}\coloneqq\{(x_{1}:x_{2}:x_{3})\in Y_{A}^{*}\,:\,x_{1},x_{2}>0\,\,{\rm and}\,\,x_{3}<0\}$. Part (3) in Theorem 1 says that $\mathcal{M}$ is the image of $Y^{*}_{A,\sigma}$ under $\phi_{(\Delta,{\bf m})}$. The MLE is $\Phi=\phi_{(\Delta,\bf m)}\circ H$.

We now come to the definitions needed for Theorem 1. Let $H=(h_{ij})$ be an $m\times(n{+}1)$ integer matrix whose columns sum to zero, i.e. $\,\sum_{i=1}^{m}h_{ij}=0$ for $j=0,\ldots,n$. We call such a matrix a Horn matrix and denote its columns by $h_{0},h_{1},\ldots,h_{n}$. The following alternating products of linear forms are rational functions of degree zero:

[TABLE]

We use the notation $v^{h}\coloneqq\prod_{i}v_{i}^{h_{i}}$ for two vectors $v,h$ of the same size. The Horn matrix $H$ is friendly if there exists a real vector $\,\lambda=(\lambda_{0},\ldots,\lambda_{n})$ with $\lambda_{i}\neq 0$ for all $i$ such that the following identity holds in the rational function field $\mathbb{R}(u_{0},u_{1},\ldots,u_{n})$:

[TABLE]

If this holds, then we call $(H,\lambda)$ a friendly pair, and we consider the rational function

[TABLE]

The friendly pair $(H,\lambda)$ is called a Horn pair if the function (4) is defined for all positive vectors, and it maps these to positive vectors. If these conditions hold then we write $\,\varphi_{(H,\lambda)}:\mathbb{R}^{n+1}_{>0}\,\to\,\mathbb{R}^{n+1}_{>0}\,$ for the restriction of (4) to the positive orthant. We call $\,\varphi_{(H,\lambda)}\,$ the Horn map associated to the Horn pair $(H,\lambda)$.

The difference between our Horn pairs and the more general pairs considered by Huh in [11] is the positivity condition we just introduced, along with the “friendliness” condition. These conditions guarantee that the image of the Horn map lies in the probability simplex, which is necessary for its interpretation as a statistical model. They also imply special properties for the Horn pair, see Propositions 22 and 23 in Section 4. The examples in Section 5 show that only a fraction of Huh’s pairs $(H,\lambda)$ are Horn pairs.

Different Horn pairs may give rise to the same Horn map. For example, the Horn pair

[TABLE]

also gives the map in Example 2. This is because the first and third rows of $H^{\prime}$ are collinear, causing the cancellation of linear factors in the Horn map. Following [3], a Horn pair $(H,\lambda)$ is minimal if the matrix $H$ has no zero rows and no pair of collinear rows.

Lemma 3**.**

Let $(H^{\prime},\lambda^{\prime})$ be a Horn pair arising from the Horn pair $(H,\lambda)$ by replacing two collinear rows $r_{k}$ and $r_{\ell}$ in $H$ such that $r_{\ell}=\mu r_{k}$ with their sum $r_{k}+r_{\ell}$ and setting

[TABLE]

Then the Horn maps $\varphi_{(H^{\prime},\lambda^{\prime})}$ and $\varphi_{(H,\lambda)}$ are equal.

Proof.

Let $w_{k}$ and $w_{\ell}$ be the linear forms associated to the rows $r_{k}$ and $r_{\ell}$ respectively. Fix a column index $j$. We have $w_{\ell}=\mu w_{k}$ and $h_{\ell j}=\mu h_{kj}$. The factors of the $j$-th coordinates of the Horn maps $\varphi_{(H,\lambda)}$ and $\varphi_{(H^{\prime},\lambda^{\prime})}$ that have changed after the operation are $\lambda_{j}w_{k}^{b_{kj}}w_{\ell}^{b_{\ell j}}=\lambda_{j}\mu^{\mu\cdot h_{kj}}w_{k}^{(1+\mu)h_{kj}}$ for $(H,\lambda)$ and $\lambda^{\prime}_{j}(w_{k}+w_{\ell})^{(1+\mu)h_{kj}}=\lambda^{\prime}_{j}(1+\mu)^{(1+\mu)h_{kj}}w_{k}^{(1+\mu)h_{kj}}$ for $(H^{\prime},\lambda^{\prime})$. Equating these two gives the desired formula. ∎

Every Horn map is represented by a unique minimal Horn pair. This follows by unique factorization, see also [3, Proposition 6.11]. To make a Horn pair minimal, while retaining the Horn map, we can use Lemma 3 repeatedly, deleting zero rows as they appear.

Example 4**.**

We illustrate the equivalence of (1) and (2) in Theorem 1 for the model described in [12, Example 3.11]. Here $n=3$ and $m=4$ and the Horn matrix equals

[TABLE]

This Horn matrix is friendly because the following vector satisfies the identity (3):

[TABLE]

The pair $(H,\lambda)$ is a Horn pair, with associated Horn map

[TABLE]

Indeed, this rational function takes positive vectors to positive vectors. The image of the map $\varphi_{(H,\lambda)}$ is a subset $\mathcal{M}$ of the tetrahedron $\Delta_{3}=\{p\in\mathbb{R}^{4}_{>0}:p_{0}+p_{1}+p_{2}+p_{3}=1\}$. We regard the subset $\mathcal{M}$ as a discrete statistical model on the state space $\{0,1,2,3\}$. The model $\mathcal{M}$ is the curve of degree $4$ inside $\Delta_{3}$ defined by the two quadratic equations

[TABLE]

As in [12, Example 3.11], one verifies that $\mathcal{M}$ has rational MLE, namely $\,\Phi=\varphi_{(H,\lambda)}$.

We next define all the terms used in part (3) of Theorem 1. Fix a matrix $A=(a_{ij})\in\mathbb{Z}^{r\times m}$ of rank $r$ that has the vector $(1,\ldots,1)$ in its row span. The connection to part (2) of Theorem 1 will be that the rows of $A$ span the left kernel of $H$. We identify the columns of $A$ with Laurent monomials in $r$ unknowns $t_{1},\ldots,t_{r}$. The associated monomial map is

[TABLE]

Here $\mathbb{R}^{*}=\mathbb{R}\backslash\{0\}$ and $\mathbb{R}\mathbb{P}^{m-1}$ denotes the real projective space of dimension $m-1$. Let $Y_{A}$ be the closure of the image of $\gamma_{A}$. This is the projective toric variety given by $A$.

Every point $x=(x_{1}:\cdots:x_{m})$ in the dual projective space $(\mathbb{R}\mathbb{P}^{m-1})^{\vee}$ corresponds to a hyperplane $H_{x}$ in $\mathbb{R}\mathbb{P}^{m-1}$. The dual variety $Y_{A}^{*}$ to the toric variety $Y_{A}$ is the closure of

[TABLE]

Here, the term singular means that the variety $\gamma_{A}^{-1}(H_{x}\cap Y_{A})$ has a singular point in $\mathbb{(}\mathbb{R}^{*})^{r}$. A general point $x$ in $Y_{A}^{*}$ hence corresponds to a hyperplane $H_{x}$ that is tangent to the toric variety $Y_{A}$ at a point $\gamma_{A}(t)$ with nonzero coordinates. We identify sign vectors $\sigma\in\{-1,+1\}^{m}$ with orthants in $\mathbb{R}^{m}$. These map in a $2$-to-$1$ manner to orthants in $\mathbb{R}\mathbb{P}^{m-1}$. If we intersect them with $Y_{A}^{*}$, then we get the orthants of the dual toric variety:

[TABLE]

One of these is the distinguished orthant in Theorem 1, part (3).

Example 5**.**

Fix $m=4$ and $r=2$. The following matrix has $(1,1,1,1)$ in its row span:

[TABLE]

As in [12, Example 3.9], the toric variety of $A$ is the twisted cubic curve in $3$-space:

[TABLE]

The dual toric variety $Y_{A}^{*}$ is a surface in $(\mathbb{R}\mathbb{P}^{3})^{\vee}$. Its points $x$ represent planes in $\mathbb{R}\mathbb{P}^{3}$ that are tangent to the curve $Y_{A}$. Such a tangent plane corresponds to a cubic $\,x_{1}t^{3}+x_{2}t^{2}+x_{3}t+x_{4}\,$ with a double root. Just as we recognize quadrics with a double root by the vanishing of the quadratic discriminant, a cubic with coefficients $(x_{1},x_{2},x_{3},x_{4})$ has a double root if and only if the following discriminant vanishes:

[TABLE]

Hence, $Y_{A}^{*}$ is the surface of degree $4$ in $(\mathbb{R}\mathbb{P}^{3})^{\vee}$ defined by $\Delta_{A}$. All eight orthants $Y_{A,\sigma}^{*}$ are non-empty. The coefficient vectors of the following eight cubics lie on different orthants:

[TABLE]

For instance, the underlined cubic corresponds to the point $\,x=(1,-1,-8,12)$ in the orthant $Y_{A,\sigma}^{*}$ associated with the sign vector $\sigma=(+1,-1,-1,+1)$.

Let $\Delta$ be a homogeneous polynomial in $m$ variables with $n+2$ monomials and $\mathbf{m}$ one of these monomials. There is a one-to-one correspondence between such pairs $(\Delta,\mathbf{m})$ and pairs $(H,\lambda)$ where $H$ is a Horn matrix of size $m\times(n+1)$ and $\lambda$ is a coefficient vector. Namely, for $k=0,\dotsc,n$ write $h_{k}^{+}$ resp. $h_{k}^{-}$ for the positive resp. negative part of the column vector $h_{k}$, so that $h_{k}=h_{k}^{+}-h_{k}^{-}$. In addition, let $\mathrm{max}_{k}(h_{k}^{-})$ be the entrywise maximum of the $h^{-}_{k}$. We pass from pairs $(H,\lambda)$ to pairs $(\Delta,\mathbf{m})$ as follows:

[TABLE]

For the converse, from pairs $(\Delta,\mathbf{m})$ to pairs $(H,\lambda)$, we divide $\Delta$ by $\mathbf{m}$ and use the same equations to determine the pair $(H,\lambda)$. Note that the polynomial $\Delta$ being homogeneous and the matrix $H$ being a Horn matrix are equivalent conditions using the equations (12). Given a pair $(\Delta,\mathbf{m})$ with associated pair $(H,\lambda)$, we define the monomial map

[TABLE]

We now present the definition that is needed for part (3) of Theorem 1.

Definition 6**.**

A discriminantal triple $(A,\Delta,{\bf m})$ consists of

an $r\times m$ integer matrix $A$ of rank $r$ having $(1,1,\ldots,1)$ in its row span, 2. 2.

an $A$-homogeneous polynomial $\Delta$ that vanishes on the dual toric variety $Y_{A}^{*}$, 3. 3.

a distinguished term ${\bf m}$ among those that occur in the polynomial $\Delta$,

such that the pair $(H,\lambda)$ associated to $(\Delta,\mathbf{m})$ is a Horn pair. Here, the polynomial $\Delta$ being $A$-homogeneous means that $Av=Aw$ for any two exponent vectors $v$ and $w$ of $\Delta$.

All definitions are now complete. We illustrate Definition 6 for our running example:

Example 7**.**

Let $A$ be the $2\times 4$ matrix in (10), $\Delta=\Delta_{A}$ its discriminant in (11), and ${\bf m}=27x_{1}^{2}x_{4}^{2}$ the special term. Then $(A,\Delta,{\bf m})$ is a discriminantal triple with associated sign vector $\sigma=(+1,-1,-1,+1)$. The orthant $Y_{A,\sigma}^{*}$, highlighted in Example 5, is a semialgebraic surface in $Y_{A}^{*}\subset\mathbb{R}\mathbb{P}^{3}$. This surface is mapped into the tetrahedron $\Delta_{3}$ by

[TABLE]

The image of this map is a curve in $\Delta_{3}$, namely the model $\mathcal{M}$ in Example 4. We verify (1) by comparing (7) with (13). The former is obtained from the latter by setting $x=Hu$.

3 Staged Trees

We consider contingency tables $u=(u_{i_{1}i_{2}\cdots i_{m}})$ of format $r_{1}\times r_{2}\times\cdots\times r_{m}$. Following [6, 15], these represent joint distributions of discrete statistical models with $n+1=r_{1}r_{2}\cdots r_{m}$ states. Namely, the contingency table $u$ represents the probability distribution $p\coloneqq u/|u|.$ For any subset $C\subset\{1,\ldots,m\}$, one considers the marginal table $u_{C}$ that is obtained by summing out all indices not in $C$. The entries of the marginal table $u_{C}$ are sums of entries in $u$. To obtain the entry $u_{I,C}$ of $u_{C}$ for any state $I=(i_{1},i_{2},\ldots,i_{m}),$ we fix the indices of the states in $C$ and sum over the indices not in $C$. For example, if $m=4$, $C=\{1,3\}$, $I=(i,j,k,l)$, then $u_{C}$ is the $r_{1}\times r_{3}$ matrix with entries

[TABLE]

Such linear forms are the basic building blocks for familiar models with rational MLE.

Consider an undirected graph $G$ with vertex set $\{1,\ldots,m\}$ which is assumed to be chordal. The associated decomposable graphical model $\mathcal{M}_{G}$ in $\Delta_{n}$ has the rational MLE

[TABLE]

where the product in the numerator is over all maximal cliques $C$ of $G$, and the product in the denominator is over all separators $S$ in a junction tree for $G$. See [15, §4.4.1]. We shall regard $G$ as a directed graph, with edge directions given by a perfect elimination ordering on the vertex set $\{1,\ldots,m\}$. This turns $\mathcal{M}_{G}$ into a Bayesian network. More generally, a Bayesian network $\mathcal{M}_{G}$ is given by a directed acyclic graph $G$. We write ${\rm pa}(j)$ for the set of parents of the node $j$. The model $\mathcal{M}_{G}$ in $\Delta_{n}$ has the rational MLE

[TABLE]

If $G$ comes from an undirected chordal graph then (14) arises from (15) by cancellations.

Example 8** ($m=4$).**

We revisit two examples from on page 36 in [6, §2.1]. The star graph $G=[14][24][34]$ is chordal. The MLE for $\mathcal{M}_{G}$ is the map $\Phi$ with coordinates

[TABLE]

The left expression is (14). The right is (15) for the directed graph $1\rightarrow 4$, $4\rightarrow 2$, $4\rightarrow 3$.

The chain graph $G=[12][23][34]$ is chordal. Its MLE is the map $\Phi$ with coordinates

[TABLE]

This is the Horn map given by the matrix $H$ in Figure 1 and $\lambda=(1,\dotsc,1).$

The formulas (14) and (15) are familiar to statisticians. Theorem 1 places them into a larger context. However, some readers may find our approach too algebraic and too general. Our aim in this section is to lay out a useful middle ground: staged tree models.

Staged trees were introduced by Smith and Anderson [17] as a generalization of discrete Bayesian networks. They furnish an intuitive representation of many situations that the above graphs $G$ cannot capture. In spite of their wide scope, staged tree models are appealing because of their intuitive formalism for encoding events. For an introduction see the textbook [4]. In what follows we study parts (1) and (2) in Theorem 1 for staged trees.

To define a staged tree model, we consider a directed rooted tree $\mathcal{T}$ with at least two edges emanating from each non-leaf vertex, a label set $S=\{s_{i}\mid i\in I\}$, and a labeling $\theta\colon\operatorname{E}(\mathcal{T})\to S$ of the edges of the tree. Each vertex of $\mathcal{T}$ has a corresponding floret, which is the multiset of edge labels emanating from it. The labeled tree $\mathcal{T}$ is a staged tree if any two florets are either equal or disjoint. Two vertices in $\mathcal{T}$ are in the same stage if their corresponding florets are the same. From now on, $F$ denotes the set of florets of $\mathcal{T}$.

Definition 9**.**

Let $J$ be the set of root-to-leaf paths in the tree $\mathcal{T}$. We set $|J|=n+1$. For $i\in I$ and $j\in J$, let $\mu_{ij}$ denote the number of times edge label $s_{i}$ appears in the $j$-th root-to-leaf path. The staged tree model $\mathcal{M}_{\mathcal{T}}$ is the image of the parametrization

[TABLE]

where the parameter space is \,\Theta:=\bigl{\{}(s_{i})_{i\in I}\in(0,1)^{|I|}:\sum_{s_{i}\in f}s_{i}=1\text{ for all florets f\in F}\bigr{\}}, and $p_{j}=\prod_{i\in I}s_{i}^{\mu_{ij}}$ is the product of the edge parameters on the $j$-th root-to-leaf path.

In the model $\mathcal{M}_{\mathcal{T}}$, the tree $\mathcal{T}$ represents possible sequences of events. The parameter $s_{i}$ associated to an edge $vv^{\prime}$ is the transition probability from $v$ to $v^{\prime}$. All parameter labels in a floret sum to $1$. The fact that distinct nodes in $\mathcal{T}$ can have the same floret of parameter labels enables staged tree models to encode conditional independence statements [17]. This allows us to represent any discrete Bayesian network or decomposable model as a staged tree model. Our first staged tree was seen in Example 2. Here is another specimen.

Example 10** ($n=15$).**

Consider the decomposable model for binary variables given by the $4$-chain $G=[12][23][34]$ as in Example 8. Figure 1 shows a realization of $\mathcal{M}_{G}$ as a staged tree model $\mathcal{M}_{\mathcal{T}}$. The leaves of $\mathcal{T}$ represent the outcome space $\{0,1\}^{4}$. Nodes with the same color have the same associated floret. The blank nodes all have different florets. The seven florets of $\mathcal{T}$ are

[TABLE]

Next we show that staged tree models have rational MLE, so they satisfy part (1) of Theorem 1. Our formula for $\Phi$ uses the notation for $I,J$ and $\mu_{ij}$ introduced in Definition 9. This formula is known in the literature on chain event graphs (see e.g. [16]).

Proposition 11**.**

Let $\mathcal{M}_{\mathcal{T}}$ be a staged tree model, and let $u=(u_{j})_{j\in J}$ be a vector of counts. For $i\in I$, let $f$ be the floret containing the label $s_{i}$, and define the estimates

[TABLE]

The rational function $\,\Phi\,$ that sends $\,(u_{j})_{j\in J}\,$ to $\,(\hat{p}_{j})_{j\in J}\,$ is the MLE of the model $\mathcal{M}_{\mathcal{T}}$.

Proof.

We prove that the likelihood function $L(p,u)$ has a unique maximum at $p=(\hat{p}_{j})_{j\in J}$. For a floret $f\in F$, we fix the vector of parameters $s_{f}=(s_{i})_{s_{i}\in f}$, and we define the local likelihood function $L_{f}(s_{f},u)=\prod_{s_{i}\in f}s_{i}^{\alpha_{i}}$, where $\alpha_{i}=\sum_{j}\mu_{ij}u_{j}$. We have

[TABLE]

Since the $L_{f}$ depend on disjoint sets of unknowns, maximizing $L$ is achieved by maximizing the factors $L_{f}$ separately. But $L_{f}$ is the likelihood function of the full model $\Delta_{|f|-1}$, given the data vector $(\alpha_{i})_{s_{i}\in f}$. The MLE of that model is $\hat{s}_{i}=\alpha_{i}/\sum_{s_{\ell}\in f}\alpha_{\ell}$, where $s_{i}\in f$. We conclude that \,\mathrm{argmax}_{s_{f}}\bigl{(}L_{f}(s_{f},u)\bigr{)}=(\hat{s}_{i})_{s_{i}\in f}\, and \,\mathrm{argmax}_{p}\bigl{(}L(p,u)\bigr{)}=(\hat{p}_{j})_{j\in J}. ∎

Remark 12**.**

Here is a method for evaluating the MLE in Proposition 11. Let $[v]\subset J$ be the set of root-to-leaf paths through a node $v$ in the tree $\mathcal{T}$ and define $u_{[v]}=\sum_{j\in[v]}u_{j}$. The ratio $\frac{u_{[v^{\prime}]}}{u_{[v]}}$ is the empirical transition probability from $v$ to $v^{\prime}$ given arrival at $v$. To obtain $\hat{s}_{i}$ we first compute the quotients $\frac{u_{[v^{\prime}]}}{u_{[v]}}$ for all edges $vv^{\prime}$ with parameter label $s_{i}$. We aggregate them by adding their numerators and denominators separately. This gives $\,\hat{s}_{i}=(\sum u_{[v^{\prime}]})/(\sum u_{[v]})$, where both sums range over all edges $vv^{\prime}$ with parameter label $s_{i}$.

Proposition 11 yields an explicit description of the Horn pair $(H,\lambda)$ associated to $\mathcal{M}_{\mathcal{T}}$.

Corollary 13**.**

Fix a staged tree model $\mathcal{M}_{\mathcal{T}}$ as above. Let $H$ be the $(|I|+|F|)\times|J|$ matrix whose rows are indexed by the set $I\sqcup F$ and entries are given by

[TABLE]

Define $\lambda\in\{-1,+1\}^{|J|}$ by $\lambda_{j}=(-1)^{\sum_{f}h_{fj}}$. Then $(H,\lambda)$ is a Horn pair for $\mathcal{M}_{\mathcal{T}}$.

Given a staged tree $\mathcal{T}$, we call the matrix $H$ in Corollary 13 the Horn matrix of $\mathcal{T}$.

Remark 14**.**

In Corollary 13, for a floret $f$, let $H_{f}$ be the submatrix of $H$ with row indices $\{i:s_{i}\in f\}\cup\{f\}$. Then $H$ is the vertical concatenation of the matrices $H_{f}$ for $f\in F$.

Example 15**.**

For the tree $\mathcal{T}$ in Example 10, the Horn matrix $H$ of $\mathcal{M}_{\mathcal{T}}$ is given in Figure 1. The vector $\lambda$ of the Horn pair $(H,\lambda)$ is the vector of ones $(1,\ldots,1)\in\mathbb{R}^{16}$. The rows of $H$ are indexed by the florets and labels

[TABLE]

Note that $(H,\lambda)$ is not minimal. Following the recipe in Lemma 3, we can delete the rows ${s_{0},s_{1},f_{2},f_{3}}$ of the matrix $H$ by summing the pairs $(s_{0},f_{2})$ and $(s_{1},f_{3})$ and deleting zero rows. The result is the minimal Horn pair $(H^{\prime},\lambda^{\prime})$, where $\lambda^{\prime}=(-1,\dotsc,-1)$.

Two staged trees $\mathcal{T}$ and $\mathcal{T}^{\prime}$ are called statistically equivalent in [9] if there exists a bijection between the sets of root-to-leaf paths of $\mathcal{T}$ and $\mathcal{T}^{\prime}$ such that, after applying this bijection, $\mathcal{M}_{\mathcal{T}}=\mathcal{M}_{\mathcal{T}^{\prime}}$ in the open simplex $\Delta_{n}$. A staged tree model may have different but statistically equivalent tree representations. In [9, Theorem 1], it is shown that statistical equivalence of staged trees can be determined by a sequence of operations on the trees, named swap and resize. One of the advantages of describing a staged tree model via its Horn pair is that it gives a new criterion to decide whether two staged trees are statistically equivalent. This is simpler to implement than the criterion given in [9].

Corollary 16**.**

Two staged trees are statistically equivalent if and only if their associated Horn pairs reduce to the same minimal Horn pair.

One natural operation on a staged tree $\mathcal{T}$ is identifying two florets of the same size. This gives a new staged tree $\mathcal{T}^{\prime}$ whose Horn matrix is easy to get from that of $\mathcal{T}$.

Corollary 17**.**

Let $\mathcal{T}^{\prime}$ be a staged tree arising from $\mathcal{T}$ by identifying two florets $f$ and $f^{\prime}$, say by the bijection $(-)^{\prime}\colon f\to f^{\prime}$. The Horn matrix $H^{\prime}$ of $\mathcal{M}_{\mathcal{T}^{\prime}}$ arises from the Horn matrix $H$ of $\mathcal{M}_{\mathcal{T}}$ by replacing the blocks $H_{f}$ and $H_{f^{\prime}}$ in $H$ by the block $H^{\prime}_{f}$ defined by

[TABLE]

Proof.

This follows from the definition of the Horn matrices for $\mathcal{M}_{\mathcal{T}}$ and $\mathcal{M}_{\mathcal{T}^{\prime}}$. ∎

Example 18**.**

Let $\mathcal{T}^{\prime}$ be the tree obtained from Example 10 by identifying florets $f_{4}$ and $f_{5}$ in $\mathcal{T}$. Then $\mathcal{M}_{\mathcal{T}^{\prime}}$ is the independence model of two random variables with four states.

Now we turn to part (3) of Theorem 1. We describe the triple $(A,\Delta,{\bf m})$ for a staged tree model $\mathcal{M}_{\mathcal{T}}$. The pair $(H,\lambda)$ was given in Corollary 13. Let $A$ be any matrix whose rows span the left kernel of $H$, set $m=|I|+|F|$, and write $s$ for the $m$-tuple of parameters $(s_{i},s_{f})_{i\in I,f\in F}$. From the Horn matrix in Corollary 13 we see that

[TABLE]

where $f$ depends on $i$, $\,{\bf m}=\operatorname{lcm}(\prod_{i}s_{f}^{\mu_{ij}}:f\in F)\,$ and $\,\epsilon_{j}={\sum_{i}\mu_{ij}}$. The sign vector $\sigma$ for the triple $(A,\Delta,\mathbf{m})$ is given by $\sigma_{i}=+1$ for $i\in I$ and $\sigma_{f}=-1$ for $f\in F$. Then $Y_{A,\sigma}^{*}$ gets mapped to $\mathcal{M}_{\mathcal{T}}$ via $\phi_{(\Delta,{\bf m})}$. Moreover, the map $\phi_{\mathcal{T}}$ from Definition 9 factors through $\phi_{(\Delta,{\bf m})}$. Indeed, if we define $\iota:\Theta\to Y_{A,\sigma}^{*}$ by $(s_{i})_{i\in I}\mapsto(s_{i},-1)_{i\in I,f\in F}$, then $\phi_{\mathcal{T}}=\phi_{(\Delta,{\bf m})}\circ\iota$. The following derivation is an extension of that in [12, Example 3.13].

Example 19**.**

Let $\mathcal{M}_{\mathcal{T}}$ be the $4$-chain model in Example 10. Here the discriminant is

[TABLE]

Our notation for the parameters matches the row labels of the Horn matrix $H$ in Figure 1. This polynomial of degree $7$ is irreducible, so it equals the $A$-discriminant: $\,\Delta=\Delta_{A}$. The underlying matrix $A$ has format $13\times 21$, and we represent it by its associated toric ideal

[TABLE]

The toric variety $Y_{A}=\mathcal{V}(I_{A})$ has dimension $12$ and degree $141$. It lives in a linear space of codimension $2$ in $\mathbb{P}^{20}$, where it is defined by eight cubics and eight quartics. The dual variety $Y_{A}^{*}=\mathcal{V}(\Delta_{A})$ is the above hypersurface of degree seven. We have ${\bf m}=f_{1}f_{2}f_{3}f_{4}f_{5}f_{6}f_{7}$, and $\sigma$ is the vector in $\{-1,+1\}^{21}$ that has entry $+1$ at the indices corresponding to the $s_{i}$ and entry $-1$ at the indices corresponding to the $f_{i}$.

It would be interesting to study the combinatorics of discriminantal triples for staged tree models. Our computations suggest that, for many such models, the polynomial $\Delta$ is irreducible and equals the $A$-discriminant $\Delta_{A}$ of the underlying configuration $A$. However, this is not true for all staged trees, as seen in equation (2) of Example 2. We close this section with a familiar class of models with rational MLE whose associated $\Delta$ factor.

Example 20**.**

The multinomial distribution encodes the experiment of rolling a $k$-sided die $m$ times and recording the number of times one observed the $j$-th side, for $j=1,\dotsc,k$. The associated model $\mathcal{M}$ is the independence model for $m$ identically distributed random variables on $k$ states. We have $n+1=\binom{k+m-1}{m}$. The Horn matrix $H$ is the $(k+1)\times(n+1)$ matrix whose columns are the vectors $(-m,i_{1},i_{2},\ldots,i_{k})^{T}$ where $i_{1},i_{2},\ldots,i_{k}$ are nonnegative integers whose sum equals $m$. Here, $\,A=(1\,\,1\,\cdots\,1)$, so the $A$-discriminant equals $\,\Delta_{A}=x_{0}+x_{1}+\cdots+x_{k}$. The following polynomial is a multiple of $\Delta_{A}$:

[TABLE]

This $\Delta$, with its marked term ${\bf m}=(-x_{0})^{m}$, encodes the MLE for the model $\mathcal{M}$:

[TABLE]

Here, $I$ ranges over all vectors in $\mathbb{N}^{k}$ that sum to $m$, and $I_{j}$ denotes the $j$-th entry of $I$.

4 Proof of the Main Theorem

In this section we prove Theorem 1. For a pair $(H,\lambda)$ consisting of a Horn matrix $H$ and a coefficient vector $\lambda$, let $\varphi$ be the rational map defined in (4). We use $\varphi$ and $\varphi_{(H,\lambda)}$ interchangeably in this section, as well as $\phi$ and $\phi_{(\Delta,\mathbf{m})}$. Recall that its $j$-th coordinate is

[TABLE]

For a fixed data vector $u\in\mathbb{N}^{n+1}$, we define the likelihood function for the image of $\varphi$:

[TABLE]

Lemma 21**.**

Let $H=(h_{ij})$ be a Horn matrix, $\lambda$ a vector satisfying (3) and $u\in\mathbb{N}^{n+1}$. Then $u$ is a critical point of its own likelihood function $L_{u}$. Furthermore, if $u^{\prime}$ is another critical point of $L_{u}$, then $\varphi(u)=\varphi(u^{\prime})$.

Proof.

We compute the partial derivatives of $L_{u}$. For $\ell=0,\dotsc,n$ we find

[TABLE]

For $v=u$, this evaluates to zero, since the sums in the fraction cancel and the $\ell$-th column of $H$ sums to zero. This shows that $u$ is a critical point.

Next, let $u^{\prime}$ be another critical point of $L_{u}$. Using terminology from [11, Theorem 1], this means that $\varphi(u^{\prime})$ is a critical point of the likelihood function $L(p,u)$ of the model $\mathcal{M}$ defined as the image of $\varphi$. The same holds for $\varphi(u)$. By the implication (ii) to (i) in [11, Theorem 1], the model $\mathcal{M}$ has ML degree one. This implies $\varphi(u)=\varphi(u^{\prime})$. ∎

We use [11] to explain the relation between models with rational MLE and Horn pairs.

Proof of Theorem 1, Equivalence of (1) and (2).

Let $\mathcal{M}$ be a model with rational MLE $\Phi$. The Zariski closure of $\mathcal{M}$ is a variety whose likelihood function has a unique critical point. By [11, Theorem 1], there is a Horn matrix $H$ and a coefficient vector $\lambda$ such that $\varphi_{(H,\lambda)}=\Phi$. Now, the required sum-to-one and positivity conditions for $\varphi_{(H,\lambda)}$ are satisfied because they are satisfied by the MLE $\Phi$. Indeed, the MLE of any discrete statistical model maps positive vectors $u$ in $\mathbb{R}^{n+1}_{>0}$ into the simplex $\Delta_{n}$. Conversely, we claim that every Horn pair $(H,\lambda)$ specifies a nonempty model $\mathcal{M}$ with rational MLE. Indeed, define $\mathcal{M}$ to be the image of $\varphi_{(H,\lambda)}$. By the defining properties of the Horn pair, we have $\mathcal{M}\subset\Delta_{n}$. Lemma 21 shows that $\varphi_{(H,\lambda)}$ is the MLE of $\mathcal{M}$. ∎

Next, we relate Horn pairs to discriminantal triples.

Proof of Theorem 1, Equivalence of (2) and (3).

We already exhibited a bijection between pairs $(H,\lambda)$ and pairs $(\Delta,\mathbf{m})$ given by Equation 12. The matrix $A$ is the left kernel of $H$ and forms the triple $(A,\Delta,\mathbf{m})$. It is a matrix of size $r\times m$ of rank $r$. When $H$ is a Horn matrix, $A$ contains $(1,\dotsc,1)$ in its row span. This implies that the polynomial $\Delta$ is homogeneous, which in turn implies that it is $A$-homogeneous by $AH=0$.

Next, we show that the pair $(H,\lambda)$ being friendly corresponds to the polynomial $\Delta$ vanishing on $Y_{A}^{*}$. This is part of the desired equivalence.

Claim**.**

The pair $(H,\lambda)$ is friendly if and only if the $A$-homogeneous polynomial $\Delta$ vanishes on the dual toric variety $Y_{A}^{*}$.

Proof of Claim.

Let $(H,\lambda)$ be friendly and $A$ as above. The Laurent polynomial $\,q:=\Delta/{\mathbf{m}}\,$ is a rational function on $\mathbb{P}^{m-1}$ that vanishes on the dual toric variety $Y_{A}^{*}$. To see this, consider the exponentiation map $\,\varphi_{2}\,\colon\,\mathbb{P}^{m-1}\to\mathbb{R}^{n+1},\,x\mapsto\lambda*x^{H}$, where $*$ is the entrywise product and $x^{H}\coloneqq(x^{h_{0}},\dotsc,x^{h_{n}})$. Let $f=1-(p_{0}+\cdots+p_{n})$. We have $q=f\circ\varphi_{2}$. By [11, Theorems 1 and 2], the function $\varphi_{2}$ maps an open dense subset of $Y_{A}^{*}$ dominantly to the closure $\overline{\mathcal{M}}$ of the image of $\varphi_{(H,\lambda)}$. Since $f=0$ on $\overline{\mathcal{M}}$, we have $f\circ\varphi_{2}=0$ on an open dense subset of $Y_{A}^{*}$, hence $q=0$ on $Y_{A}^{*}$, so $\Delta=0$ there as well.

Conversely, let $\Delta$ vanish on $Y_{A}^{*}$. We claim that $q(x)$ is zero for all $x=Hu$ in the image of the linear map $H$. We may assume $\mathbf{m}(x)\neq 0$. We only need to show that $x$ is in the dual toric variety $Y_{A}^{*}$, since $\Delta$ vanishes on it. So, let $x_{i}=\sum_{j=0}^{n}h_{ij}u_{j}$ for $i=1,\dotsc m$. We claim that $t=(1,\dotsc,1)$ is a singular point of the hypersurface

[TABLE]

First, the point $t$ lies on that hypersurface since the columns of $H$ sum to zero:

[TABLE]

For $s=1,\dotsc,r$ we have $\frac{\partial}{\partial t_{s}}t^{a_{i}}=a_{si}t^{a_{i}-e_{s}}$, with $e_{s}$ the standard basis vector of $\mathbb{Z}^{r}$, and

[TABLE]

This is zero at $t=(1,\dotsc,1)$ because $AH=0$. ∎

We now prove the rest of the equivalence. Let $(H,\lambda)$ be a Horn pair, let $\varphi$ be its Horn map and let $\phi$ be the associated monomial map. Let $\mathcal{M}$ be the statistical model with MLE $\varphi$, so $\mathcal{M}=\varphi(\mathbb{R}_{>0}^{n+1})$. We have $\varphi=\phi\circ H$. By Proposition 23, there exists a unique sign vector $\sigma$ such that $\operatorname{im}H|_{\mathbb{R}^{n+1}_{>0}}\subseteq\mathbb{R}^{m}_{\sigma}$. From the proof of the above claim we know that $\operatorname{im}H\subseteq Y_{A}^{*}$. Together, we have

[TABLE]

By [11, Theorems 1 and 2] we have $\phi(Y_{A}^{*})\subseteq\mathcal{M}^{\prime}$, where $\mathcal{M}^{\prime}$ is the real part of $\overline{\varphi(\mathbb{C}^{n+1})}$. We also have $\phi(Y_{A,\sigma}^{*})\subseteq\mathbb{R}^{n+1}_{>0}$ by definition of the orthant. Thus $\phi(Y_{A,\sigma}^{*})\subseteq\mathcal{M}^{\prime}\cap\mathbb{R}^{n+1}_{>0}$. Every element in the latter set is a fixed point of the rational function $\varphi$, by a similar argument as in Lemma 21 for complex space. Hence $\mathcal{M}^{\prime}\cap\mathbb{R}_{>0}^{n+1}=\mathcal{M}$, so $\phi(Y_{A,\sigma}^{*})\subseteq\mathcal{M}$.

Finally, if $(A,\Delta,\mathbf{m})$ is a discriminantal triple then $(H,\lambda)$ is a Horn pair by definition. This completes the proof of Theorem 1. ∎

In the next two propositions, we formulate simple criteria to decide whether the image of the map $\varphi_{(H,\lambda)}$ associated to a Horn matrix $H$ and a coefficient vector $\lambda$ is a statistical model. These are essential for constructing models with rational MLE in Algorithm 1.

Proposition 22**.**

Let $(H,\lambda)$ be a friendly pair. If there exists a vector $u_{0}\in\mathbb{R}^{n+1}$ such that $\varphi(u_{0})>0$, then we have $\varphi(u)>0$ for all $u$ in $\mathbb{R}_{>0}^{n+1}$ where it is defined.

Proof.

The function $\varphi$ is homogeneous of degree zero. It suffices to prove each coordinate of $\varphi(u)$ is a positive real number, for all vectors $u$ with positive integer entries. Indeed, every positive $u$ in $\mathbb{R}^{n+1}$ can be approximated by rational vectors, which can be scaled to be integral. The open subset $U=\varphi^{-1}(\Delta_{n})$ of $\mathbb{R}^{n+1}$ contains $u_{0}$ by our assumptions. If $U=\mathbb{R}^{n+1}$, then we are done. Else, $U$ has a nonempty boundary $\partial U$. By continuity, $\partial U\subseteq\varphi^{-1}(\partial\Delta_{n})$. The likelihood function $L_{u}$ for the data vector $u$ vanishes on $\partial U$.

We claim that $L_{u}$ has a critical point in $U$. The closed subset $\overline{U}$ is homogeneous. Seen in projective space $\mathbb{P}^{n}$, it becomes compact. The likelihood function $L_{u}$ is well defined on this compact set in $\mathbb{P}^{n}$, since it is homogeneous of degree zero, and $L_{u}$ vanishes on the boundary. Hence the restriction $L_{u}|_{U}$ is either identically zero or it has a critical point in $U$. But, since $u_{0}\in U$ is a point with $L_{u}(u_{0})\neq 0$, the second statement must be true.

Pick such a critical point $u^{\prime}$. Since $U$ is open in $\mathbb{R}^{n+1}$, the point $u^{\prime}$ is also critical point of $L_{u}$. By Lemma 21 and since $u^{\prime}\in U$, we have $\varphi(u)=\varphi(u^{\prime})>0$. ∎

Proposition 23**.**

Let $(H,\lambda)$ be a friendly pair, with no zero or collinear rows in $H$. Then $(H,\lambda)$ is a Horn pair if and only if for every row $r_{i}$ of $H$ all nonzero entries of $r_{i}$ have the same sign $\sigma_{i}$, and the sign vector $\sigma=(\sigma_{i})$ satisfies $\lambda_{j}\sigma^{h_{j}}>0$ for all columns $j$.

Proof.

Let $(H,\lambda)$ be a Horn pair. Let $\ell_{1},\dotsc,\ell_{k}$ be the linear forms corresponding to the rows in $H$ that have both positive and negative entries. Since $\ell_{1}$ has positive and negative coefficients, there exists a positive vector $u$ such that $\ell_{1}(u)=0$. Since $(H,\lambda)$ is minimal, we may choose $u>0$ such that $\ell_{1}(u)=0$ but $\ell_{k^{\prime}}(u)\neq 0$ for all $k^{\prime}\neq 1.$ The form $\ell_{1}$ appears in the numerator of some coordinate of $\varphi$, making this coordinate zero at $u$. But this contradicts the fact that $(H,\lambda)$ is a Horn pair. Therefore we cannot have rows with both positive and negative entries. The inequalities $\lambda_{j}\sigma^{h_{j}}>0$ then follow from the definition of a Horn pair by evaluating $\varphi(u)$ for some positive vector $u$.

Conversely, if the sign vector $\sigma$ is well-defined, the inequalities $\lambda_{j}\sigma^{h_{j}}>0$ imply that $\varphi(u)>0$ for all positive $u$. Hence $(H,\lambda)$ is a Horn pair. ∎

Every model with rational MLE arises from a toric variety $Y_{A}$. In some cases, the model is itself a toric variety $Y_{C}$. It is crucial to distinguish the two matrices $A$ and $C$. The two toric structures are very different. For instance, every undirected graphical model is toric [6, Proposition 3.3.3]. The toric varieties $Y_{C}$ among staged tree models $\mathcal{M}_{\mathcal{T}}$ were classified in [5]. The $4$-chain model $\mathcal{M}_{\mathcal{T}}=Y_{C}$ is itself a toric variety of dimension $7$ in $\mathbb{P}^{15}$. But it arises from a toric variety $Y_{A}$ of dimension $12$ in $\mathbb{P}^{20}$, seen in Example 19.

Toric models with rational MLE play an important role in geometric modeling [3, 7]. Given a matrix $C\in\mathbb{Z}^{r\times(n+1)}$ and a vector of weights $w\in\mathbb{R}^{n+1}_{>0}$, one considers the scaled projective toric variety $Y_{C,w}$ in $\mathbb{R}\mathbb{P}^{n}$. This is defined as the closure of the image of

[TABLE]

The set $\mathcal{M}_{C,w}$ of positive points in $Y_{C,w}$ is a statistical model in $\Delta_{n}$. There is a natural homeomorphism from the toric model $\mathcal{M}_{C,w}$ onto the polytope of $C$. This is known in geometry as the moment map. For a reference from algebraic statistics, see [6, Proposition 2.1.5]. In geometric modeling the pair $(C,w)$ defines toric blending functions [14].

It is desirable for the toric blending functions to have rational linear precision [3, 14]. The property is rare and it depends in a subtle way on $(C,w)$. Garcia-Puente and Sottile [7] established the connection to algebraic statistics. They showed that rational linear precision holds for $(C,w)$ if and only if the statistical model $\mathcal{M}_{C,w}$ has rational MLE.

Example 24**.**

The most classical blending functions with rational linear precision live on the triangle $\{x\in\mathbb{R}^{3}_{>0}:x_{1}{+}x_{2}{+}x_{3}=1\}$. They are the Bernstein basis polynomials

[TABLE]

Here $C$ is the $3\times\binom{m+1}{2}$ matrix whose columns are the vectors $(i,j,m-i-j)$. The weights are $w_{(i,j)}=\frac{m!}{i!j!(m-i-j)!}$. The toric model $\mathcal{M}_{C,w}$ is the multinomial family, where (19) is the probability of observing $i$ times $1$, $j$ times $2$ and $m-i-j$ times $3$ in $m$ trials. This model has rational MLE, as seen in Example 20. Again, notice the distinction between the two toric varieties. Here, $Y_{A}$ is a point in $\mathbb{P}^{m}$, whereas $Y_{C}$ is a surface in $\mathbb{P}^{\binom{m}{2}-1}$.

Clarke and Cox [3] raise the problem of characterizing all pairs $(C,w)$ with rational linear precision. This was solved by Duarte and Görgen [5] for pairs arising from staged trees. While the problem remains open in general, our theory in this paper offers new tools. We may ask for a characterization of discriminantal triples whose models are toric.

5 Constructing Models with Rational MLE

Part (3) in Theorem 1 allows us to construct models with rational MLE starting from a matrix $A$ that defines a projective toric variety $Y_{A}$. To carry out this construction effectively we propose Algorithm 1. In most cases, the dual variety $Y_{A}^{*}$ is a hypersurface, and we can compute its defining polynomial $\Delta_{A}$, the discriminant [8]. The polynomial $\Delta$ in a discriminantal triple can be any homogeneous multiple of $\Delta_{A}$, but we just take $\Delta=\Delta_{A}$. For all terms $\mathbf{m}$ in $\Delta_{A}$, we check whether $(A,\Delta_{A},\mathbf{m})$ is a discriminantal triple. We implemented this algorithm in Macaulay2, and our code is available online at [19].

Lines 1 and 15 of Algorithm 1 are computations with Gröbner bases. Executing Line 15 can be very slow. It may be omitted if one is satisfied with obtaining the parametric description and MLE $\Phi^{(\ell)}$ of the model $\mathcal{M}_{\ell}$. For the check in Line 14, we rely on Proposition 22 for correctness. A check based on the criterion in Proposition 23 is also possible.

Example 25** ($r=2,m=4$).**

For distinct integers $\alpha,\beta,\gamma>0$ with $\gcd(\alpha,\beta,\gamma)=1$ let

[TABLE]

We ran Algorithm 1 for all $613$ such matrices with $0<\alpha<\beta<\gamma\leq 17$. Line 1 computes the discriminant $\Delta_{A}$ of the univariate polynomial $f(t)=x_{1}+x_{2}t^{\alpha}+x_{3}t^{\beta}+x_{4}t^{\gamma}$. The number $n+2$ of terms of these discriminants equals $7927/613=12.93$ on average. Thus a total of $7927$ candidate triples $(A,\Delta_{A},{\bf m})$ were tested in Lines 12 to 21. Precisely 123 of these were found to be discriminantal triples. This is a fraction of 1.55 %. Hence, only 1.55 % of the resulting complex varieties permitted by [11] are actually statistical models.

Here is a typical model that was discovered. Take $\alpha=1,\beta=4,\gamma=7$. The discriminant

[TABLE]

has $9$ terms, so $n=7$. The term ${\bf m}$ is underlined. The associated model is a curve of degree ten in $\Delta_{7}$. Its prime ideal $I^{(\ell)}$ is generated by $18$ quadrics. Among them are $15$ binomials that define a toric surface of degree six: $49p_{1}p_{2}-48p_{0}p_{3},3p_{0}p_{4}-p_{2}^{2},\ldots,361p_{3}p_{7}-128p_{5}^{2}$. Inside that surface, our curve is cut out by three quadrics, like $\,26068p_{2}^{2}+73728p_{0}p_{5}$ $+703836p_{0}p_{6}+234612p_{2}p_{6}+78204p_{4}p_{6}+612864p_{0}p_{7}+212268p_{2}p_{7}+78204p_{4}p_{7}-8379p_{7}^{2}$.

Example 26** ($r=3,m=6$).**

For any positive integers $\alpha,\beta,\gamma,\varepsilon$, we consider the matrix

[TABLE]

The discriminant $\Delta_{A}$ is the resultant of two trinomials $\,x_{1}+x_{2}t^{\alpha}+x_{3}t^{\beta}\,$ and $\,x_{4}+x_{5}t^{\gamma}+x_{6}t^{\varepsilon}$. We ran Algorithm 1 for all 138 such matrices with $0<\alpha<\beta\leq 17,\,0<\gamma<\varepsilon\leq 17,\,\gcd(\alpha,\beta)=\gcd(\gamma,\varepsilon)=1$. The number $n+2$ of terms of these discriminants equals 2665/138 = 19.31 on average. Thus a total of 2665 candidate triples $(A,\Delta_{A},{\bf m})$ were tested in Line 13. Precisely 93 of these are discriminantal triples. This is only 3.49 %.

We now shift gears by looking at polynomials $\Delta$ that are multiples of the $A$-discriminant.

Example 27** ($r=1,m=4$).**

We saw in Examples 2 and 20 that interesting models arise from the matrix $A=(1\ 1\ \cdots\ 1)$ whose toric variety is just a point. Any homogeneous multiple $\Delta$ of the linear form $\Delta_{A}=x_{1}+x_{2}+\cdots+x_{m}$ can be the input in Line 1 of Algorithm 1. Here, taking $\Delta=\Delta_{A}$ results in the model given by the full simplex $\Delta_{m-2}$.

Let $m=4$ and abbreviate $x^{a}=x_{1}^{a_{1}}x_{2}^{a_{2}}x_{3}^{a_{3}}x_{4}^{a_{4}}$ and $|a|=a_{1}{+}a_{2}{+}a_{3}{+}a_{4}$ for $a\in\mathbb{N}^{4}$. We conducted experiments with two families of multiples. The first uses binomial multipliers:

[TABLE]

where $|a|=|b|\in\{1,2,\dots,8\}$ and $\gcd(x^{a},x^{b})=1$. This gives $1028$ polynomials $\Delta$. The numbers of polynomials of degree $2,3,4,5,6,7,8,9$ is $6,\,21,\,46,\,81,\,126,\,181,\,246,\,321$. For the second family we use the trinomial multiples

[TABLE]

where $|a|{=}|b|{=}|c|\!\in\!\{1,2,3\}$ and $\gcd(x^{a}\!,x^{b}\!,x^{c})=1$. Each list contains $4$ quadrics, $104$ cubics and $684$ quartics. We report our findings in a table:

[TABLE]

All $12$ Horn pairs in the first family represent the same model, up to permuting coordinates. All are coming from the six quadrics of the family. The model is the surface in $\Delta_{4}$ defined by the $2\times 2$ minors of the matrix $\begin{pmatrix}p_{0}&p_{1}&p_{2}\,\,\\ p_{0}{+}p_{1}{+}p_{2}&p_{3}&p_{4}\,\,\end{pmatrix}$. This is a staged tree model similar to Example 2, but now with three choices at each blue node instead of two. The eight Horn pairs in the third family represent two distinct models. Four of the eight Horn pairs represent a surface in $\Delta_{5}$ and the rest represent a surface in $\Delta_{6}.$

Our construction of models with rational MLE starts with families where $r$ and $m$ are fixed. However, as the entries of the matrix $A$ go up, the number $n+1$ of states increases. This suggests the possibility of listing all models for fixed values of $n$. Is this list finite?

Problem**.**

Suppose that $n$ is fixed. Are there only finitely many models with rational MLE in the simplex $\Delta_{n}$? Can we find absolute bounds, depending only on $n$, for the dimension, degree and number of ideal generators of the associated varieties in $\mathbb{P}^{n}$?

Algorithm 1 is a tool for studying these questions experimentally. At present, we do not have any clear answers, even for $n=3$, where the models are curves in a triangle.

Acknowledgements

The first author was supported by the Deutsche Forschungsgemeinschaft DFG under grant 314838170, GRK 2297 MathCoRe.