Motivated by previous work on moment varieties for Gaussian distributions and their mixtures, we study moment varieties for two other statistically important two-parameter distributions: the inverse Gaussian and gamma distributions. In particular, we realize the moment varieties as determinantal varieties and find their degrees and singularities. We also provide computational evidence for algebraic identifiability of mixtures, and study the identifiability degree and Euclidean distance degree.

§ INTRODUCTION

Suppose $X$ is a univariate random variable from a distribution that depends on finitely many parameters $\theta=(\theta_1,\ldots,\theta_n)$, and we want to determine $\theta$ from observed data. One approach to this problem is the method of moments. The $r$th moment of $X$ is defined as $m_r(\theta):=\mathbb{E}(X^r)$, and all the moments assemble into a moment generating function
\[M(t,\theta):=\EE(e^{tX})=\sum_{r=0}^\infty\frac{m_r(\theta)t^r}{r!}.\]

For many distributions, the moments $m_r(\theta)$ are polynomial in the parameters $\theta$ (see, e.g., Table 2 of [7]). This makes it possible to estimate $\theta$ by computing the sample moments $\widehat{m}_r:=\frac{1}{N}\sum_{i=1}^N x_i^r$ for some large number $N$ observations, and then solving the system

\begin{equation}\label{eq:sample_moments}
m_r(\theta)=\widehat{m}_r,\quad r\in[d]
\end{equation}

for some $d\in\mathbb{N}$. The law of large numbers gives that we can expect arbitrarily good approximations of $\theta$ among the solutions of this system, if $N$ is large enough.

The starting point for the algebraic-geometric perspective on the method of moments is the following construction. For a random variable $X$ as above, the $d$th moment variety $\mathcal{M}_d\subset\mathbb{P}^d$ is defined as the Zariski closure of the image of the map

\[\CC^n\to\mathbb{P}^d,\quad \theta\mapsto [m_0(\theta):m_1(\theta):\cdots:m_d(\theta)].\]

It is typically the case that if $d$ is large enough, the map $\CC^n\to\mathcal{M}_d$ is generically one-to-one (meaning that the system (<ref>) has a unique complex solution for generic $\widehat{m}\in\M_d$), and we say that we have rational identifiability. A weaker condition for the map is algebraic identifiability, meaning that it is generically finite-to-one (in the sense that the system (<ref>) has finitely many complex solutions for generic $\widehat{m}\in\M_d$). By the theorem of the dimension of fibers, we have algebraic identifiability if and only if $\dim(\mathcal{M}_{d})=n$.

An important tool for studying moment varieties is the concept of cumulants. The $r$th cumulant $k_r(\theta)$ of $X$ is defined by the cumulant generating function

\[K(t):=\log(M(t))=\sum_{r=1}^\infty\frac{k_r(\theta)t^r}{r!},\]

and the $d$th cumulant variety $\K_d$ is defined as the Zariski closure of the image of the map

\[\CC^n\to\CC^d,\quad \theta\mapsto (k_1(\theta),\ldots,k_d(\theta)).\]

The relation $K(t) = \log(M(t))$ gives rise to a nonlinear change of coordinates
\[\CC^d\to\CC^d,\quad(m_1,\ldots,m_d)\mapsto (m_1, m_2-m_1^2, m_3-3m_2m_1+2m_1^3,\ldots)\]

that expresses moments of order $\leq d$ as polynomials in cumulants of order $\leq d$, and vice versa. It restricts to an isomorphism of affine varieties $\M_d\cap\{m_0\neq 0\}\to\K_d$, and hence a birational map $\M_d\dashrightarrow\K_d$, and is an example of a Cremona tranformation in the language of [11]. Passing from moments to cumulants often gives rise to simpler varieties (see, e.g., <cit.> and the more general discussion in <cit.>), and cumulants also have the additional advantage of being additive over independent variables. However, moments behave simpler with respect to certain other statistical operations; in particular, moment varieties of mixtures are secant varieties, which is a fact that we will return to in <ref>. Furthermore, the cumulant variety only captures an open dense patch of the moment variety, and therefore might miss important geometric aspects of the moment variety such as singularities.

§.§ Previous work on moment varieties

The algebraic study of moments and cumulant varieties goes back to Pearson in 1894 [23], who studied the mixtures of two univariate Gaussians. The case of Gaussians and Gaussian mixtures was later revisited from a nonlinear algebra perspective in [2], as well as the subsequent papers [5, 1, 20] that addressed various identifiability questions, and (in the latter case) approximation methods based on numerical algebraic geometry. In [2], the authors study moments as projective varieties. The main advantage of studying moments from the perspective of projective geometry is that it allows for the study of secant varieties, which correspond to the moment varieties for mixtures of distributions. In general, algebraic and even rational identifiability for distributions listed in <cit.> is clear, but it is unknown for mixtures of these distributions.  Other distributions for which moment and cumulant varieties have been studied include uniform distributions on polytopes [19], Dirac and Pareto distributions [17], and mixtures of products [3].

§.§ Main contributions

In this paper, we study the moment varieties for the inverse Gaussian and gamma distributions, which are both two-parameter distributions used in a wide array of applications. Our main results (Theorems <ref> and <ref>) give the homogeneous prime ideals for these moment varieties. In particular, we show that these ideals are determinantal ideals, meaning that they are generated by all $t$-minors of a matrix $H$; such ideals are denoted $I_t(H)$. Interestingly, both these varieties turn out to be Cohen–Macaulay. This raises Problem <ref> on Cohen–Macaulayness of moment varieties.

Using these results, we find the degrees and Hilbert series for these moment varieties, and give Gröbner bases for their ideals. In particular, Corollary <ref> addresses a conjecture stated in [16] on the Hilbert series for the moment variety of the gamma distribution. We also identify the singular loci of the moment varieties (Propositions <ref> and <ref>).

Additionally, we study the exponential and chi-squared distributions, which are one-parameter specializations of the gamma distribution, and we show that their moment varieties are rational normal curves (Propositions <ref> and <ref>).

Using the moments-to-cumulants transformation, we prove that the cumulant varieties of the inverse Gaussian are (scaled) Veronese varieties (Proposition <ref>); this has previously been shown to be the case for the gamma distribution <cit.>. These simple geometric models, together with previous results on the projective geometry of the corresponding moment varieties, could be the starting point for addressing the question of identifiability of mixtures.

We conclude the paper by providing computational evidence in <ref> towards a conjecture on
algebraic identifiability for mixtures of inverse Gaussian distributions and mixtures of gamma distributions (Conjecture <ref>). We also report on numerical estimations of the identifiability degree (<ref>), and compute some Euclidean distance degrees (<ref>).

§.§ Applications of commutative algebra techniques

Moment varieties and other varieties arising in algebraic statistics provide interesting yet concrete algebraic structures to which one can apply techniques from commutative algebra. In this paper, the most important techniques are from the theories of determinantal ideals and Cohen–Macaulay rings, and we provide specific references to the necessary results as needed. We refer the reader to [6, 9, 13] for further general background on these topics. In particular, the Cohen–Macaulay property is key in many of our results. Cohen–Macaulay varieties can be seen as “mildly singular” varieties with good homological properties, and for which intersection multiplicity behaves well. For instance, they are equidimensional varieties with no embedded components. We also use polarization and Stanley–Reisner theory from combinatorial commutative algebra in <ref>. We hope that our proof techniques inspire further applications of commutative algebra to algebraic statistics.

§.§ Outline of the paper

The rest of the paper is organized as follows. In <ref>, we study moment varieties for the inverse Gaussian distribution. In <ref>, we study moment varieties for the gamma distribution, as well as for certain one-parameter specializations. In <ref> of the paper, we outline a number of interesting directions for future work on the inverse Gaussian and gamma distributions. In particular, we provide evidence towards algebraic identifiability for mixtures of these distributions. We then take a more computational perspective on moment varieties, and discuss questions regarding identifiability and Euclidean distance degrees.

Code for the computational experiments presented in this paper is publically available at

§.§ Acknowledgements

The authors thank Bernd Sturmfels for suggesting the problem and his guidance, and Carlos Améndola, Oliver Clarke, Francesco Galuppi, Alexandros Grosdos Koutsoumpelia, Lizzie Pratt, Kristian Ranestad, Jose Israel Rodriguez, Bernd Sturmfels, and Matteo Varbaro for helpful discussions. Part of this research was performed while the authors were visiting the Institute for Mathematical and Statistical Innovation (IMSI) for the Algebraic Statistics long program, which was supported by the National Science Foundation (NSF grant DMS-1929348), and the authors thank the other participants of the program for their feedback and comments. OH was partially supported by the Novo Nordisk project with grant reference number NNF20OC0065582. TY was partially supported by NSF grant DGE-1841052.

declare function=pdf_ig(,μ,λ)=sqrt(λ/(2*pi*^3)*exp(-λ*(-μ)^2/(2*μ^2*));
axis lines=left,
column sep=0.2em,
row sep=0.2em,
legend cell align=left,
legend entries=
$\mu=1{,}\:\:\: \lambda=5$,
$\mu=2{,}\:\:\: \lambda=5$,
$\mu=9{,}\:\:\: \lambda=2$
[smooth, domain=0:5, red]pdf_ig(x,1,5);
[smooth, domain=0:5, blue]pdf_ig(x,2,5);
[smooth, domain=0:5, teal]pdf_ig(x,1,2);

Inverse Gaussian distribution

declare function=gamma_function()=(2.506628274631*sqrt(1/) + 0.20888568*(1/)^(1.5) + 0.00870357*(1/)^(2.5) - (174.2106599*(1/)^(3.5))/25920 - (715.6423511*(1/)^(4.5))/1244160)*exp((-ln(1/)-1)*);,
declare function=pdf_gamma(,,̨θ) = ^(-̨1)*exp(-/θ)/(θ^*̨gamma_function()̨);
axis lines=left,
column sep=0.2em,
legend cell align=left,
legend entries=
$k=1{,}\:\:\: \theta=1$,
$k=2{,}\:\:\: \theta=1$,
$k=1{,}\:\:\: \theta=2$
[smooth, domain=0:5, red]pdf_gamma(x,1,1);
[smooth, domain=0:5, blue]pdf_gamma(x,2,1);
[smooth, domain=0:5, teal]pdf_gamma(x,1,2);

Gamma distribution

Probability density functions for the inverse Gaussian and gamma distribution for varying parameter values.

§ INVERSE GAUSSIAN MOMENT VARIETIES

We begin with the study of moment varieties of the inverse Gaussian distribution, which is also sometimes called the Wald distribution. We show that the ideals of its moment varieties are determinantal. Based on this, we use results from the theory of determinantal ideals to find degrees and Gröbner bases, and we compute the singular loci of the moment varieties.

Despite being less known than the classical Gaussian distribution, the inverse Gaussian has significant applications in modeling different phenomena. While the Gaussian distribution is widely known for its symmetry and versatility in modeling various data, the inverse Gaussian distribution offers a valuable alternative for situations where asymmetry and right-skewed tails are more appropriate (see Figure <ref>), such as in modeling lifetime phenomena.

As in the case of the Gaussian distribution, the inverse Gaussian is defined by two parameters, namely $\mu$ and $\lambda$. However, its density function is supported on $(0,+\infty)$ and it is given by

\[ \textstyle f_{\mu,\lambda}(x)={\sqrt {\frac {\lambda }{2\pi x^{3}}}}\,\exp {\biggl (}-{\frac {\lambda (x-\mu )^{2}}{2\mu ^{2}x}}{\biggr )}.\]

Here, $\mu$ represents the mean while $ \lambda$ is called shape paramter since it affects how peaked and right-tailed the density function is. Unlike the Gaussian distribution, $\lambda$ does not directly coincide with the variance, though these two quantities are related by the following formula:

The main motivation behind the definition of this distribution comes from Brownian motion. In this context, the inverse Gaussian works as a dual for the Gaussian distribution, whence its name. While the Gaussian models the Brownian motion's level at a fixed time, the inverse Gaussian models the first passage time of the Brownian motion to a certain fixed level.

The moment generating function of the inverse Gaussian is given by

\[M(t)=\exp\left(\tfrac{\lambda}{\mu}\left(1-\sqrt{1- \tfrac{2\mu^2t}{\lambda}}\right)\right).\]

In particular, writing down the first $d$ moments yields

\begin{equation}\label{eq:parametrization_inverse_gaussian}
\begin{array}{rl}
m_0 &= 1,\\
m_1&= \mu,\\
m_2&= \frac{\mu^2(\lambda+\mu)}{\lambda},\\
m_3&= \frac{\mu^3(\lambda^2 +3\lambda\mu+3 \mu^2)}{\lambda^2},\\
m_d&= \frac{2d-3}{\lambda} \mu^2 m_{d-1}+\mu^2 m_{d-2}.
\end{array}\end{equation}

We have rational identifiability for $d\geq 2$, with

\[\mu=m_1 \quad \text{and}\quad \lambda=\frac{m_1^3}{m_2-m_1^2}.\]

The above system of equations parameterizes the moment variety $\M_d^\mathrm{IG}$ of the inverse Gaussian of degree $d$. Using this parameterization we prove that, for fixed $d$, the moment variety has a determinantal realization and it is Cohen–Macaulay.

Let $d \geq 3$. The homogeneous prime ideal of the inverse Gaussian moment variety $\M^\mathrm{IG}_d$ is generated by   $\binom{d-1}{3}$ cubics and $\binom{d-1}{2}$ quartics, given by the maximal minors of the following $(3\times d)$-matrix

\begin{pmatrix}
m_0^2 & m_0 & m_1 & m_2 & m_3&\cdots& m_{d-2}\\
0 & m_1 & 3m_2 & 5m_3 & 7m_4 &\cdots& (2d-3)m_{d-1}\\
m_1^2 & m_2 & m_3 &m_4  &m_5& \cdots& m_d
\end{pmatrix} .

Furthermore, $\mathcal{M}^\mathrm{IG}_{d}$ is Cohen–Macaulay.

The vector $(\mu^2,\mu^2/\lambda,-1)\neq 0$ is in the left kernel of the matrix $\ker (H_d)$, so the 3-minors of $H_d$ vanish on $\M^\mathrm{IG}_d$. Thus,
$$J_d:=I_3(H_d) \subset \mathcal{I}(\M^\mathrm{IG}_d).$$

This gives that $\dim(\mathcal{V}(J_d))\geq\dim(\M^\mathrm{IG}_d)=3$ as affine varieties. With respect to the reverse lexicographic ordering, we find that

\begin{equation}\label{eq:initial_ideal_inverse_gaussian}
\operatorname{in}(J_d)\supset(m_1^4,m_2^3,\ldots,m_{d-2}^3),
\end{equation}

so $\dim(\mathcal{V}(J_d))=\dim(\mathcal{V}(\init(J_d)))\leq d+1-(d-2)= 3$. Thus, we have proved that
$$\dim (\mathcal{V}(J_d))= \dim (\M^\mathrm{IG}_d)=2$$
as projective varieties. This shows that the ring $S/J_d$, where $S=\CC[m_0,\ldots,m_d]$, has expected Krull dimension $d+1-(d-3+1)=3$, so it is Cohen–Macaulay by <cit.>.
In particular, the ideal $J_d$ has no embedded components.

We want to prove that $J_d$ is prime, so that $J_d= \mathcal{I}(\M^\mathrm{IG}_d)$. This is equivalent to proving that $R=S/J_d$ is a domain. Let $\Delta= m_0^2m_1$ be the submaximal minor in the upper left corner, and consider the localization $R_{\Delta}$.
To prove that $R$ is a domain, it is enough to prove that $\Delta$ is a non-zerodivisor in $R$, and that $R_\Delta$ is a domain.

To prove that $\Delta$ is a non-zerodivisor in $R$, we prove that $m_0$ and $m_1$ are both non-zerodivisors in $R$. Since $R$ is Cohen–Macaulay, we need to prove that $\dim (R/(m_0))= \dim (R/(m_1))= \dim R-1= 2$. Observe that

\[R/(m_0)\cong \dfrac{\CC[m_1,\ldots,m_d]}{I_3 ((H_d)\vert_{m_0=0})},\]

$$(H_d)\vert_{m_0=0}= \begin{pmatrix}
0  & 0 & m_1 & m_2&\cdots& m_{d-2}\\
0 & m_1 & 3m_2 & 5m_3 &\cdots& (2d-3)m_{d-1}\\
m_1^2 & m_2 & m_3 &m_4  & \cdots& m_d
\end{pmatrix}.$$

If we consider the 3-minor on the first three columns, we get $m_1 =0$. Shifting to the next adjacent minor we get $m_2=0$. Thus, iterating, we eventually find that $\mathcal{V}\left(I_3 (H_d) \right)$ is given by the following projective curve:

\begin{align*}
\end{align*}

So, $R/(m_0)$ has affine dimension $2$. A similar argument shows that $R/(m_1)$ is a union of two projective curves:

\begin{align*}
\end{align*}

So, $R/(m_1)$ has affine dimension $2$. This proves that $\Delta$ is a nonzero-divisor in $R$.

It remains to prove that $R_\Delta$ is a domain. Since $\Delta$ is a non-zerodivisor and $R$ is Cohen–Macaulay, we get that $\dim (R_{\Delta})=\dim (R)=3$. In fact, in order for $\dim( R_\Delta)=\dim (R)$, it is sufficient to show that $\Delta$ is not nilpotent. This is because if $\Delta$ is not nilpotent, then it is not in some maximal ideal of $R$, since $R$ is a finitely generated $\CC$-algebra and therefore the nilradical of $R$ is equal to the Jacobson radical of $R$; let $\mathfrak{m}$ denote this maximal ideal of $R$ with $\Delta\notin\mathfrak{m}$. Since $R$ is Cohen–Macaulay, the length of a maximal chain of prime ideals ending at $\mathfrak{m}$ is $\dim(R)=3$, and since $\Delta\notin\mathfrak{m}$, this chain is preserved in $R_\Delta$. Thus, $\dim(R_\Delta)=\dim(R)=3$. We now claim that $R_\Delta$ is a domain. Note that for $3\leq i\leq d$, one can inductively use the $3$-minor of $H_d$ using columns $1,2,i$ so that $\overline{m_i}\in R_\Delta$ can be expressed as

\[\overline{m_i}=\overline{f}/\overline{m_0^2m_1},\qquad \text{for some }f\in\mathbb{C}[m_0,m_1,m_2].\]

Thus, $R_\Delta$ is a finitely generated $\mathbb{C}$-algebra:

\[R_\Delta=\mathbb{C}\left[\overline{m_0},\overline{m_1},\overline{m_2},\overline{1/m_0^2m_1}\right]\subset R.\]

In particular, we get an isomorphism $\mathbb{C}[w,x,y,z]/(w^2xz-1)\cong R_\Delta$, induced by the $\CC$-algebra homomorphism

\[\varphi\from \CC[w,z,y,z]\to R_\Delta,\quad w\mapsto\overline{m_0},\quad x\mapsto\overline{m_1},\quad y\mapsto\overline{m_2},\quad z\mapsto\overline{1/m_0^2m_1}.\]

If $(w^2xz-1)\subsetneq\ker(\varphi)$, then $\dim \mathcal{V}(\ker(\varphi))\leq 2$ since $(w^2xz-1)$ is prime. This is because the height of $\ker(\varphi)$ would need to be strictly greater than the height of $(w^2xz-1)$, and $\CC[w,x,y,z]$ is Cohen–Macaulay so $\dim \mathcal{V}(I)+\text{ht}(I)=4$ for every ideal $I$. However, $\dim(R_\Delta)=3$. Thus it must be that $(w^2xz-1)=\ker(\varphi)$, and therefore $R_\Delta$ is a domain. This completes the proof.

As an immediate consequence of the previous result, we get the following corollaries.

The ideal $\mathcal{I}(\M^\mathrm{IG}_d)$ has componentwise linear (minimal) free resolution given by the Eagon–Northcott complex [22]. In particular, $\reg (\mathcal{I}(\M^\mathrm{IG}_d))=4$

The degree is given by $\deg (\mathcal{I}(\M^\mathrm{IG}_d))= (d-1)^2$.

Since the ideal is Cohen–Macaulay of expected codimension $d-2$, the degree of the variety is the elementary symmetric polynomial of degree $d-2$ in $d$ unknowns, evaluated at $e_1=2,e_2=1,\ldots,e_d=1$ where $e_i$ is the degree of the $i$th column's entries. This yields $\deg (\mathcal{I}(\M^\mathrm{IG}_d))= (d-1)^2.$ The previous formula is known as Thom–Porteous–Giambelli formula. We refer the reader to [15] for further details on the degrees of determinantal varieties.

For computational purposes, it is often useful to have a Gröbner basis of the ideal with respect to a term order. It turns out that the natural generators of $\I(\M^\mathrm{IG}_d)$ are a Gröbner basis with respect to the reverse lexicographic ordering. In order to prove this result, we will need some preliminary results that we collect in the next two lemmas.

Let $S$ be a standard graded polynomial ring and let $I \subset S$ be a homogeneous ideal. Consider the ideal $J=(\init(f_1), \ldots, \init(f_r))$, where $f_1, \ldots, f_r$ are homogeneous elements of $I$. Suppose that the following conditions hold:

* $\dim (S/I)= \dim (S/J)$;

* $J$ is unmixed, i.e., $\height (P)=\height (J)$ for all $P \in \Ass(J)$;

* $\deg(J)=\deg (I)$.

Then $J =\init(I)$.

Let $I$ be a monomial ideal and let $M$ be a minimal generator of $I$ of degree $s$. Write $I=I'+(M)$ for an ideal $I'$. Then, there is a short exact sequence of graded modules and degree 0 maps

\begin{equation}\label{SES}
0 \longrightarrow S/(I':M)[-s] \longrightarrow S/I' \longrightarrow S/I \longrightarrow 0.
\end{equation}

Since the Hilbert series is additive in short exact sequences, we get

\begin{equation}
\HS_{S/I}(t)= \HS_{S/I'}(t)- \HS_{S/(I':M)[-s]}(t).
\end{equation}

The $3\times 3$ minors of $H_d$ form a Gröbner basis for $\I(\M^\mathrm{IG}_d)$ with respect to any antidiagonal term order (for example, the reverse lexicographic ordering).

Let $I=\init(\I(\M^\mathrm{IG}_d))$ and let $J:=\left( \init([i,j,k]) \mid 1\leq i\leq j\leq k \leq d\right)$ where $[i,j,k]$ represents the 3-minor on columns $i,j$ and $k$. Then
$$J= ( m_2,\ldots,m_{d-2})^3+ m_1^2 ( m_1,\ldots,m_{d-2})^2.$$

We want to apply Lemma <ref> to conclude that $J=I$.
Note that $J$ is a primary ideal since every variable appearing in the generators of $J$ has some power in $J$ <cit.>, and $\sqrt{J}=(m_1,\ldots, m_{d-2})$. In particular, $J$ is unmixed and

$$\dim (S/J)=\dim (S/\sqrt{J})=3= \dim (S/I).$$

We are left to show that $\deg(J)=\deg (I)$. We already know that $\deg (I)=(d-1)^2$ by Corollary <ref>.
As for $\deg(J)$, we will apply Lemma <ref> to the ideal $J$ and a monomial generator $M$ of degree 4, i.e., $M \in \{m_1^4, m_1^3 m_2, m_1^2 m_2^2, \ldots, m_1^2 m_{d-2}^2 \}$.

Iteratively applying Lemma <ref>, we finally get:

\begin{equation}\label{HS-SES}
\HS_{S/J}(t)= \HS_{S/(m_2,\ldots,m_{d-2})^3}(t)- \binom{d-1}{2}\HS_{S/(m_2,\ldots,m_{d-2})[-4]}(t)
\end{equation}

where the factor $\binom{d-1}{2}$ appears in the previous identity because we are applying (<ref>) to each of the $\binom{d-1}{2}$ degree 4 monomials arising from a 3-minor of $H_d$ involving the first column. Clearly,

\begin{equation}\label{HS1}
\HS_{S/(m_2,\ldots,m_{d-2})[-4]}(t)=\frac{t^4}{(1-t)^4}.
\end{equation}

As for $S/(m_2, \ldots,m_{d-2})^3$, since it is Cohen–Macaulay we can consider its Artinian reduction
$$\overline{R}= \dfrac{S}{\left((m_2, \ldots,m_{d-2})^3+ (m_0,m_1,m_{d-1},m_d)\right)}.$$
Thus, $\overline{R}$ is a 0-dimensional ring whose Hilbert series is given by the Hilbert polynomial of $\dfrac{\CC[m_2,\ldots,m_{d-2}]}{(m_2,\ldots,m_{d-2})^3}$, that is

\begin{equation*}
\end{equation*}

Since the Hilbert series numerator remains unchanged when taking the Artinian reduction, we get

\begin{equation}\label{HS2}
\end{equation}

Hence, by identity (<ref>),

$$ \HS_{S/J}(t)=\dfrac{1+(d-3)t+\binom{d-1}{2}-\binom{d-1}{2} t^4}{(1-t)^4}=\dfrac{1+(d-2)t+ \binom{d-1}{2}t^2+ \binom{d-1}{2}t^3}{(1-t)^3}.$$

This yields
$$\deg(J)= 1+(d-2)+ 2\binom{d-1}{2}=(d-1)^2=\deg (I).$$

Applying Lemma <ref>, we get that $J=\init(\I(\M^\mathrm{IG}_d))$ and we are done.

We point out that in proving this result, we also found the Hilbert series of the coordinate ring of the moment variety of the inverse Gaussian. In fact, since the Hilbert series does not change when taking the initial ideal, we have that the Hilbert series of $\I(\M^\mathrm{IG}_d)$ is the same as the Hilbert series of $J=\init(\I(\M^\mathrm{IG}_d))$.

The Hilbert series of $S/\I(\M^\mathrm{IG}_d)$ is given by

$$ \dfrac{1+(d-2)t+ \binom{d-1}{2}t^2+ \binom{d-1}{2}t^3}{(1-t)^3}.$$

For $d=3$, we have the following principal homogeneous prime ideal:

\begin{align*}
\mathcal{I}(\M^\mathrm{IG}_3)&=(-m_{0}^{2} m_{1} m_{3}+3 m_{0}^{2} m_{2}^{2}-3 m_{0} m_{1}^{2} m_{2}+m_{1}^{4}
\end{align*}

Its Hilbert series is $\dfrac{1+t+t^2+t^3}{(1-t)^3}$, and its degree is $4$.

For $d=4$, we have the following  homogeneous prime ideal:

\begin{align*}
\mathcal{I}(\M^\mathrm{IG}_4)&=\left(\begin{array}{l}
-m_{0}^{2} m_{1} m_{3}+3 m_{0}^{2} m_{2}^{2}-3 m_{0} m_{1}^{2} m_{2}+m_{1}^{4},
\\
-3 m_{0}^{2} m_{2} m_{4}+5 m_{0}^{2} m_{3}^{2}-5 m_{1}^{3} m_{3}+3 m_{1}^{2} m_{2}^{2},
\\
-m_{0}^{2} m_{1} m_{4}+5 m_{0}^{2} m_{2} m_{3}-5 m_{0} m_{1}^{2} m_{3}+m_{1}^{3} m_{2},
\\
-3 m_{0} m_{2} m_{4}+5 m_{0} m_{3}^{2}+m_{1}^{2} m_{4}-6 m_{1} m_{2} m_{3}+3 m_{2}^{3}
\end{array}\right).
\end{align*}

Its Hilbert series is $\dfrac{1+2t+3t^2+3t^3}{(1-t)^3}$, and its degree is $9$.

We end our discussion about the algebraic properties of $\M^\mathrm{IG}_d$ by computing its singular locus.

The singular locus of $\M^\mathrm{IG}_d$ is given by the projective curve $m_0=m_1=\cdots=m_{d-2}=0$ and the projective point $m_1=m_2= \cdots=m_d=0$.

We begin by noting that the open affine patch $\M^\mathrm{IG}_d\cap\{m_0m_1\neq 0\}$ is included in the smooth locus. To see this, note that the 3-minors involving the first two columns of $H_d$ together give rational expressions for $m_3,m_4,\ldots,m_d$ in $m_0,m_1,m_2$, with denominators that are monomials in $m_0$ and $m_1$, which gives an isomorphism of varieties $(\CC^*)^2\times\CC\cong \M^\mathrm{IG}_d\cap\{m_0m_1\neq 0\}$.
The complement of this affine patch in $\M_d^\mathrm{IG}$ is the union of two projective lines:

$$\mathcal{L}_1:\quad m_0=m_1=\cdots=m_{d-2}=0,$$

$$\mathcal{L}_2:\quad m_1=m_2=\cdots=m_{d-1}=0.$$

The singular locus of $\M^\mathrm{IG}_d$ must therefore be contained in $\mathcal{L}_1\cup\mathcal{L}_2$.

Let $\mathcal{J}$ be the $\binom{d}{3}\times(d+1)$ Jacobian of the maximal minors of $H_d$. The singular locus is precisely the set of points where the rank of $\mathcal{J}$ is less than $\codim(\I(\M^\mathrm{IG}_d))=d-2$.

If we evaluate $\mathcal{J}$ at $\mathcal{L}_1$, we get a matrix with just $d-3$ nonzero entries. To see this, note that the only minors that can give a nonzero contribution to the Jacobian are those that have a term that is at most linear in the variables $m_0,\ldots,m_{d-2}$. The only such maximal minor arises from picking column indices $\{i,d-1,d\}$ for $i\in\{2,\ldots,d-2\}$, and its only contribution to the Jacobian will come from the term $(2d-3)m_{d-1}^2m_{i-2}$.

If, on the other hand, we evaluate $\mathcal{J}$ at $\mathcal{L}_2\setminus\mathcal{L}_1=\mathcal{L}_2\cap\{m_0\neq 0\}$, the only minors of $H_d$ that make a nonzero contribution to the Jacobian, are those that have total degree at most 1 in the variables $m_1,m_2,\ldots,m_{d-1}$. This corresponds to the column indices $\{1,i,d\}$ for $i\not\in\{1,d\}$ (which gives a minor containing the term $(2i-3)m_0^2m_im_d$) or $\{2,i,d\}$ for $i\not\in\{1,2,d\}$ (which gives a minor containing the term $(2i-3)m_0m_im_d$). From this, we see that $\mathcal{J}$ evaluated at $\mathcal{L}_2\cap\{m_0\neq 0\}$ has rank $d-2$ for $m_d\neq 0$, and vanishes completely for $m_d=0$.

Based on Corollary  <ref>, it might be tempting to believe that $\M^\mathrm{IG}_d$ is a Roman surface (a generic projection of a Veronese variety). However, the above proposition shows that this cannot be the case, since the singular locus of a Roman surface consists of 3 points, while the singular locus of $\M^\mathrm{IG}_d$ is given by a curve and a point.

Before closing the subsection, we briefly turn our attention to the cumulants of the inverse Gaussian distribution. The cumulant generating function is

\[K(t)=\log(M(t))=\frac{\lambda}{\mu}\left(1-\sqrt{1- \frac{2\mu^2t}{\lambda}}\right),\]
from which we obtain the following formula for the cumulants:

\[\kappa_r=\frac{(2r-3)!!\mu^{2r-1}}{\lambda^{r-1}}.\]

As for the moment variety, the ideal of the cumulant variety $\K^\mathrm{IG}_d\subset\CC^d$  can be realized as a determinantal ideal.

The ideal $\mathcal{I}(\K^\mathrm{IG}_d)$ is generated by the $\binom{d-1}{2}$ quadrics given by the 2-minors of the matrix

\[K_d=\begin{pmatrix}
-\kappa_1 & \kappa_2 & 3\kappa_3 & \cdots & (2d-3)\kappa_{d-1}\\
\kappa_2 & \kappa_3 & \kappa_4 & \cdots & \kappa_d
\end{pmatrix}.\]

Furthermore, $\mathcal{I}(\K^\mathrm{IG}_d)$ is Cohen–Macaulay, its degree is $d-1$, and its Gröbner basis is given by the $2$-minors of $K_d$.

The recursive relation $\kappa_r=\frac{\mu^2}{\lambda}(2r-3)\kappa_{r-1}$ gives that $(\mu^2/\lambda,-1)\in\ker(K_d)$. We conclude that $J_d:=I_2(K_d)\subset\mathcal{I}(\K^\mathrm{IG}_d)$.
The matrix $H_d$ is a 1-generic Hankel matrix, so $J_d=I_2(K_d)$ is prime with the expected dimension $d-((d-1)-2+1)=2$ [13]. In particular, it is Cohen–Macaulay. Since we have a parameterization $(\CC^*)\times\CC\to\K^\mathrm{IG}_d$ with algebraic identifiability, we also know that $\mathcal{I}(\K^\mathrm{IG}_d)$ is prime of dimension 2. We conclude that $\mathcal{I}(\K^\mathrm{IG}_d)=I_2(K_d)$. The rest follows from known results on Hankel determinantal rings.

Since the ideal obtained above is homogeneous, we note that the cumulant variety is a ruled surface.
In particular, it is a (scaled) rational normal curve when viewed as a projective curve in $\PP^{d-1}$.
An interesting question is whether any features of $\K^\mathrm{IG}_d$ can be pulled back by the birational moments-to-cumulants map $\M_d^\mathrm{IG}\dashrightarrow\K_d^\mathrm{IG}$. For example, we hope that this simple geometric model given by the cumulants could help prove nondefectiveness of secants of moment varieties, which corresponds in statistics to algebraic identifiability of mixtures of univariate inverse gaussian distributions (see <ref> for further details).

§ GAMMA MOMENT VARIETIES

In this section, we study moment varieties $\M_d^\Gamma$ for the gamma distribution. We find their defining ideals and use this to find their degrees and singular loci. We also discuss statistically important special cases of the gamma distribution whose moment varieties are rational normal curves embedded in $\M_d^\Gamma$.

§.§ The gamma distribution

The gamma distribution is commonly used to model physical and economic processes, especially in relation to arrival or waiting times. The density function is supported on $(0,+\infty)$, and involves the gamma function $\Gamma$ (see Figure <ref>). The distribution has two parameters, and there are two commonly used parameterizations:

* The shape-scale parameterization is given by a shape parameter $k>0$ and a scale parameter $\theta>0$. The density function is given by

\[f(x)=\frac{1}{\Gamma(k)\theta^k}x^{k-1}e^{-x/\theta}.\]

* The other parameterization is given by a shape parameter $\alpha>0$ and a rate parameter $\beta>0$. This parameterization is related to the previous one via $\alpha=k$ and $\beta=1/\theta$.

The moment generating function is
$M(t)=(1-\theta t)^{-k}$, and the moments are given by

$$m_i=\prod \limits_{j=0}^{i-1} (k+j)\theta^i= m_{i-1}(k+(i-1))\theta.$$

Let $\M^\Gamma_d$ denote the $d$th moment variety of the gamma distribution. The affine part $\M^\Gamma_d\cap\{m_0\neq0\}$ has previously been studied in <cit.>. Here, we study the full projective variety. We begin by finding its defining ideal. The following proof is similar to the proof of the defining ideal for the inverse Gaussian moment variety (Theorem <ref>), and so we omit some details.

Let $d \geq 3$. The homogeneous prime ideal of the gamma moment variety $\M^\Gamma_d$ is  generated by the $\binom{d}{3}$ cubics given by the maximal minors of the following $(3\times d)$-matrix:

\begin{pmatrix}
0  & m_1 & 2m_2 & 3m_3&\cdots& (d-1)m_{d-1}\\
m_0 & m_1 & m_2 & m_3 &\cdots& m_{d-1}\\
m_1 & m_2 & m_3 &m_4  & \cdots& m_d
\end{pmatrix}.

Furthermore, $\M^\Gamma_d$ is Cohen–Macaulay.

Notice that the vector $(k\theta,\theta, -1)\neq 0$ is in the left kernel of the matrix $\ker (H_d)$, so the 3-minors of $H_d$ vanish on $\M^{\Gamma}_d$. Thus,
$$J_d:=I_2(H_d) \subset \mathcal{I}(\M^{\Gamma}_d).$$

Then, $\mathcal{V}(J_d) \supset \M^{\Gamma}_d$ and $\dim(\mathcal{V}(J_d)) \geq \dim(\M^{\Gamma}_d)=2$ as projective varieties.
On the other hand, if we fix an antidiagonal term order, then

$$ \left(m_i m_j m_k \mid 0<i \leq j <k < d\right) \subset \init(J_d).$$

We observe that $ \sqrt{\left(m_i m_j m_k \mid 0<i \leq j <k < d\right)}= (m_i m_j \mid i \neq j, 0<i,j< d)$, so

$$\dim (\mathcal{V}(J_d)) =\dim (\mathcal{V}(\init(J_d)))\leq  \dim (\mathcal{V} \left(m_i m_j \mid i \neq j, 0<i,j< d\right))=2.$$

Thus, we have proved that $\dim (\mathcal{V}(J_d))= \dim (\M^{\Gamma}_d)=2$ as projective varieties.
This shows that the ring $S/J_d$ has expected Krull dimension $d+1-(d-3+1)=3$, so it is Cohen–Macaulay <cit.>.
In particular, the ideal $J_d$ has no embedded components.

We now show that $J_d$ is prime. Let $\Delta= m_0m_1$ be the submaximal minor in the upper left corner, and consider the localization $R_{\Delta}$. As in the inverse Gaussian case, it is sufficient to prove that $\Delta$ is a non-zero-divisor in $R$, and that $R_\Delta$ is a domain.

To prove that $\Delta$ is a non-zerodivisor in $R$, we prove that $m_0$ and $m_1$ are both non-zero-divisors in $R$. Since $R$ is Cohen–Macaulay, we need to prove that $\dim R/(m_0)= \dim R/(m_1)= \dim R-1= 2$. Observe that

$$R/(m_0)\cong \dfrac{\mathbb{C}[m_1,\ldots,m_d]}{I_3 ((H_d)\vert_{m_0=0})},$$

$$(H_d)\vert_{m_0=0}= \begin{pmatrix}
0  & m_1 & 2m_2 & 3m_3&\cdots& (d-1)m_{d-1}\\
0 & m_1 & m_2 & m_3 &\cdots& m_{d-1}\\
m_1 & m_2 & m_3 &m_4  & \cdots& m_d
\end{pmatrix}.$$

If we consider the 3-minor on the first three columns, we get $m_1^2 m_2=0$. Thus we have two possibilities, either $m_1=0$ or $m_2=0$. If $m_1=0$, then from the next adjacent minor we get $m_2^2 m_3=0$. Iterating this argument, we eventually find that $\mathcal{V}\left(I_3 (H_d) \right)$ is a union of $d-1$ curves given by

\[m_0=m_1=\cdots=\widehat{m_i}=\cdots=m_{d-1}=0,\qquad i=1,\ldots,d-1.\]

Thus, $R/(m_0)$ has affine dimension $2$. A similar argument shows that $R/(m_1)$ has affine dimension $1$. This shows that $\Delta$ is a non-zero-divisor.

The proof that $R_\Delta$ is a domain is similar to the inverse Gaussian case. One can use the $3$-minors of $H_d$ to see that $R_\Delta=\CC[\overline{m_0},\overline{m_1},\overline{m_2},\overline{1/m_0m_1}]\subset R$, and so $R_\Delta\cong\CC[w,x,y,z]/(wxz-1)$. Since $(wxz-1)$ is prime, this shows that $R_\Delta$ is a domain, and this completes the proof.

It is worth noticing that both $\M_d^{\mathrm{IG}}$ and $\M_d^{\Gamma}$ are Cohen–Macaulay. Moreover, to the best of our knowledge, there is no known distribution for which the moment varieties fail to be Cohen–Macaulay. This leads to the following problem:

Classify all Cohen–Macaulay moment varieties from polynomial distributions.

By Theorem <ref>, we have that $\mathcal{I}(M^\Gamma_d)$ is a Cohen–Macaulay ideal generated by the maximal minors of a matrix with linear entries. Therefore, $\mathcal{I}(M^\Gamma_d)$ has a linear minimal free resolution given by the Eagon–Northcott complex; we refer the reader to <cit.> for background on such complexes. From this, one can deduce many invariants of the ideal. In particular, <cit.> gives the degree.

The gamma moment surface $\M^\Gamma_d$ has degree $\binom{d}{2}$ in $\PP^d$.

We now show that the $3$-minors of $H_d$ form a Gröbner basis for $\mathcal{I}(\M^\Gamma_d)$ with respect to the reverse lexicographic order. The outline of the proof is as in the inverse Gaussian case. Let $J$ denote the ideal generated by the initial terms of the minors. We show in the following two lemmas that $J$ is unmixed, and that $\deg(J)=\binom{d}{2}$. We then apply Lemma <ref> to conclude that $J$ is indeed the initial ideal of $\mathcal{I}(\M^\Gamma_d)$.

To show that $J$ is unmixed, we apply the technique of polarization to obtain a squarefree monomial ideal (see [14] for background on polarization). We then use Stanley–Reisner theory to study the associated primes of the polarization.

The ideal $J$ generated by initial terms of the $3$-minors of $H_d$ is unmixed.

Define the ideal $J'$ of $\CC[m_1,\ldots,m_{d-1}]$ by eliminating $m_0$ and $m_d$:

\[J'= J\cap \mathbb{C}[m_1,\ldots,m_{d-1}].\]

If $J'$ is unmixed, then $J$ is as well. We show that every associated prime of $J'$ has height $d-2$.
Since $J'$ is a monomial ideal, its associated primes are also monomial. We therefore show that if $P\in\Ass(J')$, then $P$ is of the form

\[P=(m_{i_1},\ldots,m_{i_{d-2}}\mid 1\leq i_1<\cdots<i_{d-2}\leq d-1).\]

Let $\mathcal{P}(J')$ be the polarization of $J'$, so

\begin{align*}
\mathcal{P}(J')&=(m_{i,1}m_{i,2}m_{j,1},\,m_{i,1}m_{j,1}m_{k,1}\mid 1\leq i<j<k\leq d-1]\\
&\subset S=\mathbb{C}[m_{i,1},m_{i,2},m_{d-1,1}\mid 1\leq i\leq d-2].
\end{align*}

Suppose $Q$ is an associated prime of $\mathcal{P}(J')$, and define the set $A\subset[d-1]$ to be

\[A=\{i\mid m_{i,j}\in Q\}.\]

By the correspondence between associated primes of an ideal and those of its polarization, to show that $J'$ is unmixed, we therefore want to show that $\#A=d-2$ <cit.>.

Let $\Delta$ denote the Stanley–Reisner complex associated to $J'$. Then a prime ideal $Q\subset S$ is an associated prime of $\mathcal{P}(J')$ if and only if the variables of $S$ that are not generators of $Q$ form a facet of $\Delta$.

First, suppose for contradiction that $\#A=d-1$. Let $F$ be the facet corresponding to $Q$. Then $F$ must only involve variables $m_{i,j}$ with $1\leq i\leq d-2$, and it cannot contain both $m_{i,1}$ and $m_{i,2}$ for any $i$. Since $F$ is a facet, the number of $i$ such that $m_{i,1}\in F$ cannot be greater than $2$; otherwise, this would correspond to the generator of $\mathcal{P}(J_d')$ given by $m_{i,1}m_{j,1}m_{k,1}$. We then have the following cases.

* Suppose $m_{i,1},m_{j,1}\in F$ with $i<j$. This implies that

\[F=\{m_{i,1},m_{j,1},m_{k_1,2},\ldots,m_{k_{d-3},2}\mid k_\ell\in[d-2]\setminus\{i,j\}\}.\]

Then, $F$ cannot be a facet, as $F\cup\{m_{j,2}\}$ is a face of $\Delta$.

* If there is at most one element of the form $m_{i,1}$ in $F$, then $F$ cannot be a facet, as it is contained in $F\cup\{m_{d-1,1}\}$, which is also a face of $\Delta$.

Now suppose for contradiction that $\#A<d-2$, and let $F$ be the corresponding facet to $Q$. Then one of the following must be true:

* $m_{i,1},m_{i,2},m_{d-1,1}\in F$ for some $1\leq i\leq d-2$, or

* $m_{i,1}m_{i,2},m_{j,1},m_{j,2}\in F$ for some $1\leq i<j\leq d-2$.

In both cases, we obtain a contradiction that $F$ is a face of $\Delta$, since $\mathcal{P}(J')$ has generators of the form $m_{i,1}m_{i,2}m_{j,1}$ for all $1\leq i<j\leq d-1$.

We therefore see that any associated prime of $J'$ has height $d-2$, and so $J$ is unmixed.

We now prove that $J$ is indeed the initial ideal by computing the degree of $J$; the argument is similar to the proof of Proposition <ref>.

The $3$-minors of $H_d$ form a Gröbner basis for the homogeneous prime ideal of the gamma moment variety with respect to the reverse lexicographic term order.

Let $J$ be the ideal generated by initial terms of $3$-minors of $H_d$. We have that $J\subset\init(\mathcal{I}(\mathcal{M}^\Gamma_d))$, and Lemma <ref> shows that $J$ is unmixed. Therefore, to show equality of the ideals, it suffices by Lemma <ref> to show that the degrees of the two ideals are equal.

We iteratively apply Lemma <ref>. Each time, we remove a minimal generator of $J$ of the form $m_i^2m_j$, with $1\leq i<j\leq d-1$; there are $\binom{d-1}{2}$ such generators. Then the colon ideal is always of the form $(m_{i_1},\ldots,m_{i_{d-3}})$, and the final $J'$ ideal that we end up with is

\[J'=(m_im_jm_k\mid 1\leq i<j<k\leq d-1).\]

Thus, the Hilbert series of $J$ is given by

\begin{align*}
\text{HS}_{S/J}(t)&=\text{HS}_{S/J'}(t)-\binom{d-1}{2}t^3\text{HS}_{S/(m_1,\ldots,m_{d-3})}(t)\\
\end{align*}

Substituting $t=1$ in the numerator, we see that

\[\deg(J)=1+(d-2)+\binom{d-1}{2}=\binom{d}{2}.\]

By Corollary <ref>, this is also the degree of $\mathcal{I}(\mathcal{M}^\Gamma_d)$, and so this completes the proof.

The above proof gives the following corollary, which addresses <cit.>.

The Hilbert series of $S/\mathcal{I}(\mathcal{M}^\Gamma_d)$ is given by

\[\frac{1+(d-2)t+\binom{d-1}{2}t^2}{(1-t)^3}.\]

For $d=3$, we have the following principal homogeneous prime ideal:

\begin{align*}
\I(\M_3^\Gamma)=(-m_0 m_1 m_3 + 2 m_0 m_2^2 - m_1^2 m_2),
\end{align*}

Its Hilbert series is $\HS_{S/\I(\M_3^\Gamma)}(t)=\dfrac{1+t+t^{2}}{\left(1-t\right)^{3}}$, and its degree is $3$.

For $d=4$, we have the following  homogeneous prime ideal:

\begin{align*}
\I(\M_4^\Gamma)=\begin{pmatrix}
-m_0 m_1 m_3 + 2 m_0 m_2^2 - m_1^2 m_2\\ -m_0 m_1 m_4 + 3 m_0 m_2 m_3 - 2 m_1^2 m_3\\ -2 m_0 m_2 m_4 + 3 m_0 m_3^2 - m_1 m_2 m_3\\ -m_1 m_2 m_4 + 2 m_1 m_3^2 - m_2^2 m_3
\end{pmatrix}.
\end{align*}

Its Hilbert series is $\HS_{S/\I(\M_4^\Gamma)}(t)=\dfrac{1+2t+3t^{2}}{(1-t)^3}$, and its degree is $6$.

We now describe the singular locus of $ \M^\Gamma_d$.

The singular locus of $\M^\Gamma_d$ is given by two projective points, defined by the ideals $(m_0,m_1,\ldots,m_{d-1})$ and $(m_1,m_2,\ldots,m_d)$.

Analogously to the proof of Proposition <ref>, one can construct an isomorphism $(\CC^*)^2\times\CC\cong\M^{\Gamma}_d\cap\{m_0m_1\neq 0\}$ by considering the maximal minors involving the first two columns of $H_d$. The complement of this affine open patch in $\M_d^\Gamma$ is the union of $2(d-1)$ projective curves:

* $m_0=m_1=\cdots \widehat{m_{i}} =\cdots = m_{d-2}=m_{d-1}=0$ for $i\in\{1,\ldots,d-1\}$, and

* $m_1=m_2=\cdots \widehat{m_{i}} =\cdots = m_{d-1}=m_{d}=0$ for $i\in\{2,\ldots,d\}$,

where the hat denotes omission. Let $\mathcal{J}$ denote the Jacobian of the maximal minors of $H_d$. If we evaluate $\mathcal{J}$ at one of the curves of the first kind described above, then the only minors that make a contribution to $J$ are those that contain a term that is at most linear in $m_0,\ldots,\widehat{m_i},\ldots,m_{d-1}$. This happens precisely for column indices $\{i,i+1,j\}$ for $j\not\in\{i,i+1\}$, which gives a minor with a term $im_i^2m_{j+1}$, and $\{i+1,j,d\}$ for $j\not\in\{i+1,d\}$, which gives a minor with a term divisivble by $m_im_jm_d$. Hence, the resulting evaluated matrix $\mathcal{J}$ vanishes completely when $m_i=0$, and has rank at least $d-2$ if $m_i\neq 0$. The situation is analogous when $\mathcal{J}$ is evaluated on a curve of the second kind.

We conclude this section by pointing out that, as for the inverse Gaussian, the cumulant variety of the gamma distribution is a ruled surface. This result has been first proved in <cit.>. We include it here for the sake of completeness.

The ideal $\mathcal{I}(\K^\Gamma_d)$ is generated by the $\binom{d-1}{2}$ quadrics given by the 2-minors of the matrix

\[K_d=\begin{pmatrix}
\kappa_1 & \kappa_2 & \dfrac{\kappa_3}{2} & \cdots & \dfrac{\kappa_{d-1}}{(d-2)!}\\
\kappa_2 & \dfrac{\kappa_3}{2} & \dfrac{\kappa_4}{3!} & \cdots & \dfrac{\kappa_d}{(d-1)!}
\end{pmatrix}.\]

Furthermore, $\mathcal{I}(\K^\Gamma_d)$ is Cohen–Macaulay, its degree is $d-1$, and its Gröbner basis is given by the $2$-minors that generate the ideal.

As for the inverse Gaussian, having this simpler geometric model given by the cumulants could help in understanding the geometry of secants of moment varieties of the gamma distribution.

§.§ One-dimensional special cases of the gamma distribution

We now consider two important special cases of the gamma distribution: the exponential and chi-squared distributions. Both distributions are given by one parameter, and have one-dimensional moment varieties.
We show that the moment varieties are in fact rational normal curves.

In [16], the author studies moment varieties for the Poisson distribution, which is also a one-parameter distribution. These varieties were shown to be rational normal curves. It would be interesting to classify projective curves that can appear as moment varieties of one-parameter distributions.

§.§.§ Exponential distribution

The exponential distribution is the distribution of the time between events in a process, where events occur continuously and independently at a constant average rate $\lambda>0$. The only parameter is $\lambda$, and the moment generating function is

\[M(t)=(1-\lambda t)^{-1}.\]

The moments are given by

\[m_i=i! \lambda^i= i \lambda m_{i-1}.\]

This is a specialization of the gamma distribution with shape $k=1$ and rate $\theta=\lambda$ under the shape-scale parameterization.

Let $\M^{\exp}_d$ denote the $d\text{th}$ moment variety of the exponential distribution. The following result shows that the variety is a rational normal curve.

The homogeneous ideal $\I(\M^{\exp}_d)$ is  generated by the maximal minors of

\begin{pmatrix}
m_0 & 2m_1 & 3m_2 & 4m_3&\cdots& dm_{d-1}\\
m_1 & m_2 & m_3 & m_4 &\cdots& m_{d}
\end{pmatrix}.$$

In particular, $\M^{\exp}_d$ is a rational normal curve.

Notice that the vector $(\lambda,-1)\neq 0$ is in $\ker (H_d)$, so the 2-minors of $H_d$ vanish on $\M^{\exp}_d$. Thus,

$$J_d:=I_2(H_d) \subset \mathcal{I}(\M^{\exp}_d).$$

Moreover, $H_d$ is a 1-generic Hankel matrix in the sense of Eisenbud [13], so $J_d=I_2(H_d)$ is prime of expected codimension $d-1$ <cit.>. In particular, it defines a Cohen–Macaulay ring of dimension $2$ <cit.>. As projective varieties,

$$\dim (\mathcal{V}(J_d))= 1= \dim (\M^{\exp}_d).$$

Hence, $\mathcal{V}(J_d)=\M^{\exp}_d$. Since $J_d$ is prime, we get $J_d= \mathcal{I}(\M^{\exp}_d)$.

§.§.§ Chi-squared distribution

The chi-squared distribution is the distribution of a sum of the squares of $k\geq 1$ independent standard Gaussian variables. The only parameter is $k$, and the moment generating function is

\[M(t)=(1-2t)^{-k/2}.\]

The moments are given by

\[m_i=k(k+2)\cdots (k+2i-2)=m_{i-1} (k+2i-2)=m_{i-1} k + m_{i-1} (2i-2).\]

This is a specialization of the gamma distribution with shape $k/2$ and scale $2$.

Let $\mathcal{M}^{\chi}_{d}$ denote the moment variety. Following the proof of Proposition <ref>, we attempt to find generators of the ideal $\I(\mathcal{M}^{\chi}_{d})$ by using the recursive formula given above, which gives that the maximal minors of

\begin{pmatrix}
0  & 2m_1 & 4m_2 & 6m_3&\cdots& (2d-2)m_{d-1}\\
m_0 & m_1 & m_2 & m_3 &\cdots& m_{d-1}\\
m_1 & m_2 & m_3 &m_4  & \cdots& m_d
\end{pmatrix}

vanish on $\mathcal{M}^{\chi}_{d}$.
This matrix is similar to the one for the gamma distribution, and the dimension of the projective variety $\mathcal{V}(I_3(H_d))$ is $(d+1)-(d-3+1)-1=2$. However, the dimension of the projective variety $\mathcal{M}^{\chi}_{d}$ is $1$. Thus, the $3$-minors of $H_d$ are not sufficient to generate the ideal of $\mathcal{M}^{\chi}_{d}$.

Instead, we follow the technique used in <cit.>. Here, the author uses a change in coordinates from moments to powers of the parameter to obtain the defining ideal of the Poisson moment variety. Let $s_{d,i}$ denote the (signed) Stirling numbers of the first kind, and $S_{d,i}$ denote the Stirling numbers of the second kind. Recall that the Stirling numbers of the first kind satisfy

\[x(x-1)(x-2)\cdots(x-n+1)=\sum_{k=0}^n s_{n,k}x^k,\]

and that they are related to the Stirling numbers of the second kind by the following property. Let $(f_n)_{n=0}^\infty$ and $(g_n)_{n=0}^\infty$; then

\[g_d=\sum_{i=0}^d S_{d,i} f_i=\sum_{i=1}^d S_{d,i} f_i \quad \Leftrightarrow \quad f_d=\sum_{i=0}^d s_{d,i} g_i=\sum_{i=1}^d s_{d,i}g_i.\]

Note that $s_{n,0}=S_{n,0}=0$ for all $n\geq 1$.

Using the formula for the $d$th moment, we see that

\[m_d=\sum_{i=0}^d (-2)^{d-i}s_{d,i}k^i,\]

and so applying the relationship between Stirling numbers of the first and second kind with $f_d=m_d/(-2)^d$ and $g_d=(-k/2)^d$, we have that

\[k^d=\sum_{i=0}^d (-2)^{d-i}S_{d,i}m_i.\]

Let $d\geq 2$. The homogeneous prime ideal of the chi-squared moment variety $\mathcal{M}^{\chi}_{d}$ is  generated by the $2\times 2$ minors of the $(2\times d)$-matrix

\[H_d=\begin{pmatrix}
m_0 & m_1 & m_2-2m_1 & \cdots & \sum_{i=0}^{d-1} (-2)^{d-i-1}S_{d-1,i}m_i\\
m_1 & m_2-2m_1 & m_3 - 6m_2 + 4m_1 & \cdots & \sum_{i=0}^d (-2)^{d-i}S_{d,i}m_i
\end{pmatrix}.\]

In particular, $\mathcal{M}^{\chi}_{d}$ is a rational normal curve.

Consider the algebra homomorphism $\varphi\from\CC[m_0,\ldots,m_d]\to\CC[x_0,\ldots,x_d]$ defined by the linear map

\[m_j\mapsto\begin{cases}
\sum_{i=0}^j (-2)^{j-i}s_{j,i}x_i&\quad j\geq 1,\\
x_0&\quad j=0.
\end{cases}\]

This linear map can be given by an upper triangular matrix and is invertible, with inverse

\[x_j\mapsto\begin{cases}
\sum_{i=0}^j (-2)^{j-i}S_{j,i}m_i & \quad j\geq 1,\\
m_0&\quad j=0.
\end{cases}\]

Thus, $\varphi$ is an algebra isomorphism, and induces an isomorphism of varieties $\Phi\from\mathbb{P}^d\to\mathbb{P}^d$, where the coordinates of the domain are given by the $x_i$'s, and the coordinates of the codomain are given by the $m_i$'s.

Under this map, a point $[1:m_1:\cdots:m_d]\in\mathcal{M}^{\chi}_{d}\cap\{m_0=1\}$ is mapped to $[1:x_1:x_1^2:\cdots:x_1^d]$. Thus, the homogenous prime ideal of the projective variety $\Phi(\mathcal{M}^{\chi}_{d})$ is  generated by the $2\times 2$ minors of the matrix

\[H'_d=\begin{pmatrix} x_0 & x_1 & x_2 & \cdots & x_{d-1} \\
x_1 & x_2 & x_3 & \cdots & x_d\end{pmatrix},\]

and we conclude that the the homogenous prime ideal of $\mathcal{M}^{\chi}_{d}$ is given by $\varphi^{-1}(I_2(H'_d))=I_2(H_d)$.

§ FUTURE DIRECTIONS

declare function=pdf_ig(,μ,λ)=sqrt(λ/(2*pi*^3)*exp(-λ*(-μ)^2/(2*μ^2*));
axis lines=left,
column sep=0.2em,
row sep=0.2em,
legend cell align=left,
nodes=scale=0.88, transform shape,
legend entries=
Component with $\mu_1=1$, $\lambda_1=5$,
Component with $\mu_2=2$, $\lambda_2=20$,
Mixture with $\alpha_1=0.4$,
[smooth, domain=0:5, red, dashed]pdf_ig(x,1,5);
[smooth, domain=0:5, blue, dashed]pdf_ig(x,2,20);
[smooth, domain=0:5, purple]0.4*pdf_ig(x,1,5)+0.6*pdf_ig(x,2,20);

Probability density for a mixture of two inverse Gaussians.

It is often the case that real-world datasets exhibit multimodal behavior or heterogeneity. To model these datasets, statisticians use mixtures of distributions which, thanks to their inherent flexibility, allow modelling complex and diverse patterns within data. Mixtures of distributions are simply convex combinations of probability distributions. Geometrically, the moment variety of any probability distribution that is the mixture of $k$ independent observations of the same distribution corresponds to the $k$th secant variety of $\M_d$, and it is denoted by $\Sec_k(\M_d)$.

§.§ Algebraic identifiability of mixtures

Let $\M_d$ be either  $\M^\mathrm{IG}_d$ or $\M^{\Gamma}_d$. We are interested in studying the $k$th secant variety of  $\M_d$. This variety is defined as the projective closure in $\PP^d$ of the image of the parametrization

\begin{equation}\label{eq:parametrization_secant_variety}
\begin{array}{l}\Psi\from \CC^k\times\CC^k\times\CC^{k-1}\to\CC^d\\\hspace{2em}(\lambda,\mu,\alpha)\mapsto \left(\textstyle\sum_{i=1}^{k-1}\alpha_im_r(\lambda_i,\mu_i)+(1-\textstyle\sum_{i=1}^{k-1}\alpha_i)m_r(\lambda_k,\mu_k)\right)_{r=1,\ldots,d}\,.\end{array}
\end{equation}

Clearly, the dimension of  $\Sec_k(\M_d)$ is  bounded above by the number of parameters, so we have

$$ \dim \left(\Sec_k (\M_d)\right) \leq \min\{d, 3k-1\}.$$

If equality holds we say that $\Sec_k(\M_d)$ is nondefective; otherwise, we say that $\M_d$ is $k$-defective. In particular, if $\Sec_k(\M_d)$ is nondefective and $d\geq 3k-1$, then the parameters of each component, as well as the mixing coefficients, are identifiable from the $3k-1$ first moments. This motivates the study of nondefectiveness of $\Sec_k(\M_d)$.

The dimension of $\Sec_k(\M_d)$ is bounded below by the rank of the Jacobian of the parametrization (<ref>). By computing this rank for randomly chosen points with rational entries in exact arithmetic in , we are able to computationally verify nondefectiveness for all $d,k\leqslant 100$.
Thus, it seems natural to conjecture the following.

$\Sec_k (\M^\mathrm{IG}_d)$ and $\Sec_k (\M^{\Gamma}_d)$ are always nondefective.

A helpful step towards a proof of this conjecture is
Terracini's classification of $k$-defective surfaces, which establishes that the only $k$-defective surfaces are either cones or quadratic Veronese embeddings of a rational normal surface in $\PP^k$ [10]. Following [5] and using the above classification, we get the following.

Let $\M_d$ be either $\M^\mathrm{IG}_d$ or $\M^{\Gamma}_d$. If $\M_d$ is $k$-defective, then it is contained in a cone over a curve, with apex a linear space of dimension $\leq k-2$.

Assume for a contradiction that $\M_d^\mathrm{IG}$ (resp. $\M_d^\Gamma$) is not contained in such a cone. Then, by Terracini's lemma, it should be the quadratic Veronese embedding of a rational normal surface in $\PP^k$. But this is a contradiction, since the singular locus of $\M_d^\mathrm{IG}$ (resp. $\M_d^\Gamma$) contains a line by Proposition <ref> (resp. Proposition <ref>).

However, ruling out the possibility of $\M_d$ being contained within a cone is not straightforward and requires further investigation. This is an ongoing work with Kristian Ranestad.

§.§ Identifiability degree

Suppose $\dim\left(\Sec_k(\M_d)\right)=\min\{3k-1,d\}$ for both the inverse Gaussian and gamma distribution. Then, for $d=3k-1$, we have generically finite fibers of the parameterization (<ref>).
The generic cardinality of these fibers is the identifiability degree, which measures the complexity of solving the system (<ref>). Identifiability degrees for mixtures of Gaussians have previously been studied with numerical algebraic geometry tools [4], and here we take a similar approach.

Table <ref> contains bounds for the identifiability degree up to the label-swapping symmetry for $d=3k-1$.
The lower bounds were computed using monodromy and certification in . Note that the certification step constitutes a proof that the lower bounds are true.
The three values for the Gaussian distribution were computed in the works [23], [2] and [4], respectively.

$k$     2     3      4

9     225     $\geq10\,350$

Gamma     $\geq 9$      $\geq 271$     $\geq 288\,091$

Inverse Gaussian     $\geq 24$      $\geq 1637$     $> 1\,200\,000$

Identifiability degrees for $k$-mixtures of three distributions, up to label-swapping symmetry. The monodromy calculation for $k=4$ for the inverse Gaussian was manually terminated, and is therefore not expected to be sharp.

Many open questions in this direction remain. Apart from computing the very large identifiability degree for $k=4$ in the inverse Gaussian case, one can also try to fix some of the parameters or impose equality of some of the parameters of the components of the mixtures to get more manageable systems, for which it is easy to find explicit formulas. This has been a successful line of research in the Gaussian case [20, 1].

Another possible direction is to better understand the smaller identifiability degrees already found in Table <ref> from the point of view of [23].

Find explicit univariate polynomials of degree 9 and 24, analogous to the “Person polynomial” found in [23], that account for the identifiability degrees for $k=2$ in Table <ref>, e.g., using a similar computer-algebra strategy as in <cit.>.

Finally, we also note that numerical experiments for $k\leq 4$ indicate that a single additional polynomial is enough to decrease the generic number of solutions given in Table <ref> to 1. This would be in line with <cit.> in the Gaussian case, and we therefore propose the following.

For $k$-mixtures of the inverse Gaussian distribution or the gamma distribution, we have rational identifiability from the first $3k$ moments.

A natural first step could be to prove the more modest claim of rational identifiability from $3k+2$ moments, for instance by attempting to adapt the techniques previously employed to prove this in the Gaussian setting <cit.>.

§.§ Euclidean distance degree

For noisy sample moments, the system (<ref>) might not have any statistically meaningful solutions. In that case, a common approach is to instead solve the optimization problem

\begin{equation}\label{eq:optimization_problem}
\min_{\theta\in\CC^n} \Vert m(\theta)-\widehat{m}\Vert^2,
\end{equation}

where $m(\theta)=(m_1(\theta),\ldots,m_d(\theta))$ is the first $d$ moments as funcitons of $\theta$, $\widehat{m}\in\RR^d$ is the first $d$ sample moments, and $d>n$ is such that $\dim(\M_d)=n$. Here, $\Vert{\cdot}\Vert$ denotes the Euclidean norm, but we remark that the computations done in this section can be adapted also to weighted norms $\Vert{\cdot}\Vert_W$ with $\Vert x\Vert_W^2=x^\top Wx$ for some (positive semidefinite) matrix $W\in\RR^{d\times d}$. Such norms are often considered in the more general framework of the generalized method of moments [18].

Compared to maximum likelihood estimation (MLE), which instead minimizes the Kullback–Liebler divergence, the distance-based approach studied here has the advantage of being a polynomial optimization problem. Hence, it has a finite optimization degree, which is the number of complex critical points for generic sample moments $\widehat{m}_i$. Understanding this degree is an open problem, which we in this section take some first steps towards exploring, from the lense of the related notion of the Euclidean distance degree (ED degree).

The ED degree is the generic number of critical points of the optimization problem

\begin{equation}\label{eq:ED_problem}
\min_{m\in\M_d\cap\{m_0=1\}}\Vert m-\widehat{m}\Vert^2.
\end{equation}

This notion was introduced in [12], and has been shown to have a rich geometric meaning in terms of, e.g., Euler characteristics [21], which might make it more approachable to compute than the optimization degree of (<ref>). Since the optimization degree is bounded below by the product of the indentifiability degree and the ED degree of $\M_d$, understanding the ED degree is a natural subproblem.

In this work, we take a first step in this direction by computing the ED degree of $\M_3$ for the inverse Gaussian and gamma distribution. For comparison purposes, we also do this for the Gaussian distribution, which to our knowledge has not been studied from the ED degree perspective before. In all three cases, the affine patch $\{m_0=1\}$ is a hypersurface of the form $\mathcal{V}(f)\subset\CC^3$, where $f\in\CC[m_1,m_2,m_3]$ is given in Examples <ref> and <ref>, and <cit.>, respectively. The ED degree is the number of complex roots of the system

\begin{equation}
\label{eq:ED_system_hypersurface}
f(m)=\nabla f(m)-a(m-\widehat{m})=ab-1=0,\quad m\in\CC^3,\quad a,b\in\CC
\end{equation}

for generic parameters $\widehat{m}\in\CC^3$. A standard Gröbner basis calculation over the field $\CC(\widehat{m}_1,\widehat{m}_2,\widehat{m}_3)$ of rational functions in the sample moments proves the following result.

For the inverse Gaussian distribution, the gamma distribution, and the Gaussian distribution, the ED degree is given by

\[\operatorname{EDdegree}(\M_3^\Gamma)=12\,,\quad \operatorname{EDdegree}(\M_3^\mathrm{IG})=10,\quad \operatorname{EDdegree}(\M_3^\mathrm{Gaussian})=7. \]

This means that in all three cases, $f$ is generic enough to satisfy the mixed volume bounds given in  <cit.>, but not generic enough to satisfy the bound of <cit.>, which evaluates to 52, 21, and 21, respectively. The determinantal realizations of the three hypersurfaces are also not generic enough in the family of $3\times 3$ matrices to give the value $1$ obtained from <cit.>.

A natural future direction is to compute and analyze the ED degrees for higher values of $d$, e.g., using numerical algebraic geometry and ideas from the emerging field of metric algebraic geometry. It would also be interesting to understand the ED degree of $\Sec_k(\M_d)$; since the implicitization problem for these secant varieties turns out to be very hard in our experience, this is a considerable challenge, even with numerical techniques.

[1]
Daniele Agostini, Carlos Améndola, and Kristian Ranestad.
Moment identifiability of homoscedastic gaussian mixtures.
Foundations of Computational Mathematics, 21:695–724, 2021.

[2]
Carlos Améndola, Jean-Charles Faugère, and Bernd Sturmfels.
Moment varieties of Gaussian mixtures.
J. Algebr. Stat., 7(1):14–28, 2016.

[3]
Yulia Alexandr, Joe Kileel, and Bernd Sturmfels.
Moment varieties for mixtures of products.
In Proceedings of the International Symposium on Symbolic
& Algebraic Computation (ISSAC 2023), pages 53–60. ACM, New York,

[4]
Carlos Améndola, Julia Lindberg, and Jose Israel Rodriguez.
Solving parameterized polynomial systems with decomposable
projections, 2016.

[5]
Carlos Améndola, Kristian Ranestad, and Bernd Sturmfels.
Algebraic identifiability of Gaussian mixtures.
Int. Math. Res. Not. IMRN, (21):6556–6580, 2018.

[6]
Wilfred Bruns, Aldo Conca, Claudiu Raicu, and Matteo Varbaro.
Determinants, Gröbner Bases and Cohomology.
Springer Monographs in Mathematics, 2022.

[7]
Mikhail Belkin and Kaushik Sinha.
Polynomial learning of distribution families.
SIAM Journal on Computing, 44(4):889–911, 2015.

[8]
Paul Breiding, Frank Sottile, and James Woodcock.
Euclidean distance degree and mixed volume.
Foundations of Computational Mathematics, 22(6):1743–1765,

[9]
Winfried Bruns and Udo Vetter.
Determinantal rings, volume 1327 of Lecture Notes in
Springer-Verlag, Berlin, 1988.

[10]
Luca Chiantini and Ciro Ciliberto.
Weakly defective varieties.
Trans. Amer. Math. Soc., 354(1):151–178, 2002.

[11]
Ciro Ciliberto, Maria Angelica Cueto, Massimiliano Mella, Kristian Ranestad,
and Piotr Zwiernik.
Cremona linearizations of some classical varieties.
In Gianfranco Casnati, Alberto Conte, Letterio Gatto, Livia Giacardi,
Marina Marchisio, and Alessandro Verra, editors, From Classical to
Modern Algebraic Geometry, pages 375–407, Cham, 2016. Springer
International Publishing.

[12]
Jan Draisma, Emil Horobeţ, Giorgio Ottaviani, Bernd Sturmfels, and Rekha R.
The Euclidean distance degree of an algebraic variety.
Found. Comput. Math., 16(1):99–149, 2016.

[13]
David Eisenbud.
On the resiliency of determinantal ideals.
In Commutative algebra and combinatorics (Kyoto, 1985),
volume 11 of Adv. Stud. Pure Math., pages 29–38. North-Holland,
Amsterdam, 1987.

[14]
Sara Faridi.
Monomial ideals via square-free monomial ideals.
In Commutative algebra, volume 244 of Lect. Notes Pure
Appl. Math., pages 85–114. Chapman & Hall/CRC, Boca Raton, FL, 2006.

[15]
William Fulton and Piotr Pragacz.
Schubert Varieties and Degeneracy Loci.
Springer Berlin, Heidelberg, 2006.

[16]
Alexandros Grosdos Koutsoumpelias.
Algebraic Methods for the Estimation of Statistical
PhD thesis, Osnabrück University, December 2020.

[17]
Alexandros Grosdos Koutsoumpelias and Markus Wageringel.
Moment ideals of local Dirac mixtures.
SIAM J. Appl. Algebra Geom., 4(1):1–27, 2020.

[18]
Lars Peter Hansen.
Large sample properties of generalized method of moments estimators.
Econometrica, 50(4):1029–1054, 1982.

[19]
Kathlén Kohn, Boris Shapiro, and Bernd Sturmfels.
Moment varieties of measures on polytopes.
Ann. Sc. Norm. Super. Pisa Cl. Sci. (5), 21:739–770, 2020.

[20]
Julia Lindberg, Carlos Améndola, and Jose Israel Rodriguez.
Estimating Gaussian mixtures using sparse polynomial moment
systems, 2021.

[21]
Laurentiu G. Maxim, Jose I. Rodriguez, and Botong Wang.
Euclidean distance degree of the multiview variety.
SIAM Journal on Applied Algebra and Geometry, 4(1):28–48,

[22]
Uwe Nagel and Tim Römer.
Criteria for componentwise linearity.
Communications in Algebra, 43(3):935–952, 2015.

[23]
Karl Pearson.
Contributions to the mathematical theory of evolution.
Philosophical Transactions of the Royal Society of London. A,
185:71–110, 1894.

Authors' addresses:

Oskar Henriksson, University of Copenhagen<EMAIL_ADDRESS>
Lisa Seccia, MPI-MiS Leipzig<EMAIL_ADDRESS>
Teresa Yu, University of Michigan<EMAIL_ADDRESS>