# How Much is Unseen Depends Chiefly on Information About the Seen

Seongmin Lee  Marcel Böhme

###### Abstract

It might seem counter-intuitive at first: We find that, _in expectation_ , the
proportion of data points in an unknown population—that belong to classes that
do _not_ appear in the training data—is almost entirely determined by the
number $f_{k}$ of classes that _do_ appear in the training data the same
number of times. While in theory we show that the difference of the induced
estimator decays exponentially in the size of the sample, in practice the high
variance prevents us from using it directly for an estimator of the sample
coverage. However, our precise characterization of the dependency between
$f_{k}$’s induces a large search space of different _representations_ of the
expected value, which can be deterministically instantiated as estimators.
Hence, we turn to optimization and develop a genetic algorithm that, given
only the sample, searches for an estimator with minimal mean-squared error
(MSE). In our experiments, our genetic algorithm discovers estimators that
have a substantially smaller MSE than the state-of-the-art Good-Turing
estimator. This holds for over 96% of runs when there are at least as many
samples as classes. Our estimators’ MSE is roughly 80% of the Good-Turing
estimator’s.

Good-Turing frequency estimation, Multinomial probability estimation, Unseen
events, Missing mass, Probability mass

##  1 Introduction

Suppose, we are drawing balls with replacement from an urn with colored balls.
What is the proportion of balls in that urn that have a color not observed in
the sample; or equivalently, what is the probability that the next ball has a
previously unobserved color? What is the distribution of rarely observed
colors in that urn? These questions represent a fundamental problem in machine
learning: How can we extrapolate from properties of the training data to
properties of the unseen, underlying distribution of the data?

###  1.1 Background

Consider a _multinomial distribution_ $p=\langle p_{1},\cdots p_{S}\rangle$
over a support set $\mathcal{X}$ where support size $S=|\mathcal{X}|$ and
probability values are _unknown_. Let $X^{n}=\langle X_{1},\cdots
X_{n}\rangle$ be a set of independent and identically distributed random
variables representing the sequence of elements observed in $n$ samples from
$p$. Let $N_{x}$ be the number of times element $x\in\mathcal{X}$ is observed
in the sample $X^{n}$. For $k:0\leq k\leq n$, let $\Phi_{k}$ be the number of
elements appearing exactly $k$ times in $X^{n}$.

$\displaystyle N_{x}=\sum_{i=1}^{n}\mathbbm{1}(X_{i}=x)\text{\quad and
\quad}\Phi_{k}=\sum_{x\in\mathcal{X}}\mathbbm{1}(N_{x}=k)$

Let $f_{k}(n)$ be the expected value of $\Phi_{k}$ (Good, 1953), i.e.,

$\displaystyle f_{k}(n)={n\choose
k}\sum_{x\in\mathcal{X}}p_{x}^{k}(1-p_{x})^{n-k}=\mathbb{E}\left[\Phi_{k}\right]$

Estimating the multinomial. Suppose, we want to estimate $p$. We cannot expect
all elements to exist in $X^{n}$. While the empirical estimator
$\hat{p}_{x}^{\text{Emp}}=N_{x}/n$ is generally unbiased,
$\hat{p}_{x}^{\text{Emp}}$ distributes the entire probability mass only over
the observed elements. This leaves a “missing probability mass” over the
unobserved elements. In particular, $\hat{p}_{x}^{\text{Emp}}$ _given that_
$N_{x}>0$ overestimates $p_{x}$, i.e., for observed elements

$\displaystyle\mathbb{E}\left[\left.\frac{N_{x}}{n}\ \right|\ N_{x}>0\right]$
$\displaystyle=\frac{p_{x}}{1-(1-p_{x})^{n}}.$

We notice that the bias increases as $p_{x}$ decreases. Bias is maximized for
the rarest observed element.

Missing mass, realizability, and natural estimation. Good and Turing (GT)
(Good, 1953) discovered that the expected value of the probability
$M_{k}=\sum_{x\in\mathcal{X}}p_{x}\mathbbm{1}(N_{x}=k)$ that the $(n+1)$-th
observation $X_{n+1}$ is an element that has been observed exactly $k$ times
in $X^{n}$ (incl. $k=0$) is a function of the expected number of colors
$f_{k+1}(n+1)$ that will be observed $k+1$ times in an enlarged sample
$X^{n}\cup X_{n+1}$.

$\mathbb{E}\left[M_{k}\right]=\frac{k+1}{n+1}f_{k+1}(n+1).$ (1)

We also call $M_{k}$ as _total probability mass_ over the elements that have
been observed exactly $k$ times. Since our sample $X^{n}$ is only of size $n$,
GT suggested to estimate $f_{k+1}(n+1)$ using $\Phi_{k+1}$. Concretely,
$\hat{M}_{k}^{G}=\frac{k+1}{n}\Phi_{k+1}$.

For $k=0$, $M_{0}$ gives the _“missing” (probability) mass_ over the elements
not in the sample. In genetics and biostatistics, the complement $1-M_{0}$
measures _sample coverage_ , i.e., the proportion of individuals in the
population belonging to a species _not_ observed in the sample (Chao & Jost,
2012). In the context of supervised machine learning, assuming the training
data is a random sample, the sample coverage of the training data gives the
proportion of all data (seen or unseen) with labels not observed in the
training data.

Using an estimate $\hat{M}_{k}$ of $M_{k}$, we estimate the probability
$p_{x}$ of an element $x\in\mathcal{X}$ that appears $k$ times as
$\hat{M}_{k}/\Phi_{k}$. The estimation of $p$ in the presence of unseen
elements $x\not\in X^{n}$ is a _fundamental problem in machine learning_
(Orlitsky et al., 2003; Orlitsky & Suresh, 2015; Painsky, 2022; Acharya et
al., 2013; Hao & Li, 2020). For instance, in natural language processing the
estimation of the probability of a given sequence of words occurring in a
sentence is the main challenge of language models, particularly in the
presence of sequences that appear in the training data rarely or not at all.
Different smoothing and backoff techniques have been developed to tackle this
data sparsity challenge.

A _natural estimator_ of $p_{x}$ assigns the same probability to all elements
$x$ appearing the same number of times in the sample $X^{n}$ (Orlitsky &
Suresh, 2015). For $k>0$, $\hat{p}_{x}=M_{k}/\Phi_{k}$ gives the
hypothetical111The best natural estimator is also called oracle-aided
estimator for its knowledge about $p_{x}$ (Orlitsky & Suresh, 2015) but cannot
actually be used for estimation. _best natural estimator_ of $p_{x}$ for every
element $x$ that has been observed $k$ times.

Bias of Good-Turing (GT). In terms of bias, Juan and Lo (Juang & Lo, 1994)
observe that the GT estimator $\hat{M}_{k}^{G}=\frac{k+1}{n}\Phi_{k+1}$ is an
unbiased estimate of $M_{k}(X^{n-1})$, i.e., where the $n$-th sample was
_deleted_ from $X^{n}$ and find:

$\displaystyle\left|\mathbb{E}\left[\hat{M}_{k}^{G}-M_{k}\right]\right|$
$\displaystyle=\left|\mathbb{E}\left[M_{k}(X^{n-1})-M_{k}(X^{n})\right]\right|$
$\displaystyle\leq\frac{k+2}{n+1}=\mathcal{O}\left(\frac{1}{n}\right).$

Convergence of GT. McAllester and Schapire (McAllester & Schapire, 2000)
analyzed the _convergence_. With high probability,

$\displaystyle|\hat{M}_{k}^{G}-M_{k}(n)|$
$\displaystyle=_{\delta}\begin{cases}\mathcal{O}\left(1/\sqrt{n}\right)&\text{for
$k=0$ (missing mass)}\\\ \mathcal{O}\left(\log(n)/\sqrt{n}\right)&\text{for
small $k$ (rare elements)}\\\ \mathcal{O}\left(k/\sqrt{n}\right)&\text{for
large $k$ (abundant elements).}\end{cases}$

This result was improved by Drukh and Mansour (Drukh & Mansour, 2004) and more
recently by Painsky (Painsky, 2022) who showed that GT estimator converges at
a rate of $\mathcal{O}(1/\sqrt{n})$ for all $k$ based on worst-case mean
squared error analysis.

Competitiveness of GT. Using the Poisson approximation, Orlitsky and Suresh
(Orlitsky & Suresh, 2015) showed that natural estimators from GT, i.e.,
$\hat{p}_{x}^{G}=\hat{M}_{N_{x}}^{G}/\Phi_{N_{x}}$, performs close to the best
natural estimator. Regret, the metric of the competitiveness of an estimator
against the best natural estimator, is measured as KL divergence between the
estimate $\hat{p}$ and the actual distribution $p$, $D_{KL}(\hat{p}||p)$.
Orlitsky and Suresh also showed that finding the best natural estimator for
$p$ is same as finding the best estimator for $M=\\{M_{k}\\}_{k=0}^{n}$.

Poisson approximation. The Poisson approximation with parameter
$\lambda_{x}=p_{x}n$ has often been used to tackle a major challenge in the
formal analysis of the missing mass and natural estimators (Orlitsky & Suresh,
2015; Orlitsky et al., 2016; Acharya et al., 2013; Efron & Thisted, 1976;
Valiant & Valiant, 2016; Good, 1953; Good & Toulmin, 1956; Hao & Li, 2020).
The challenge is the _dependencies between frequencies_ $N_{x}$ for different
elements $x\in\mathcal{X}$. In this Poisson Product model, a continuous-time
sampling scheme with $S=|\mathcal{X}|$ independent Poisson distributions is
considered where the frequency $N_{x}$ of an element $x$ is represented as a
Poisson random variable with mean $p_{x}n$. In other words, in the Poisson
approximation, the frequencies $N_{x}$ are modelled as independent random
variables. Consequently, GT estimator is an unbiased estimator for the Poisson
Product model (Orlitsky et al., 2016), yet it is biased in the multinomial
distribution (Juang & Lo, 1994). Hence, we tackle the dependencies between
frequencies analytically, without approximation via the Poisson Product model.

###  1.2 Contribution of the Paper

In this paper, we reinforce the foundations of multinomial distribution
estimation with a precise characterization of the the _dependencies_ between
$N_{x}=\sum_{i=1}^{n}\mathbbm{1}(X_{i}=x)$ across different $x\in\mathcal{X}$
(rather than assuming independence and using the Poisson approximation). The
theoretical analysis is based on the _expected value_ of the frequency of
frequencies $\mathbb{E}\left[\Phi_{k}\right]=f_{k}(n)$ between different $k$
and $n$, which is

$\frac{f_{k}(n)}{{n\choose k}}=\frac{f_{k}(n+1)}{{n+1\choose
k}}-\frac{f_{k+1}(n+1)}{{n+1\choose k+1}}.$ (2)

Exploring this new theoretical tool, we bring two contributions to the
estimation of the total probability mass $M_{k}$ for any $k:0\leq k\leq n$.
Firstly, we show _exactly to what extent_ $\mathbb{E}\left[M_{k}\right]$ can
be estimated from the sample $X^{n}$ and _how much remains_ to be estimated
from the underlying distribution $p$ and the number of elements
$|\mathcal{X}|$. Specifically, we show the following.

######  Theorem 1.1.

$\displaystyle\mathbb{E}\left[M_{k}\right]={n\choose
k}\left[\sum_{i=1}^{n-k}(-1)^{i-1}f_{k+i}(n)\left/{n\choose
k+i}\right.\right]+R_{n,k}$

where $R_{n,k}={n\choose k}(-1)^{n-k}f_{n+1}(n+1)$ is the remainder.

This decomposition shows that the GT estimator is the _first term_ of
$\mathbb{E}\left[M_{k}\right]$ using the plug-in estimator $\Phi_{1}$ for
$f_{1}(n)$. Hence, it gives the _exact bias_ of the GT estimator in the
multinomial setting (which would incorrectly be identified as _unbiased_ using
the Poisson approximation).

Our second contribution is the development of a class of natural estimators.
We start by defining a nearly unbiased estimator $\hat{M}_{k}^{B}={n\choose
k}\left[\sum_{i=1}^{n-k}(-1)^{i-1}\Phi_{k+i}\left/{n\choose
k+i}\right.\right]$ that uses the plug-in estimator $\Phi_{l}$ for $f_{l}(n)$
in Theorem 1.1 and drops $R_{n,k}$. While the bias of $\hat{M}_{k}^{B}$ decays
exponentially, the variance of $\hat{M}_{k}^{B}$ is too high to be practical.

We observe that we can cast the estimation of the expected total mass as an
optimization problem. From Theorem 1.1 and Eqn. 2, we can see that there are
many _representations_ of $\mathbb{E}\left[M_{k}\right]$, all of which suggest
different estimators for $\mathbb{E}\left[M_{k}\right]$. We introduce a
deterministic method to construct a unique estimator from a representation,
and show how to estimate the mean squared error (MSE) for such an estimator.
Equipped with a large _search space_ of representations and a _fitness
function_ to estimate the MSE of a candidate estimator, we develop a genetic
algorithm that finds an estimator with a minimal MSE (and at most N terms).

We compare the performance of the minimal-bias estimator $\hat{M}_{k}^{B}$ and
the minimal-MSE estimators discovered by our genetic algorith to the that of
the widely used GT estimator on a variety of multinomial distributions used
for evaluation in previous work. Our results show that 1) the minimal-bias
estimator has a substantially smaller bias than the GT estimator by thousands
of order of magnitude, 2) Our genetic algorithm can produce estimators with
MSE smaller than the GT estimator over 96% of the time when there are at least
as many samples as classes; their MSE is roughly 80% of the GT estimator. We
also publish all data and scripts to reproduce our results.

##  2 Dependencies Between Frequencies $N_{x}$

Firstly, we propose a new, distribution-free222A _distribution-free analysis_
is free of assumptions about the shape of the probability distribution
generating the sample. In this case, we make no assumptions about parameters
$p$ or $n$. methodology for reasoning about properties of estimators of the
missing and total probability masses for multinomial distributions. The _main
challenge_ for the statistical analysis of $M_{k}$ has been reasoning in the
presence of dependencies between frequencies $N_{x}$ for different elements
$x\in\mathcal{X}$. As discussed in Section 1.1, a Poisson approximation with
parameter $\lambda_{x}=p_{x}n$ is often used to render these frequencies as
independent (Orlitsky & Suresh, 2015; Orlitsky et al., 2016; Acharya et al.,
2013; Efron & Thisted, 1976; Valiant & Valiant, 2016; Good, 1953; Good &
Toulmin, 1956; Hao & Li, 2020). In the following, we tackle this challenge by
formalizing these dependencies between frequencies. Thus, we establish a link
between the expected values of the corresponding total probability masses.

###  2.1 Dependency Among Frequencies

Recall that the expected value $f_{k}(n)$ of the number of elements $\Phi_{k}$
that appear exactly $k$ times in the sample $X^{n}$ is defined as
$f_{k}(n)=\sum_{x\in\mathcal{X}}\binom{n}{k}p_{x}^{k}(1-p_{x})^{n-k}$. For
convencience, let $g_{k}(n)=f_{k}(n)/{n\choose k}$. We notice the following
relationship among $k$ and $n$:

$\displaystyle g_{k}(n+1)$
$\displaystyle=\sum_{x=1}^{S}p_{x}^{k}(1-p_{x})^{n-k}\cdot(1-p_{x})$
$\displaystyle=g_{k}(n)-g_{k+1}(n+1)$ (3)
$\displaystyle=\sum_{i=0}^{n-k}(-1)^{i}g_{k+i}(n)+(-1)^{n-k+1}g_{n+1}(n+1)$

We can now write the expected value $\mathbb{E}\left[M_{k}\right]$ of the
total probability mass in terms of the frequencies with which different
elements $x\in\mathcal{X}$ have been observed in the sample $X^{n}$ of size
$n$ as follows

$\displaystyle\mathbb{E}\left[M_{k}\right]$
$\displaystyle=\sum_{x\in\mathcal{X}}\binom{n}{k}p_{x}^{k+1}(1-p_{x})^{n-k}$
$\displaystyle={n\choose k}g_{k+1}(n+1)$ (4) $\displaystyle={n\choose
k}\left[\sum_{i=1}^{n-k}(-1)^{i-1}g_{k+i}(n)\right]+R_{n,k}$ (5)

where $R_{n,k}={n\choose k}(-1)^{n-k}f_{n+1}(n+1)$ is a remainder term. This
demonstrates Theorem 1.1.

Figure 1: Lower triangle matrix of $g_{k}(n)$.

Figure 1 illustrates the relationship between the expected frequency of
frequencies $f_{k}(n)=g_{k}(n)/{n\choose k}$, the frequency $k$, and the
sample size $n$. The y- and x-axis represents the sample size $n$ and the
frequency $k$, respectively. As per Eqn. (3), for every 2x2 lower triangle
matrix, the value of the lower left cell ($g_{k}(n+1)$) is value of the upper
left cell ($g_{k}(n)$) minus the value of the lower right cell
($g_{k+1}(n+1)$).

We can use this visualization to quickly see how to rewrite $g_{k}(n)$ as an
alternating sum of values of the cells in the upper row, starting from the
cell in the same column to the rightmost cell, and adding/subtracting the
value of the rightmost cell in the current row. For instance, the value
$g_{0}(13)$ in the bottom red cell (row $n=13$, column $k=0$) in the figure
can be equally calculated as the various linear combinations of its
surrounding cells: (1) with $g_{0}(12)$ and $g_{1}(13)$ (blue colored), (2)
with $g_{0}(11)$, $g_{1}(11)$, $\cdots$, $g_{4}(13)$ (purple colored), or (3)
with $g_{0}(11)$, $g_{1}(11)$, $\cdots$, $g_{8}(13)$ (green colored).

Missing Mass. The missing probability mass $M_{0}$ gives the proportion of
_all possible observations_ for which the elements $x\in\mathcal{X}$ have
_not_ been observed in $X^{n}$. By Eqn. (4) and (5), the expected value of
$M_{0}$ is

$\displaystyle\mathbb{E}\left[M_{0}\right]$ $\displaystyle=g_{1}(n+1)$
$\displaystyle=\left[\sum_{k=1}^{n}(-1)^{k-1}g_{k}(n)\right]+(-1)^{n}f_{n+1}(n+1).$

The values in the second column of Figure 1 ($k=1$) represents the expected
values of missing mass; $\mathbb{E}\left[M_{0}\right]$ being the cumulative
sum of $(-1)^{k-1}g_{k}(n)$ is intuitively clear from the figure (the red cell
in the row $n=7$). It is here where we observe that
$\mathbb{E}\left[M_{0}\right]=g_{1}(n+1)$ is _almost entirely determined_ by
the $g_{*}(n)$, the expected frequencies of frequencies in the sample $X^{n}$,
and _not_ by the number of elements $|\mathcal{X}|$ or their underlying
distribution $p$. In fact, the influence of $p$ in the remainder term decays
exponentially, i.e.,
$f_{n+1}(n+1)=\sum_{x\in\mathcal{X}}p_{x}^{n+1}\leq\sum_{x\in\mathcal{X}}\left(e^{1-p_{x}}\right)^{-n-1}$
which is dominated by the discovery probability of the most abundant element
$\max(p)$.

Total Mass. Similarly, the expected value of the total probability mass
$\mathbb{E}\left[M_{k}\right]$ (the red cell in the row $n=10$), which is
equal to ${n\choose k}g_{k+1}(n+1)$, is almost entirely determined by the
expected frequencies of the sample $X^{n}$ with remainder $R_{n,k}={n\choose
k}\sum_{x\in\mathcal{X}}p_{x}^{n+1}$.

##  3 A Large Class of Estimators

From the representation of $\mathbb{E}\left[M_{k}\right]$ in terms of
frequencies in Eqn. (4) and the relationship across frequencies in Eqn. (3),
we can see that there is a large number of representations of the expected
total probability mass $\mathbb{E}\left[M_{k}\right]$. Each representation
might suggest different estimators.

We start by defining the minimal bias estimator $\hat{M}_{k}^{B}$ from the
representation in Eqn. (5) and explore its properties.

###  3.1 Estimator with Exponentially Decaying Bias

Let

$\hat{M}_{k}^{B}=-{n\choose
k}\sum_{i=1}^{n-k}\frac{(-1)^{i}\Phi_{k+i}}{{n\choose k+i}}$

Bias. For some constant $k:0\leq k\leq n$ and some constant $\mathbf{c}>1$,
the bias of $\hat{M}_{k}^{B}$ is in the order of
$\mathcal{O}(n^{k}\mathbf{c}^{-n})$, i.e.,

$\displaystyle\mathbb{E}\left[\hat{M}_{k}^{B}-M_{k}\right]$
$\displaystyle=-R_{n,k}=(-1)^{n-k-1}{n\choose
k}\sum_{x\in\mathcal{X}}p_{x}^{n+1}$
$\displaystyle\left|\mathbb{E}\left[\hat{M}_{k}^{B}-M_{k}\right]\right|$
$\displaystyle\leq{n\choose k}\sum_{x\in\mathcal{X}}\mathbf{c}_{x}^{-n}\leq
n^{k}\sum_{x\in\mathcal{X}}\mathbf{c}_{x}^{-n}$

where $\mathbf{c}_{x}>1$ for all $x\in\mathcal{X}$ are constants.

Variance. The variance of $\hat{M}_{k}^{B}$ is given by the variances and
covariances of the frequencies $\Phi_{k+i}$ for $i=1,\cdots,n-k$. Under the
certain conditions, the variance of $\hat{M}_{k}^{B}$ also decays
exponentially in $n$.

######  Theorem 3.1.

The variance of $\hat{M}_{k}^{B}$ decreases exponentially with $n$ if
$p_{\max}<0.5$ or $\frac{(1-p_{\max})(1-p_{\min})}{p_{\max}}<1$, where
$p_{\max}=\max_{x\in\mathcal{X}}p_{x}$ and
$p_{\min}=\min_{x\in\mathcal{X}}p_{x}$.

The proof is postponed to Appendix B.

Comparison to Good-Turing (GT). The bias of $\hat{M}_{k}^{B}$ not only decays
exponentially in $n$ but is also _smaller_ than that of GT estimator
$\hat{M}_{k}^{G}$ by an exponential factor. For a simpler variant of GT
estimator, $\hat{M}_{k}^{G^{\prime}}=\frac{k+1}{n-k}\Phi_{k+1}$ (suggested in
(McAllester & Schapire, 2000)), which corresponds to the first term in the
expected total probability mass $\mathbb{E}\left[M_{k}\right]$ in Eqn. (5), we
show that its bias is _larger by an exponential factor_ than the absolute bias
of $\hat{M}_{k}^{B}$. To see this, we provide bounds on the individual sums
and then on the bias ratio:

$\displaystyle\text{\emph{Bias}}_{G^{\prime}}=\mathbb{E}\left[\hat{M}_{k}^{G^{\prime}}-M_{k}\right]$
$\displaystyle\geq{n\choose k}p_{\min}^{k+2}(1-p_{\min})^{n-k-1}$ (6)
$\displaystyle|\text{\emph{Bias}}_{B}|=\left|\mathbb{E}\left[\hat{M}_{k}^{B}-M_{k}\right]\right|$
$\displaystyle\leq{n\choose k}Sp_{\max}^{n+1}$ (7)

where $S=|\mathcal{X}|$, such that

$\displaystyle\left|\frac{\text{\emph{Bias}}_{G^{\prime}}}{\text{\emph{Bias}}_{B}}\right|$
$\displaystyle\geq\frac{p^{k+2}_{\min}}{Sp^{k+2}_{\max}}\left(\frac{1-p_{\min}}{p_{\max}}\right)^{n-k-1}.$

Notice that $(1-p_{\min})/p_{\max}>1$ for all distributions over
$\mathcal{X}$, except where $S=2$ and $p=\\{0.5,0.5\\}$. The same can be shown
for the original GT estimator $\hat{M}_{k}^{G}=\frac{k+1}{n}\Phi_{k+1}$ for a
sufficiently large sample size (see Appendix A).

Example (Missing mass for the uniform). Suppose, we seek to estimate the
missing mass from a sequence of elements $X^{n}$ observed in $n$ samples from
the uniform distribution; $p_{x}=1/S$ for all $x\in\mathcal{X}$.
$\hat{M}_{0}^{G}$ overestimates $M_{0}$ on the average by $(S-1)^{n-1}/S^{n}$
while our new estimator $\hat{M}_{0}^{B}$ has a bias of $(-1)^{n}/S^{n}$.
Hence, for the uniform distribution, our estimator exhibits a bias that is
lower by a factor of $1/(S-1)^{n-1}$.

While the bias of our estimator $\hat{M}_{k}^{B}$ is lower than that of
$\hat{M}_{0}^{G}$ by an exponential factor, the variance is higher. The
variance of $\hat{M}_{k}^{B}$ depends on the variances of and covariances
between $\Phi_{k+1},\Phi_{k+2},\ldots,\Phi_{n}$, i.e.,

$\displaystyle\begin{split}\mathrm{Var}\left(\hat{M}_{k}^{B}\right)&=\sum_{i=1}^{n-k}c_{i}^{2}\mathrm{Var}\left(\Phi_{k+i}\right)\\\
&+\sum_{i\neq
j}(-1)^{i+j}c_{i}c_{j}\mathrm{Cov}\left(\Phi_{k+i},\Phi_{k+j}\right),\end{split}$
(8)

where $c_{i}=\left.\binom{n}{k}\right/\binom{n}{k+i}$. In contrast, the
variance of $\hat{M}_{k}^{G}$ depends only on the variance of $\Phi_{k+1}$. In
the empirical study, we investigate the difference of the two estimators in
terms of the bias and the variance.

###  3.2 Estimation with Minimal MSE as Search Problem

There are many representations of $\mathbb{E}\left[M_{k}\right]={n\choose
k}g_{k+1}(n+1)$ that can be constructed by recursively rewriting terms
according to the dependency among frequencies we identified (cf. Eqn. (3 &
4)). The representation used to construct our minimal-bias estimator
$\hat{M}_{k}^{B}$ was one of them. However, we notice that the variance of
$\hat{M}_{k}^{B}$ is too high to be practical. To find a representation from
which an estimator with a minimal mean squared error (MSE) can be derived, we
cast the efficient estimation of $M_{k}$ as an _optimization problem_. To
efficiently navigate the large search space of representations of
$\mathbb{E}\left[M_{k}\right]$, we develop a genetic algorithm.

Search space. Let $\mathbb{E}\left[M_{k}\right]$ be _represented_ by a
suitable choice of coefficients $\alpha_{i,j}$ such that

$\displaystyle\mathbb{E}\left[M_{k}\right]=\sum_{i=1}^{n+1}\sum_{j=i}^{n+1}\alpha_{i,j}g_{i}(j).$
(9)

One representation of $\mathbb{E}\left[M_{k}\right]={n\choose k}g_{k+1}(n+1)$
is

$\displaystyle r_{0}=\left\\{\alpha_{i,j}=\begin{cases}\binom{n}{k}&\text{for
$i=k+1$ and $j={n+1}$}\\\ 0&\text{otherwise.}\hfill\end{cases}\right\\}$ (10)

Mutation. Given any representation $r$ of $\mathbb{E}\left[M_{k}\right]$, we
can construct a new representation $r^{\prime}$ of
$\mathbb{E}\left[M_{k}\right]$, s.t. Eqn. (9) holds by _recursively_
considering the following identities:

$\displaystyle\alpha_{i,j}\cdot g_{i}(j)$
$\displaystyle=\alpha_{i,j}\cdot((1-\delta)g_{i}(j)+\delta g_{i}(j))$ (11)
$\displaystyle=\alpha_{i,j}\cdot(g_{i}(j+1)+g_{i+1}(j+1))$ (12)
$\displaystyle=\alpha_{i,j}\cdot(g_{i}(j-1)-g_{i+1}(j))$ (13)
$\displaystyle=\alpha_{i,j}\cdot(g_{i-1}(j-1)+g_{i-1}(j))$ (14)

for any choice of $\delta:0\leq\delta\leq 1$. _Importantly_ , after applying
these identities, we must work out the new coefficients accordingly. For
instance, applying Eqn. (11) with $\delta=0.5$ and Eqn. (13) to $r_{0}$, we
obtain the following representation $r_{1}$ of $\mathbb{E}\left[M_{k}\right]$:

$\displaystyle
r_{1}=\left\\{\alpha_{i,j}=\begin{cases}\left.\binom{n}{k}\right/2&\text{for
$i=k+1$ and $j={n+1}$}\\\ \left.\binom{n}{k}\right/2&\text{for $i=k+1$ and
$j={n}$}\\\ \left.-\binom{n}{k}\right/2&\text{for $i=k+2$ and $j={n+1}$}\\\
0&\text{otherwise.}\hfill\end{cases}\right\\}$

Algorithm 1 Genetic Algorithm

0: Target frequency $k$, Sample $X^{n}$

0: Iteration limit $G$, mutant size $m$

1: Population $P_{0}=\\{r_{0}\\}$

2: Fitness $f^{\text{best}}=f_{0}=\text{\emph{fitness}}(r_{0})$

3: Limit $G_{L}=G$

4: for $g$ from $1$ to $G_{L}$ do

5: $P=\text{\emph{selectTopM}}(P_{g-1},m)$

6: $P^{\prime}=\text{\emph{lapply}}(P,\text{\emph{mutate}})$

7: $P_{g}=P^{\prime}\cup\\{r_{0}\\}\cup\text{\emph{selectTopM}}(P_{g-1},3)$

8: $f_{g}=\min(\text{\emph{lapply}}(P_{g},\text{\emph{fitness}}))$

9: if $(g=G_{L})\wedge((f_{g}=f_{0})\vee(f^{\text{best}}>0.95\cdot f_{g}))$
then

10: $G_{L}=G_{L}+G$

11: $f^{\text{best}}=f_{g}$

12: end if

13: end for

14: Estimator
$\hat{M}_{k}^{\text{Evo}}=\text{\emph{instantiate}}(\text{\emph{selectTopM}}(P_{G_{L}},1))$

14: Minimal-MSE Estimator $\hat{M}_{k}^{\text{Evo}}$

Estimator instantiation. To construct a unique estimator $\hat{M}_{k}^{r}$ of
$M_{k}$ from a representation $r$ of $\mathbb{E}\left[M_{k}\right]$, we
propose a deterministic method. But first, we define our random variables on
subsamples of $X^{n}$. For any $m\leq n$, let $N_{x}(m)$ be the number of
times element $x\in\mathcal{X}$ is observed in the subsample $X^{m}=\langle
X_{1},\cdots X_{m}\rangle$ of $X^{n}$. Let $\Phi_{k}(m)$ be the number of
elements appearing exactly $k$ times in $X^{m}$, i.e.,

$\displaystyle N_{x}(m)$ $\displaystyle=\sum_{i=1}^{m}\mathbbm{1}(X_{i}=x)$
$\displaystyle\text{Note that }N_{x}=N_{x}(n).$ $\displaystyle\Phi_{k}(m)$
$\displaystyle=\sum_{x\in\mathcal{X}}\mathbbm{1}(N_{x}(m)=k)$
$\displaystyle\text{Note that }\Phi_{k}=\Phi_{k}(n).$

Hence, given a representation $r$, we can construct $\hat{M}_{k}^{r}$ as

$\displaystyle\hat{M}_{k}^{r}=\left[\sum_{i=1}^{n}\sum_{j=i}^{n}\frac{\alpha_{i,j}}{{j\choose
i}}\Phi_{i}(j)\right]+\left[\sum_{i=1}^{n}\frac{\alpha_{i,n+1}}{{n+1\choose
i}}\Phi_{i}\right]$

Notice that $\Phi_{i}(j)\left/{j\choose i}\right.$ is just the plug-in
estimator for $g_{i}(j)$.

Fitness function. To define the quantity to optimize, any meta-heuristic
search requires a fitness function. Our _fitness function_ takes a candidate
representation $r$ and returns an estimate of the MSE of the corresponding
estimator $\hat{M}_{k}^{r}$. We decompose the MSE as the sum of its variance
and squared bias. For convenience, let $g_{n+1}(n)=0$.

$\displaystyle\text{MSE}(\hat{M}_{k}^{r})$
$\displaystyle=\left[\sum_{i=1}^{n+1}\alpha_{i,n+1}[g_{i}(n+1)-g_{i}(n)]\right]^{2}$
(15)
$\displaystyle+\sum_{i=1}^{n}\sum_{j=i}^{n}\left(\frac{\alpha_{i,j}}{{j\choose
i}}\right)^{2}\mathrm{Var}\left(\Phi_{i}(j)\right)$
$\displaystyle+\sum_{i=1}^{n}\sum_{j=i}^{n}\sum_{l=1\atop l\neq
i}^{n}\sum_{m=l\atop m\neq j}^{n}\frac{\alpha_{i,j}}{{j\choose
i}}\frac{\alpha_{l,m}}{{m\choose
l}}\mathrm{Cov}\left(\Phi_{i}(j),\Phi_{l}(m)\right)$

We expand on the computation of the MSE in Appendix C.

(a) $M_{0}$ ($S=1000$).

(b) $M_{0}$ ($S=1000$, cropped).

(c) $M_{k}$ ($S=1000$, $n=2000$).

$n$ | $\hat{M}_{0}$ | Bias | Var | MSE
---|---|---|---|---
100 | GT | 3.6973e-003 | 2.3372e-03 | 2.3508e-03
ours | 1.0000e-200 | 2.3515e-03 | 2.3515e-03
500 | GT | 6.6369e-005 | 1.1430e-05 | 1.1434e-05
ours |  $<$ 1.00e-200 | 1.1445e-05 | 1.1445e-05
1000 | GT | 4.3607e-007 | 4.3439e-08 | 4.3439e-08
ours |  $<$ 1.00e-200 | 4.3441e-08 | 4.3441e-08

(d) $M_{0}$ (uniform, $S=100$).

Figure 2: Results for our minimal-bias estimator $\hat{M}_{k}^{B}$. For our
and baseline estimators, we show the logarithm of the absolute bias (a&b) as a
function of $n$ for $k=0$ and (c) as a function of $k$ for $n=2000$. We also
show (d) the bias, variance, and MSE of our and baseline estimator for three
values of $n$. More plots can be found in Appendix E.

Since the underlying distribution $\\{p_{x}\\}_{x\in\mathcal{X}}$ is _unknown_
, we can only _estimate_ the MSE. For any element $x$ that has been observed
exactly $k>0$ time in the sample $X^{n}$, we use
$\hat{p}_{x}=\hat{M}_{k}^{G}/\Phi_{k}$ as natural estimator of $p_{x}$, where
$\hat{M}_{k}^{G}$ is the GT estimator. To handle unobserved elements ($k=0$),
we first estimate the number of unseen elements
$\mathbb{E}\left[\Phi_{0}\right]=f_{0}(n)$ using Chao’s nonparamteric species
richness estimator $\hat{f}_{0}=\frac{n-1}{n}\frac{\Phi_{1}^{2}}{2\Phi_{2}}$
(Chao, 1984), and then estimate the probability of each such unseen element as
$\hat{p}_{y}=\hat{M}_{0}^{G}/\hat{f}_{0}$, where $\hat{M}_{0}^{G}$ is the GT
estimator. Finally, we plug these estimates into Eqn. (15) to estimate the
MSE. It is interesting to note that it is precisely the GT estimator whose MSE
our approach is supposed to improve upon.

Genetic algorithm. With the required concepts in place, we are ready to
introduce our _genetic algorithm_ (GA) (Mitchell, 1998). Algorithm 1 sketches
the general procedure. Given a target frequency $k$ (incl. $k=0$), the sample
$X^{n}$, an iteration limit G, and the number $m$ of candidate representations
to be mutated in every iteration, the algorithm produces an estimator
$\hat{M}_{k}^{\text{Evo}}$ with minimal MSE.

Starting from the _initial representation_ $r_{0}$ (Eqn. (10); Line 1), our GA
iteratively improves a population of candidate representations $P_{g}$, called
_individuals_. For every generation $g$ (Line 4), our GA selects the $m$
fittest individuals from the previous generation $P_{g-1}$ (Line 5), mutates
them (Line 6), and creates the current generation $P_{g}$ by adding the
initial representation $r_{0}$ and the Top-3 individuals from the previous
generation (Line 7). The initial and previous Top-3 individuals are added to
mitigate the risk of convergence to a local optimum. To _mutate_ a
representation $r$, our GA (i) chooses a random term $r$, (ii) applies Eqn.
(11) where $\delta$ is chosen uniformly at random, (iii) applies a random
identity from Eqn. (12–14), and (iv) adjusts the coefficients for the
resulting representation $r^{\prime}$ accordingly. The iteration limit $G_{L}$
is increased if the current individuals do _not_ improve on the initial
individual $r_{0}$ or _substantially_ improve on those discovered recently
(Line 9–12).

Distribution-free. While our approach itself is _distribution-free_ , the
output is _distribution-specific_ , i.e., the discovered estimator has a
minimal MSE on the specific, unknown distribution. We explore this property in
our experiments.

##  4 Experiment

We design experiments to evaluate the performance (i) of our minimal-bias
estimator $\hat{M}_{k}^{B}$ and (ii) of our the minimal-MSE estimator
$\hat{M}_{k}^{\text{Evo}}$ that is discovered by our genetic algorithm against
the performance of the widely-used Good-Turing estimator $\hat{M}_{k}^{G}$
(Good, 1953).

Distibutions. We use the same six multinomial distributions that are used in
previous evaluations (Orlitsky & Suresh, 2015; Orlitsky et al., 2016; Hao &
Li, 2020): a uniform distribution (uniform), a half-and-half distribution
where half of the elements have three times of the probability of the other
half (half&half), two Zipf distributions with parameters $s=1$ and $s=0.5$
(zipf-1, zipf-0.5, respectively), and distributions generated by Dirichlet-1
prior and Dirichlet-0.5 prior (diri-1, diri-0.5, respectively).

Open Science and Replication. For scrutiny and replicability, we publish all
our evaluation scripts at:

https://anonymous.4open.science/r/Better-Turing-157F.

###  4.1 Evaluating our Minimal-Bias Estimator

  * •

RQ1. _How does our estimator for the missing mass $\hat{M}_{0}^{B}$ compare to
the Good-Turing estimator $\hat{M}_{0}^{G}$ in terms of bias as a function of
sample size $n$?_

  * •

RQ2. _How does our estimator for the total mass $\hat{M}_{k}^{B}$ compare to
the Good-Turing estimator $\hat{M}_{k}^{G}$ in terms of bias as a function of
frequency $k$?_

  * •

RQ3. _How do the estimators compare in terms of variance and mean-squared
error?_

We focus specifically on the _bias_ of $\hat{M}_{k}^{B}$, i.e., the average
difference between the estimate and the expected value
$\mathbb{E}\left[M_{k}\right]$. We expect that the bias of the _missing mass_
estimate $\hat{M}_{0}^{B}$ as a function of $n$ across different distributions
provides empirical insight for our claim that how much is unseen chiefly
depends on information about the seen.

Table 1:  The MSE of the best evolved estimator $M_{0}^{\text{Evo}}$ and GT
estimator $\hat{M}_{0}^{G}$ for the missing mass $M_{0}$, the success rate
$\hat{A}_{12}$, and the ratio (Ratio,
$\mathit{MSE}(\hat{M}_{0}^{\text{Evo}})/\mathit{MSE}(M_{0}^{G})$) for three
sample sizes $n$ and six distributions with support size $S=200$.

Dist. | $n=S/2$ | $n=S$ | $n=2S$
---|---|---|---
$\mathit{MSE}(\hat{M}_{0}^{G})$ | $\mathit{MSE}(M_{0}^{\text{Evo}})$ | $\hat{A}_{12}$ | Ratio | $\mathit{MSE}(\hat{M}_{0}^{G})$ | $\mathit{MSE}(M_{0}^{\text{Evo}})$ | $\hat{A}_{12}$ | Ratio | $\mathit{MSE}(\hat{M}_{0}^{G})$ | $\mathit{MSE}(M_{0}^{\text{Evo}})$ | $\hat{A}_{12}$ | Ratio
uniform | 3.32e-03 | 2.04e-03 | 0.95 | 61% | 1.17e-03 | 8.90e-04 | 0.99 | 76% | 2.01e-04 | 1.70e-04 | 0.93 | 84%
half&half | 3.33e-03 | 1.97e-03 | 0.96 | 59% | 1.09e-03 | 8.58e-04 | 0.99 | 78% | 2.11e-04 | 1.72e-04 | 1.00 | 81%
zipf-1 | 2.32e-03 | 2.41e-03 | 0.74 | 103% | 8.16e-04 | 7.24e-04 | 0.88 | 88% | 2.39e-04 | 2.11e-04 | 0.96 | 88%
zipf-0.5 | 3.23e-03 | 2.29e-03 | 0.89 | 71% | 1.09e-03 | 8.52e-04 | 0.97 | 78% | 2.30e-04 | 1.93e-04 | 1.00 | 83%
diri-1 | 2.99e-03 | 2.36e-03 | 0.85 | 78% | 8.88e-04 | 6.65e-04 | 1.00 | 74% | 1.96e-04 | 1.65e-04 | 0.96 | 84%
diri-0.5 | 2.55e-03 | 1.81e-03 | 0.94 | 71% | 6.88e-04 | 4.86e-04 | 0.98 | 70% | 1.61e-04 | 1.31e-04 | 0.93 | 81%
Avg. |  |  | 0.88 | 74% |  |  | 0.96 | 77% |  |  | 0.96 | 84%

RQ.1. 2(a) illustrates how fast our estimator $\hat{M}_{k}^{B}$ and the
baseline estimator $\hat{M}_{k}^{G}$ (GT) approach the expected missing mass
$\mathbb{E}\left[M_{0}\right]$ as a function of sample size $n$. As it might
difficult for the reader to discern differences across distributions for the
baseline estimator, we refer to 2(b), where we zoom into a relevant region.

The _magnitude_ of our estimator’s bias is significantly smaller than the
magnitude of GT’s bias _for all distributions_ (by thousands of orders of
magnitude).333Recall that the plot shows the _logarithm_ of the absolute bias.
Figure 2(a) also nicely illustrates the _exponential decay_ of our estimator
in terms of $n$ and how our estimator is less biased than GT by an exponential
factor. In Figure 2(b), we can observe that GT’s bias also decays
exponentially, although not nearly at the rate of our estimator.

In terms of distributions, a closer look at the performance differences
confirms our suspicion that the bias of our estimator is strongly influenced
by the probability $p_{\max}$ of the most abundant element while the bias of
GT is strongly influenced by the probability $p_{\min}$ of the rarest element.
In fact, by Eqn. (7) the absolute bias of our estimator is minimized when
$p_{\max}$ is minimized. By Eqn. (6), GT’s bias is minimized if $p_{\min}$ is
maximized. Since both is true for the uniform, both estimators exhibit the
lowest bias for the uniform across all six distributions. GT performs similar
on all distributions apart from the uniform (where bias seems minimal) and
zipf-1 (where bias is maximized). For our estimator, if we ranked the
distributions by values of $p_{\max}$ with the smallest value first
$\langle\text{{\small{uniform}}, {\small{half\&half}}, {\small{zipf-0.5}},
{\small{zipf-1}}}\rangle$,444diri-1 and diri-0.5 are not considered in the
order because multiple distributions are sampled from the Dirichlet prior. we
would arrive at the same ordering in terms of performance of our estimator as
shown in 2(a).

RQ2. Figure 2(c) illustrates for both estimators of the _total mass_ $M_{k}$
how the bias behaves as $k$ varies between $0$ and $n=2000$ when $S=1000$. The
trend is clear; the bias of our estimator is strictly smaller than the bias of
GT for all $k$ and all the distributions. The difference is the most
significant for rare elements (small $k$) and gets smaller as $k$ increases.
The bias of our estimator is maximized when $k=1000=0.5n$, the bias for GT
when $k=0$.

RQ3. Table 2(d) shows variance and MSE of both estimators for the missing mass
$M_{0}$ for the uniform and three values of $n$. As we can see, the MSE of our
estimator is approximately the same as that of GT. The reason is that the MSE
is dominated by the variance. We make the same observation for all other
distributions (see Appendix E). The MSEs of both estimators are comparable.

###  4.2 Evaluating our Estimator Discovery Algorithm

  * •

RQ1 (Effectiveness). _How does our estimator for the missing mass
$\hat{M}_{0}^{\text{Evo}}$ compare to the Good-Turing estimator
$\hat{M}_{0}^{G}$ in terms of MSE?_

  * •

RQ2 (Efficiency). _How long does it take for our genetic algorithm to generate
an estimator $M_{k}^{\text{Evo}}$ given a sample?_

  * •

RQ3 (Distribution-awareness). _How well does an estimator discovered from a
sample from one distribution perform on another distribution in terms of MSE?_

To handle the randomness in our evaluation, we repeat each experiment 100
times: 10 runs of the GA with 10 different samples $X^{n}$.555 For diri-1,
diri-0.5, each of the ten samples $X^{n}$ is sampled from 10 distributions
sampled from the Dirichlet prior with the same parameter $\alpha=0.5,1$,
respectively.  More details about our experimental setup can be found in
Appendix D.

RQ.1 (Effectiveness). Table 1 shows average MSE of the estimator
$M_{0}^{\text{Evo}}$ discovered by our genetic algorithm and that of the GT
estimator $\hat{M}_{0}^{G}$ for the missing mass $M_{0}$ across three sample
sizes. We measure effect size using _Vargha-Delaney_ $\hat{A}_{12}$ (Vargha &
Delaney, 2000) (_success rate_), i.e., the probability that the MSE of the
estimator discovered by our genetic algorithm has a smaller MSE than the GT
estimator (_larger is better_). Moreover, we measure the MSE of our estimator
as a proprtion of the MSE of GT, called _ratio_ (_smaller is better_). Results
for other support sizes $S$ can be found in Appendix E.

Overall, the estimator discovered by our GA performs _significantly_ better
than GT estimator in terms of MSE (avg. $\hat{A}_{12}>0.9$; $ratio<85\%$). The
performance difference increases with sample size $n$. When the sample size is
twice the support size ($n=2S$), in 96% of runs our discovered estimator
performs better. The average MSE of our estimator is somewhere between 70% and
88% of the MSE of GT. The high success rate and the low ratio of the MSE shows
that the GA is effective in finding the estimator with the minimal MSE for the
missing mass $M_{0}$. A Wilcoxon signed-rank test shows that all performance
differences are statistically significant at $\alpha<10^{-9}$.

In terms of distributions, the performance of our estimator is similar across
all distributions, showing the generality of our algorithm. The only exception
is the zipf-1, where the success rate is lower than for other distribution for
$n=S/2$ and $S$, and the average ratio is 103% (yet, the median ratio is 85%)
for $n=S/2$. The potential reason for this is due to the overfitting to the
approximated distribution $\hat{p}_{x}$. Since the zipf-1 is the most skewed
distribution, there are more elements unseen in the sample than in other
distributions, which makes the approximated distribution $\hat{p}_{x}$ less
accurate. Yet, the performance already improves and become similar to other
distributions for $n=S$ and $n=2S$.

RQ.2 (Efficiency). The time it takes to discover the estimator is reasonable.
To compute an estimator fo Table 1, it takes about seven (7) minutes on
average and five (5) minutes on median. The average time for each iteration is
1.25s (median: 0.92).

Figure 3: The MSE of an estimator discovered using a sample from one
distribution (individual boxes) applied to another target distribution
(clusters of boxes).

RQ.3 (Distribution-awareness). Figure 3 shows the performance of an estimator
discovered from a sample from one distribution (source) when applied to
another distribution (target). Applying an estimate from the zipf-1 on the
zipf-1 gives the optimal MSE (right-most red box). However, applying an
estimator from the zipf-1 on the uniform (left red box) yields a huge increase
in variance. In terms of effect size, we measure a Vargha Delaney
$\hat{A}_{12}>0.84$ between the “home” and “away” estimator. While each of the
uniform and half&half also shows that the home estimator performs best on the
home distribution ($\hat{A}_{12}=0.68$ (medium), $0.58$ (small),
respectively), the difference between the estimators from uniform, half&half,
and zipf-0.5 is less significant. Perhaps unsurprisingly, an estimator
performs optimal when the source of the samples is similar to the target
distribution.

Summary. To summarize, our GA is effective in finding the estimator with the
minimal MSE for the missing mass $M_{0}$ with the smaller MSE than GT
estimator $\hat{M}_{0}^{G}$ for all distributions and sample sizes. The effect
is substantial and significant and the average decrease of the MSE is roughly
one fifth against GT estimator $\hat{M}_{0}^{G}$.

##  5 Discussion

Beyond the General Estimator. In this study, we propose a meta-level
estimation methodology that can be applied to a set of samples from a specific
unknown distribution. The conventional approach is to develop _an estimator_
for an arbitrary distribution. Yet, each distributions has its own
characteristics, and, because of that, the manner of the (frequencies of)
frequencies of the classes in the sample can be differ from, for example, the
uniform distribution to the Zipf distribution. In contrast to the conventional
approach, we propose a _distribution-free_ methodology to discover the a
_distribution-specific_ estimator with low MSE (given only the sample). Note
that, while we use the genetic algorithm to discover the estimator, any
optimization method can be used to discover the estimator, for instance, a
constrained optimization solver.

Extrapolating the Future Sampling. Estimating the number of unseen species is
a well-known problem in many scientific fields, such as ecology, linguistics,
and machine learning. Given $n$ samples, the expected number of hitherto
unseen species that would be uncovered if $t$ times more samples were taken is
$\mathbb{E}\left[U(t)\right]=f_{0}(n)-f_{0}(n+nt)$. Good & Toulmin (1956)
proposed a seminal estimator using the frequencies of frequencies $\Phi_{k}$,
similar to the Good-Turing estimator. Until recently, various subsequent
studies have been conducted to improve the estimator (Efron & Thisted, 1976;
Orlitsky et al., 2016; Hao & Li, 2020), while most of them still relies on the
Poisson approximation to design the estimator. We believe that our analysis
can be extended to the Good-Toulmin estimator seeking more accurate estimators
for $U(t)$.

## References

  * Acharya et al. (2013) Acharya, J., Jafarpour, A., Orlitsky, A., and Suresh, A. T.  Optimal probability estimation with applications to prediction and classification.  In Shalev-Shwartz, S. and Steinwart, I. (eds.), _Proceedings of the 26th Annual Conference on Learning Theory_ , volume 30 of _Proceedings of Machine Learning Research_ , pp. 764–796, Princeton, NJ, USA, 12–14 Jun 2013. PMLR.
  * Chao (1984) Chao, A.  Nonparametric estimation of the number of classes in a population.  _Scandinavian Journal of statistics_ , pp. 265–270, 1984.
  * Chao & Jost (2012) Chao, A. and Jost, L.  Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size.  _Ecology_ , 93 12:2533–47, 2012.
  * Drukh & Mansour (2004) Drukh, E. and Mansour, Y.  Concentration bounds for unigrams language model.  In Shawe-Taylor, J. and Singer, Y. (eds.), _Learning Theory_ , pp. 170–185, Berlin, Heidelberg, 2004. Springer Berlin Heidelberg.  ISBN 978-3-540-27819-1.
  * Efron & Thisted (1976) Efron, B. and Thisted, R.  Estimating the number of unsen species: How many words did shakespeare know?  _Biometrika_ , 63(3):435–447, 1976.  ISSN 00063444.  URL http://www.jstor.org/stable/2335721.
  * Good (1953) Good, I. J.  The population frequencies of species and the estimation of population parameters.  _Biometrika_ , 40(3-4):237–264, 1953.
  * Good & Toulmin (1956) Good, I. J. and Toulmin, G. H.  The number of new species, and the increase in population coverage, when a sample is increased.  _Biometrika_ , 43(1/2):45–63, 1956.  ISSN 00063444.  URL http://www.jstor.org/stable/2333577.
  * Hao & Li (2020) Hao, Y. and Li, P.  Optimal prediction of the number of unseen species with multiplicity.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 8553–8564. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/618790ae971abb5610b16c826fb72d01-Paper.pdf.
  * Juang & Lo (1994) Juang, B.-H. and Lo, S.  On the bias of the turing-good estimate of probabilities.  _IEEE Transactions on signal processing_ , 2(42):496–498, 1994.
  * McAllester & Schapire (2000) McAllester, D. and Schapire, R. E.  On the convergence rate of good-Turing estimators.  In _Proceedings of the 13th Annual Conference on Computational Learning Theory_ , pp. 1–6. Morgan Kaufmann, San Francisco, 2000.
  * Mitchell (1998) Mitchell, M.  _An introduction to genetic algorithms_.  MIT press, 1998.
  * Orlitsky & Suresh (2015) Orlitsky, A. and Suresh, A. T.  Competitive distribution estimation: Why is good-turing good.  In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 28. Curran Associates, Inc., 2015.  URL https://proceedings.neurips.cc/paper/2015/file/d759175de8ea5b1d9a2660e45554894f-Paper.pdf.
  * Orlitsky et al. (2003) Orlitsky, A., Santhanam, N. P., and Zhang, J.  Always good turing: Asymptotically optimal probability estimation.  _Science_ , 302(5644):427–431, 2003.  doi: 10.1126/science.1088284.
  * Orlitsky et al. (2016) Orlitsky, A., Suresh, A. T., and Wu, Y.  Optimal prediction of the number of unseen species.  _Proceedings of the National Academy of Sciences_ , 113(47):13283–13288, 2016.  doi: 10.1073/pnas.1607774113.  URL https://www.pnas.org/doi/abs/10.1073/pnas.1607774113.
  * Painsky (2022) Painsky, A.  Convergence guarantees for the good-turing estimator.  _Journal of Machine Learning Research_ , 23(279):1–37, 2022.  URL http://jmlr.org/papers/v23/21-1528.html.
  * Valiant & Valiant (2016) Valiant, G. and Valiant, P.  Instance optimal learning of discrete distributions.  In _Proceedings of the Forty-Eighth Annual ACM Symposium on Theory of Computing_ , STOC ’16, pp. 142–155, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450341325.  doi: 10.1145/2897518.2897641.  URL https://doi.org/10.1145/2897518.2897641.
  * Vargha & Delaney (2000) Vargha, A. and Delaney, H. D.  A critique and improvement of the cl common language effect size statistics of mcgraw and wong.  _Journal of Educational and Behavioral Statistics_ , 25(2):101–132, 2000.

##  Appendix A Comparing the Bias of the Estimators

In Section 3.1, we have shown that the bias of a simpler variant of GT,
$\hat{M}_{k}^{G^{\prime}}=\frac{k+1}{n-k}\Phi_{k+1}$, is larger by an
exponential factor than the absolute bias of our minimal bias estimator
$\hat{M}_{k}^{B}$. In this section, we show that the bias of the original GT
estimator $\hat{M}_{k}^{G}=\frac{k+1}{n}\Phi_{k+1}$ is also larger by an
exponential factor than the absolute bias of $\hat{M}_{k}^{B}$ for a
sufficiently larger sample size. Recall that

$\text{\emph{Bias}}_{G^{\prime}}=\mathbb{E}\left[\hat{M}_{k}^{G}-M_{k}\right]=\frac{k+1}{n-k}f_{k+1}(n)-{n\choose
k}g_{k+1}(n+1)=\sum_{x}{n\choose k}p_{x}^{k+2}(1-p_{x})^{n-k-1},\\\ $ (16)

and

$\displaystyle\text{\emph{Bias}}_{G}=\mathbb{E}\left[\hat{M}_{k}^{G}-M_{k}\right]$
$\displaystyle={n\choose k}g_{k+2}(n+1)-\frac{k(k+1)}{n(n-k)}f_{k+1}(n)$ (17)
$\displaystyle={n\choose k}g_{k+2}(n+1)-\binom{n-1}{k-1}g_{k+1}(n)$ (18)
$\displaystyle=\sum_{x}p_{x}^{k+2}(1-p_{x})^{n-k-1}\left(\binom{n}{k}-\frac{1}{p_{x}}\cdot\binom{n-1}{k-1}\right)$
(19)
$\displaystyle=\sum_{x}\binom{n}{k}p_{x}^{k+2}(1-p_{x})^{n-k-1}\left(1-\frac{k}{n\cdot
p_{x}}\right)$ (20) $\displaystyle\geq\left(1-\frac{k}{n\cdot
p_{\max}}\right)\text{\emph{Bias}}_{G^{\prime}},$ (21)

where $1-\frac{k}{n\cdot p_{min}}>0$ when $n$ is sufficiently large. Above
inequality leads to the following:

$\frac{\text{\emph{Bias}}_{G}}{\text{\emph{Bias}}_{G^{\prime}}}\geq\left(1-\frac{k}{n\cdot
p_{\max}}\right)\text{,}\quad\text{while}\quad\frac{\text{\emph{Bias}}_{B}}{\text{\emph{Bias}}_{G^{\prime}}}\leq\frac{Sp_{\max}^{k+2}}{p_{\min}^{k+2}}\left(\frac{1-p_{\min}}{p_{\max}}\right)^{-n+k+1},$
(22)

which proves our claim.

##  Appendix B Bounding the Variance and the MSE of $\hat{M}_{k}^{B}$

The MSE of an estimator $\hat{e}$ for an estimand $e$ is defined as
$MSE(\hat{e})=\mathbb{E}[(\hat{e}-e)^{2}]=\mathrm{Var}\left(\hat{e}\right)+\text{Bias}^{2}(\hat{e})$.
As we have shown the bias of $\hat{M}_{k}^{B}$ in Section 3.1, the remaining
part to compute the MSE of $\hat{M}_{k}^{B}$ is to compute its variance.

The variance of the linear combination of random variables is given by

$\mathrm{Var}\left(\sum_{i}c_{i}X_{i}\right)=\sum_{i}c_{i}^{2}\mathrm{Var}\left(X_{i}\right)+\sum_{i\neq
j}c_{i}c_{j}\mathrm{Cov}\left(X_{i},X_{j}\right).$ (23)

Therefore, the variance and the covariance of $\Phi_{i}(n)$s are the missing
pieces to compute the variance of $\hat{M}_{k}^{B}$.

######  Theorem B.1.

Given the multinomial distribution $p=(p_{1},\dots,p_{S})$ with support size
$S$, the variance of $\Phi_{i}=\Phi_{i}(n)$ from $n$ samples $X^{n}$ is given
by

$\mathrm{Var}\left(\Phi_{i}(n)\right)=\begin{cases}f_{i}(n)-f_{i}(n)^{2}+\sum_{x\neq
y}\frac{n!}{i!^{2}(n-2i)!}p_{x}^{i}p_{y}^{i}(1-p_{x}-p_{y})^{n-2i}&\text{if
$2i\leq n$,}\\\ f_{i}(n)-f_{i}(n)^{2}&\text{otherwise.}\end{cases}$ (24)

###### Proof.

$\displaystyle\mathrm{Var}\left(\Phi_{i}\right)$
$\displaystyle=\mathbb{E}\left[\Phi_{i}^{2}\right]-\mathbb{E}\left[\Phi_{i}\right]^{2}$
(25) $\displaystyle\mathbb{E}\left[\Phi_{i}^{2}\right]$
$\displaystyle=\mathbb{E}\left[\left(\sum_{x}\mathbbm{1}(N_{x}=i)\right)^{2}\right]$
(26) $\displaystyle=\mathbb{E}\left[\sum_{x}\mathbbm{1}(N_{x}=i)+\sum_{x\neq
y}\mathbbm{1}(N_{x}=i\land N_{y}=i)\right]$ (27)
$\displaystyle=\begin{cases}f_{i}(n)+\sum_{x\neq
y}\frac{n!}{i!^{2}(n-2i)!}p_{x}^{i}p_{y}^{i}(1-p_{x}-p_{y})^{n-2i}&\text{if
$2i\leq n$,}\\\ f_{i}(n)&\text{otherwise.}\end{cases}$ (28)
$\displaystyle\therefore\mathrm{Var}\left(\Phi_{i}\right)$
$\displaystyle=\begin{cases}f_{i}(n)+\sum_{x\neq
y}\frac{n!}{i!^{2}(n-2i)!}p_{x}^{i}p_{y}^{i}(1-p_{x}-p_{y})^{n-2i}-f_{i}(n)^{2}&\text{if
$2i\leq n$,}\\\ f_{i}(n)-f_{i}(n)^{2}&\text{otherwise.}\end{cases}$ (29)

∎

Now we compute the upper bound of the variance of $\hat{M}_{k}^{B}$.

######  Lemma B.2.

$\mathrm{Var}\left(\Phi_{i}\right)\begin{cases}\leq
Sf_{i}(n)-f_{i}(n)^{2}&\text{if $2i\leq n$.}\\\
=f_{i}(n)-f_{i}(n)^{2}&\text{otherwise.}\end{cases}$ (30)

###### Proof.

From Theorem B.1,

$\displaystyle\mathbb{E}\left[\Phi_{i}^{2}\right]$
$\displaystyle=f_{i}(n)+\mathbb{E}\left[\sum_{x\neq y}\mathbbm{1}(N_{x}=i\land
N_{y}=i)\right]$ (31) $\displaystyle\leq
f_{i}(n)+(S-1)\mathbb{E}\left[\sum_{x}\mathbbm{1}(N_{x}=i)\right]$ (32)
$\displaystyle=Sf_{i}(n)\quad\text{(if $2i\leq n$)}.$ (33)

The lemma directly follows from the above inequality. ∎

######  Lemma B.3.

$g_{i}(n)\leq S\cdot\beta_{\min}^{-n}o_{\max}^{i},$

where $S=|\mathcal{X}|$, $p_{\max}=\max_{x\in\mathcal{X}}p_{x}$,
$\beta_{\min}=\frac{1}{1-p_{\min}}$, and
$o_{\max}=\frac{p_{\max}}{1-p_{\max}}$.

###### Proof.

$\frac{1}{1-x}$ and $\frac{x}{1-x}$ are increasing functions for $x\in(0,1)$.
Therefore,

$g_{i}(n)=\sum_{x\in\mathcal{X}}p_{x}^{i}(1-p_{x})^{n-i}=\sum_{x\in\mathcal{X}}\left(\frac{1}{1-p_{x}}\right)^{-n}\left(\frac{p_{x}}{1-p_{x}}\right)^{i}\leq|\mathcal{X}|\cdot\beta_{\min}^{-n}o_{\max}^{i}.$

∎

######  Theorem B.4.

The variance of the estimator $\hat{M}_{k}^{B}$ is bounded as follows:

$\mathrm{Var}(\hat{M}_{k}^{B})\leq c_{1}\cdot n^{2k+1}\cdot c_{2}^{-n},$

where $c_{1}=S\cdot\left(\frac{e}{k}\right)^{2k}$,
$c_{2}=\min\left(\frac{1}{1-p_{\min}},\frac{1-p_{\max}}{p_{\max}(1-p_{\min})}\right)$.
In other words,
$\mathrm{Var}(\hat{M}_{k}^{B})=\mathcal{O}(n^{2k+1}\cdot\beta_{\min}^{-n}\cdot\max(1,o_{\max}^{n}))$.

###### Proof.

From Lemma B.2 in the supplementary, we have $\mathrm{Var}(\Phi_{i})\leq
Sf_{i}(n)-f_{i}(n)^{2}$. Thus,

$\displaystyle\mathrm{Var}\left(\frac{\Phi_{k+i}}{\binom{n}{k+i}}\right)$
$\displaystyle\leq\frac{Sf_{k+i}(n)-f_{k+i}(n)^{2}}{\binom{n}{k+i}^{2}}\leq
S\cdot\frac{g_{k+i}(n)}{\binom{n}{k+i}}\leq\frac{S^{2}\cdot\beta_{\min}^{-n}o_{\max}^{k+i}}{\binom{n}{k+1}}\quad\text{(by
Lemma~{}\ref{lemma:bound})}$
$\displaystyle\binom{n}{k}^{2}\mathrm{Var}\left(\frac{\Phi_{k+i}}{\binom{n}{k+i}}\right)$
$\displaystyle\leq
S\beta_{\min}^{-n}\frac{\binom{n}{k}^{2}}{\binom{n}{k+1}}o_{\max}^{k+i}\leq
S\beta_{\min}^{-n}\frac{\left(\frac{e^{2}n^{2}}{k^{2}}\right)^{k}}{\left(\frac{n}{k+i}\right)^{k}}o_{\max}^{k+i}\leq
S\beta_{\min}^{-n}\left(\frac{e^{2}n(k+i)}{k^{2}}\right)^{k}o_{\max}^{k+i}$
$\displaystyle\leq
S\beta_{\min}^{-n}\left(\frac{e^{2}n^{2}}{k^{2}}\right)^{k}O_{M}=S\beta_{\min}^{-n}\left(\frac{en}{k}\right)^{2k}O_{M},\quad\text{where
}O_{M}=\max(o_{\max}^{k},o_{\max}^{n}),$
$\displaystyle\mathrm{Var}(\hat{M}_{k}^{B})$
$\displaystyle=\binom{n}{k}^{2}\mathrm{Var}\left(\sum_{i=1}^{n-k}(-1)^{i-1}\frac{\Phi_{k+1}}{\binom{n}{k+1}}\right)$
$\displaystyle\leq(n-k)\binom{n}{k}^{2}\mathrm{Var}\left(\frac{\Phi_{k+1}}{\binom{n}{k+1}}\right)$
(35) $\displaystyle\leq
S(n-k)\left(\frac{en}{k}\right)^{2k}\cdot\beta_{\min}^{-n}O_{M}=\mathcal{O}(n^{2k+1})\cdot\beta_{\min}^{-n}O_{M},$

where (35) follows from Cauchy-Schwarz inequality
($\mathrm{Var}(\sum_{j=1}^{M}X_{i})\leq
M\cdot\sum_{j=1}^{M}\mathrm{Var}(X_{i})$). The proof follows from dividing the
variance of the estimator into two cases: $o_{\max}<1$ and $o_{\max}>1$: If
$o_{\max}<1$, $O_{M}=o_{\max}^{k}$, and
$\mathrm{Var}(\hat{M}_{k}^{B})=\mathcal{O}(n^{2k+1}\beta_{\min}^{-n})$. If
$o_{\max}>1$, $O_{M}=o_{\max}^{n}$, and,
$\mathrm{Var}(\hat{M}_{k}^{B})=\mathcal{O}(n^{2k+1}\beta_{\min}^{-n}o_{\max}^{n})$.
∎

Therefore, the variance exponentially decreases with $n$ if $p_{\max}<0.5$ or
$\frac{1-p_{\max}}{p_{\max}(1-p_{\min})}<1$.

######  Corollary B.5.

There exists a constant $c>1$ such that

$\mathrm{MSE}(\hat{M}_{k}^{B})\leq\mathcal{O}(n^{2k+1}c^{-n}).$

###### Proof.

From Equ. (7) in the manuscript, the bias
$|\mathbb{E}(\hat{M}_{k}^{B})-M_{k}|\leq S\cdot p_{\max}\cdot n^{k}\cdot
p_{\max}^{n}$. The proof follows from the fact that
$\mathrm{MSE}=\mathrm{Var}+\mathrm{Bias}^{2}$ and the bound of the variance
and the bias. ∎

##  Appendix C Computing the Variance and the MSE of the Evolved Estimators

Same as $\hat{M}_{k}^{B}$, the evolved estimators from the genetic algorithm
are also linear combinations of $\Phi_{k}(n)$s (while varying both $k$ and $n$
unlike $\hat{M}_{k}^{B}$). Given the evolved estimator
$\hat{M}_{k}^{\text{Evo}}=\sum_{i}c_{i}\Phi_{k_{i}}(n_{i})$, the expected
value of $\hat{M}_{k}^{\text{Evo}}$ is given by substituting $\Phi_{k}(n)$
with $f_{k}(n)$:

$\mathbb{E}(\hat{M}_{k}^{\text{Evo}})=\sum_{i}c_{i}f_{k_{i}}(n_{i}).$ (36)

Given the multinomial distribution $p$, the covariance between $\Phi_{k}(n)$
and $\Phi_{k^{\prime}}(n^{\prime})$, which is needed to compute the variance
of $\hat{M}_{k}^{\text{Evo}}$ as Equ. (23), can be computed as follows:

######  Theorem C.1.

Given the multinomial distribution $p=(p_{1},\dots,p_{S})$ with support size
$S$, let $X^{n_{total}}$ be the set of $n_{total}$ samples from $p$. Let
$X^{n}$ and $X^{n^{\prime}}$ be the first $n$ and $n^{\prime}$ samples from
$X^{n_{total}}$, respectively; WLOG, we assume $1\leq n^{\prime}\leq n\leq
n_{total}$. Then, the covariance of $\Phi_{k}(n)=\Phi_{k}(X^{n})$ and
$\Phi_{k^{\prime}}(n^{\prime})=\Phi_{k^{\prime}}(X^{n^{\prime}})$ ($1\leq
k\leq n$, $1\leq k^{\prime}\leq n^{\prime}$) is given by following:

$\displaystyle\mathrm{Cov}\left(\Phi_{k}(n),\Phi_{k^{\prime}}(n^{\prime})\right)$
$\displaystyle=\mathbb{E}\left[\Phi_{k}(n)\cdot\Phi_{k^{\prime}}(n^{\prime})\right]-f_{k}(n)\cdot
f_{k^{\prime}}(n^{\prime})$ (37)
$\displaystyle=\mathbb{E}\left[\left(\sum_{x}\mathbbm{1}(N_{x}=k)\right)\cdot\left(\sum_{x^{\prime}}\mathbbm{1}({N^{\prime}}_{x^{\prime}}=k^{\prime})\right)\right]-f_{k}(n)\cdot
f_{k^{\prime}}(n^{\prime})$ (38)
$\displaystyle=\sum_{x}\sum_{x^{\prime}}\mathbb{E}\left[\mathbbm{1}(N_{x}=k\land{N^{\prime}}_{x^{\prime}}=k^{\prime})\right]-f_{k}(n)\cdot
f_{k^{\prime}}(n^{\prime}),$ (39)

where ${N^{\prime}}_{x^{\prime}}$ is the number of occurrences of $x^{\prime}$
in $X^{n^{\prime}}$. Depending on the values of $n,n^{\prime},k,k^{\prime},x,$
and $x^{\prime}$, the
$\mathbb{E}\left[\mathbbm{1}(N_{x}=k\land{N^{\prime}}_{x^{\prime}}=k^{\prime})\right]$
can be computed as follows:

$\forall n,n^{\prime}$ _s.t._ |  $\forall x,x^{\prime}$ _s.t._ |  $\forall k,k^{\prime}$ _s.t._ | $\mathbb{E}\left[\mathbbm{1}(N_{x}=k\land{N^{\prime}}_{x^{\prime}}=k^{\prime})\right]$
---|---|---|---
$n=n^{\prime}$ | $x=x^{\prime}$ | $k=k^{\prime}$ | $\binom{n}{k}p_{x}^{k}(1-p_{x})^{n-k}$
$k\neq k^{\prime}$ |  $0$ (_infeasible_)
$x\neq x^{\prime}$ | $k+k^{\prime}\leq n$ | $\frac{n!}{k!k^{\prime}!(n-k-k^{\prime})!}p_{x}^{k}p_{x^{\prime}}^{k^{\prime}}(1-p_{x}-p_{x^{\prime}})^{n-k-k^{\prime}}$
$k+k^{\prime}>n$ |  $0$ (_infeasible_)
$n\neq n^{\prime}$ | $x=x^{\prime}$ | $k^{\prime}\leq k$ | $\binom{n^{\prime}}{k^{\prime}}p_{x}^{k^{\prime}}(1-p_{x})^{n^{\prime}-k^{\prime}}\cdot\binom{n-n^{\prime}}{k-k^{\prime}}p_{x}^{k-k^{\prime}}(1-p_{x})^{(n-n^{\prime})-(k-k^{\prime})}$
$k^{\prime}>k$ |  $0$ (_infeasible_)
$x\neq x^{\prime}$ | $k+k^{\prime}\leq n$ | $\sum_{i=\max(0,k-(n-n^{\prime}))}^{\min(k,n-k^{\prime})}\frac{n^{\prime}!}{k^{\prime}!i!(n^{\prime}-k^{\prime}-i)!}\frac{(n-n^{\prime})!}{(k-i)!((n-n^{\prime})-(k-i))!}p_{x^{\prime}}^{k^{\prime}}p_{x^{\prime}}^{k^{\prime}}(1-p_{x}-p_{x^{\prime}})^{n^{\prime}-k^{\prime}-i}p_{x}^{k}(1-p_{x})^{(n-n^{\prime})-(k-i)}$
$k+k^{\prime}>n$ |  $0$ (_infeasible_)

###### Proof.

The proof is straightforward from the definition of $N_{x}$ and
${N^{\prime}}_{x^{\prime}}$. ∎

Given the expected value and the variance of $\hat{M}_{k}^{\text{Evo}}$, the
MSE of $\hat{M}_{k}^{\text{Evo}}$ naturally follows.

##  Appendix D Details of the Hyperparameters of the Evolutionary Algorithm

For evaluating Algorithm 1, we use the following hyperparameters:

  * •

Same as the Orlitsky’s study (Orlitsky & Suresh, 2015), which assess the
performance of the Good-Turing estimator, we use the hybrid estimator
$\hat{p}$ of the empirical estimate and the Good-Turing estimate to
approximate the underlying distribution $\\{p_{x}\\}_{x\in\mathcal{X}}$ for
estimating the MSE of the evolved estimator. The hybrid estimator $\hat{p}$ is
defined as follows: If $N_{x}=k$,

$\hat{p}_{x}=\begin{cases}c\cdot\frac{k}{N}&\text{if }k<\Phi_{k+1},\\\
c\cdot\frac{\hat{M}_{k}^{G}}{\Phi_{k}}&\text{otherwise,}\end{cases}$

where $c$ is a normalization constant such that
$\sum_{x\in\mathcal{X}}\hat{p}_{x}=1$.

  * •

The number of generations $G=100$. To avoid the algorithm from converging to a
local minimum, we limit the maximum number of generations to be $2000$.

  * •

The mutant size $m=40$.

  * •

When selecting the individuals for the mutation, we use _tournament_ selection
with tournament size $t=3$, i.e., we randomly choose three individuals with
replacement and select the best one, and repeat this process $m$ times.

  * •

When choosing the top three individuals when constructing the next generation,
we use _elitist_ selection, i.e., choosing the top three individuals with the
smallest fitness values.

  * •

To avoid the estimator from being too complex, we limit the maximum number of
terms in the estimator to be $20$.

The actual script implementing Algorithm 1 can be found at the publically
available repository

https://anonymous.4open.science/r/Better-Turing-157F.

##  Appendix E Additional Experimental Results

Table 2: The MSE of the Good-Turing estimator $\hat{M}_{0}^{G}$, minimal bias estimator $\hat{M}_{0}^{B}$, and the best evolved estimator $\hat{M}_{0}^{\text{Evo}}$ and for the missing mass $M_{0}$, the success rate $\hat{A}_{12}$ of the evolved estimator ($X_{2}$) against the Good-Turing estimator ($X_{1}$), and the ratio (Ratio, $\mathit{MSE}(M_{0}^{\text{Evo}})/\mathit{MSE}(\hat{M}_{0}^{G})$) for two support sizes $S=100$ and $200$, three sample sizes $n$ and six distributions. $S$ | $n/S$ | Distribution | $\mathit{MSE}(\hat{M}_{0}^{B})$ | $\mathit{MSE}(\hat{M}_{0}^{G})$ | $\mathit{MSE}(\hat{M}_{0}^{\text{evo}})$ | $\hat{A}_{12}$ | Ratio
---|---|---|---|---|---|---|---
100 | 0.5 | uniform | 6.834e-03 | 6.681e-03 | 6.267e-03 | 62% | 93%
half&half | 6.821e-03 | 6.694e-03 | 4.489e-03 | 85% | 67%
zipf-0.5 | 6.676e-03 | 6.565e-03 | 3.311e-03 | 94% | 50%
zipf-1 | 4.995e-03 | 4.943e-03 | 3.065e-03 | 98% | 62%
diri-1 | 6.166e-03 | 6.086e-03 | 3.202e-03 | 96% | 52%
diri-0.5 | 5.223e-03 | 5.167e-03 | 2.708e-03 | 100% | 52%
1.0 | uniform | 2.365e-03 | 2.351e-03 | 1.905e-03 | 88% | 81%
half&half | 2.200e-03 | 2.190e-03 | 1.439e-03 | 97% | 65%
zipf-0.5 | 2.207e-03 | 2.194e-03 | 1.982e-03 | 75% | 90%
zipf-1 | 1.713e-03 | 1.704e-03 | 1.705e-03 | 75% | 100%
diri-1 | 1.787e-03 | 1.778e-03 | 1.066e-03 | 100% | 60%
diri-0.5 | 1.388e-03 | 1.381e-03 | 8.747e-04 | 97% | 63%
2.0 | uniform | 4.047e-04 | 4.028e-04 | 3.428e-04 | 89% | 85%
half&half | 4.237e-04 | 4.221e-04 | 3.035e-04 | 97% | 71%
zipf-0.5 | 4.580e-04 | 4.561e-04 | 3.321e-04 | 95% | 72%
zipf-1 | 4.826e-04 | 4.810e-04 | 3.633e-04 | 98% | 75%
diri-1 | 3.946e-04 | 3.932e-04 | 2.473e-04 | 100% | 62%
diri-0.5 | 3.276e-04 | 3.264e-04 | 2.587e-04 | 87% | 79%
200 | 0.5 | uniform | 3.361e-03 | 3.323e-03 | 2.044e-03 | 95% | 61%
half&half | 3.357e-03 | 3.326e-03 | 1.968e-03 | 96% | 59%
zipf-0.5 | 3.254e-03 | 3.227e-03 | 2.293e-03 | 89% | 71%
zipf-1 | 2.335e-03 | 2.324e-03 | 2.410e-03 | 74% | 103%
diri-1 | 3.011e-03 | 2.992e-03 | 2.359e-03 | 85% | 78%
diri-0.5 | 2.563e-03 | 2.550e-03 | 1.813e-03 | 94% | 71%
1.0 | uniform | 1.172e-03 | 1.169e-03 | 8.900e-04 | 99% | 76%
half&half | 1.092e-03 | 1.090e-03 | 8.584e-04 | 99% | 78%
zipf-0.5 | 1.091e-03 | 1.088e-03 | 8.525e-04 | 97% | 78%
zipf-1 | 8.185e-04 | 8.165e-04 | 7.244e-04 | 88% | 88%
diri-1 | 8.898e-04 | 8.876e-04 | 6.652e-04 | 100% | 74%
diri-0.5 | 6.900e-04 | 6.882e-04 | 4.861e-04 | 98% | 70%
2.0 | uniform | 2.017e-04 | 2.012e-04 | 1.702e-04 | 93% | 84%
half&half | 2.113e-04 | 2.109e-04 | 1.716e-04 | 100% | 81%
zipf-0.5 | 2.307e-04 | 2.302e-04 | 1.929e-04 | 100% | 83%
zipf-1 | 2.390e-04 | 2.387e-04 | 2.109e-04 | 96% | 88%
diri-1 | 1.961e-04 | 1.958e-04 | 1.648e-04 | 96% | 84%
diri-0.5 | 1.609e-04 | 1.607e-04 | 1.315e-04 | 93% | 81%


Therefore, higher-order multiples are neglected in BiPoS1, and the success of
BiPoS1 with reproducing the observed binary populations lends credence to the
initial assumptions on which BiPoS1 is based, namely the shape of the initial
period distribution function of the binaries, and the initial binarity of 100
percent.

For the treatise of higher-order multiples, Hamers et al. (2021) recently
published the Multiple Stellar Evolution code (MSE). MSE takes a single
higher-order multiple system as its input, and evolves it, including its
stellar evolution and post-Newtonian terms, besides evolution of its sub-
binaries. However, MSE does not, unlike BiPoS1, perform dynamical processing
in embedded clusters of an entire birth population of stars. It will be an
important task in the future to combine both codes to improve the population
synthesis of galactic field populations.

It is also important to study how changing the initial conditions affect the
outcome of BiPoS1. In particular, it will be very interesting to investigate
degeneracies: for example, a larger birth half-mass radius might be
compensated by a larger star formation efficiency, which leads to the same
dynamically processed binary-star population. In this case, the initial
conditions would be ´dynamically equivalent’ (Belloni et al., 2018).

Note that there are also other, earlier successful attempts to explain binary
populations with BiPoS1, although it was never properly introduced until this
paper. This concerns the papers by Marks & Kroupa (2011, 2012), Marks et al.
(2014, 2015) and Marks et al. (2017), even though adapted versions of BiPoS
were used for some of them.

There are two caveats regarding the ranges in which BiPoS1 is applied. The
first caveat is for star clusters with masses $M_{\rm ecl}>10^{6}\,{\rm
M}_{\odot}$ in the beginning. This is because GCs with masses $M\apprge
10^{6}\,{\rm M}_{\odot}$ today show evidence that their IMF was not canonical,
but top-heavy (Dabringhausen et al., 2009, 2012; Marks et al., 2012). The IMF
in BiPoS1 is however fixed to the canonical slopes. The second caveat is for
stars with masses $m\apprle 0.2\,{\rm M}_{\odot}$. In this mass range, an
additional brown-dwarf-like mass function with different binary properties
comes in besides the star-like mass function, which is the single mass
function above $m=0.2\,{\rm M}_{\odot}$ (Thies et al., 2015). However, BiPoS1
assumes that all stars belong to the star-like population (see Marks et al.
2015).

A few examples of how to operate with BiPoS1 are also given in this paper,
namely 1.) a comparison between the birth population of binaries with the
initial population after internal binary evolution, or pre-main sequence
eigenevolution, 2.) a comparison of the binary survival fractions of
progenitor models for globular clusters over the range of their embedded
stellar masses, 3.) a discussion of the high binary fraction at low apparent
separation discovered by Duchêne et al. (2018) in the Orion Nebula Cluster,
and 4.) a demonstration of the galactic field mode in BiPoS1 with two
associations of O- and B-stars.

Future developements of BiPoS might include adding brown dwarfs, hardening of
binaries, a variable IMF, and an update of the pre-main sequence
eigenevolution, or internal binary evolution (Belloni et al., 2017).

## Acknowledgements

Jörg Dabringhausen and Pavel Kroupa acknowledge support from the Grant Agency
of the Czech Republic under grant number 20-21855S.

## Data availability

All data in this paper is either available in the cited literature, or was
generated with the computer programme BiPoS1, using the commands detailed in
this Paper. BiPoS1 can be downloaded at GitHub under the web adress
https://github.com/JDabringhausen/BiPoS1.

## References

  * Aarseth (1999) Aarseth S. J., 1999, PASP, 111, 1333
  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations
  * Banerjee & Kroupa (2014) Banerjee S., Kroupa P., 2014, ApJ, 787, 158
  * Bastian et al. (2010) Bastian N., Covey K. R., Meyer M. R., 2010, ARA&A, 48, 339
  * Bate (2009) Bate M. R., 2009, MNRAS, 392, 590
  * Bate (2012) Bate M. R., 2012, MNRAS, 419, 3115
  * Beccari et al. (2017) Beccari G., Petr-Gotzens M. G., Boffin H. M. J., Romaniello M., Fedele D., Carraro G., De Marchi G., de Wit W. J., Drew J. E., Kalari V. M., Manara C. F., Martin E. L., Mieske S., Panagia N., Testi L., Vink J. S., Walsh J. R., Wright N. J., 2017, A& A, 604, A22
  * Beck et al. (2018) Beck P. G., Mathis S., Gallet F., Charbonnel C., Benbakoura M., García R. A., do Nascimento J. D., 2018, MNRAS, 479, L123
  * Belloni et al. (2017) Belloni D., Askar A., Giersz M., Kroupa P., Rocha-Pinto H. J., 2017, MNRAS, 471, 2812
  * Belloni et al. (2018) Belloni D., Kroupa P., Rocha-Pinto H. J., Giersz M., 2018, MNRAS, 474, 3740
  * Blaauw (1964) Blaauw A., 1964, ARA&A, 2, 213
  * Bonnell & Bastien (1992) Bonnell I., Bastien P., 1992, in McAlister H. A., Hartkopf W. I., eds, IAU Colloq. 135: Complementary Approaches to Double and Multiple Star Research Vol. 32 of Astronomical Society of the Pacific Conference Series, The Formation of Non-Equal Mass Binary and Multiple Systems. p. 206
  * Bressert et al. (2010) Bressert E., Bastian N., Gutermuth R., Megeath S. T., Allen L., Evans Neal J. I., Rebull L. M., Hatchell J., Johnstone D., Bourke T. L., Cieza L. A., Harvey P. M., Merin B., Ray T. P., Tothill N. F. H., 2010, MNRAS, 409, L54
  * Briceño et al. (2019) Briceño C., Calvet N., Hernández J., Vivas A. K., Mateu C., Downes J. J., Loerincs J., Pérez-Blanco A., Berlind P., Espaillat C., Allen L., Hartmann L., Mateo M., Bailey John I. I., 2019, AJ, 157, 85
  * Chabrier (2005) Chabrier G., 2005, The Initial Mass Function: From Salpeter 1955 to 2005. p. 41
  * Dabringhausen (2019) Dabringhausen J., 2019, MNRAS, 490, 848
  * Dabringhausen et al. (2010) Dabringhausen J., Fellhauer M., Kroupa P., 2010, MNRAS, 403, 1054
  * Dabringhausen et al. (2008) Dabringhausen J., Hilker M., Kroupa P., 2008, MNRAS, 386, 864
  * Dabringhausen et al. (2009) Dabringhausen J., Kroupa P., Baumgardt H., 2009, MNRAS, 394, 1529
  * Dabringhausen et al. (2012) Dabringhausen J., Kroupa P., Pflamm-Altenburg J., Mieske S., 2012, ApJ, 747, 72
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
  * Duchêne et al. (2018) Duchêne G., Lacour S., Moraux E., Goodwin S., Bouvier J., 2018, MNRAS, 478, 1825
  * Duchêne et al. (1999) Duchêne G., Monin J. L., Bouvier J., Ménard F., 1999, A& A, 351, 954
  * Duquennoy & Mayor (1991) Duquennoy A., Mayor M., 1991, A& A, 500, 337
  * Elmegreen & Efremov (1997) Elmegreen B. G., Efremov Y. N., 1997, ApJ, 480, 235
  * Evans et al. (2005) Evans N. R., Carpenter K. G., Robinson R., Kienzle F., Dekas A. E., 2005, AJ, 130, 789
  * Fischer & Marcy (1992) Fischer D. A., Marcy G. W., 1992, ApJ, 396, 178
  * Fontanot et al. (2017) Fontanot F., De Lucia G., Hirschmann M., Bruzual G., Charlot S., Zibetti S., 2017, MNRAS, 464, 3812
  * Goodman et al. (1993) Goodman A. A., Benson P. J., Fuller G. A., Myers P. C., 1993, ApJ, 406, 528
  * Goodwin & Kroupa (2005) Goodwin S. P., Kroupa P., 2005, A& A, 439, 565
  * Goodwin et al. (2007) Goodwin S. P., Kroupa P., Goodman A., Burkert A., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V The Fragmentation of Cores and the Initial Binary Population. p. 133
  * Halbwachs et al. (2003) Halbwachs J. L., Mayor M., Udry S., Arenou F., 2003, A& A, 397, 159
  * Hamers et al. (2021) Hamers A. S., Rantala A., Neunteufel P., Preece H., Vynatheya P., 2021, MNRAS, 502, 4479
  * Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
  * Hillenbrand (1997) Hillenbrand L. A., 1997, AJ, 113, 1733
  * Hillenbrand & Hartmann (1998) Hillenbrand L. A., Hartmann L. W., 1998, ApJ, 492, 540
  * Hut et al. (1992) Hut P., McMillan S., Romani R. W., 1992, ApJ, 389, 527
  * Jerabkova et al. (2019) Jerabkova T., Beccari G., Boffin H. M. J., Petr-Gotzens M. G., Manara C. F., Prada Moroni P. G., Tognelli E., Degl’Innocenti S., 2019, A& A, 627, A57
  * Jeřábková et al. (2018) Jeřábková T., Hasani Zonoozi A., Kroupa P., Beccari G., Yan Z., Vazdekis A., Zhang Z. Y., 2018, A&A, 620, A39
  * Jeřábková et al. (2017) Jeřábková T., Kroupa P., Dabringhausen J., Hilker M., Bekki K., 2017, A&A, 608, A53
  * Joncour et al. (2018) Joncour I., Duchêne G., Moraux E., Motte F., 2018, A& A, 620, A27
  * Joy (1945) Joy A. H., 1945, ApJ, 102, 168
  * Kenyon & Hartmann (1995) Kenyon S. J., Hartmann L., 1995, ApJS, 101, 117
  * Kobulnicky & Fryer (2007) Kobulnicky H. A., Fryer C. L., 2007, ApJ, 670, 747
  * Köhler et al. (2000) Köhler R., Kunkel M., Leinert C., Zinnecker H., 2000, A& A, 356, 541
  * Kohler & Leinert (1998) Kohler R., Leinert C., 1998, A& A, 331, 977
  * Krause et al. (2020) Krause M. G. H., Offner S. S. R., Charbonnel C., Gieles M., Klessen R. S., Vázquez-Semadeni E., Ballesteros-Paredes J., Girichidis P., Kruijssen J. M. D., Ward J. L., Zinnecker H., 2020, SSRv, 216, 64
  * Kroupa (1995a) Kroupa P., 1995a, MNRAS, 277, 1491
  * Kroupa (1995b) Kroupa P., 1995b, MNRAS, 277, 1522
  * Kroupa (1995c) Kroupa P., 1995c, MNRAS, 277, 1507
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kroupa (2002) Kroupa P., 2002, Binary Systems, Star Clusters and the Galactic-Field Population, Habilitation Thesis, University of Kiel Library
  * Kroupa (2008) Kroupa P., 2008, Initial Conditions for Star Clusters, LNP 760. p. 181
  * Kroupa (2011) Kroupa P., 2011, in Alves J., Elmegreen B. G., Girart J. M., Trimble V., eds, Computational Star Formation Vol. 270, The universality hypothesis: binary and stellar populations in star clusters and galaxies. pp 141–149
  * Kroupa et al. (2001) Kroupa P., Aarseth S., Hurley J., 2001, MNRAS, 321, 699
  * Kroupa & Bouvier (2003) Kroupa P., Bouvier J., 2003, MNRAS, 346, 369
  * Kroupa & Jerabkova (2018) Kroupa P., Jerabkova T., 2018, arXiv e-prints, p. arXiv:1806.10605
  * Kroupa et al. (2018) Kroupa P., Jeřábková T., Dinnbier F., Beccari G., Yan Z., 2018, A& A, 612, A74
  * Kroupa et al. (1993) Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545
  * Kroupa & Weidner (2003) Kroupa P., Weidner C., 2003, ApJ, 598, 1076
  * Kroupa et al. (2013) Kroupa P., Weidner C., Pflamm-Altenburg J., Thies I., Dabringhausen J., Marks M., Maschberger T., 2013, The Stellar and Sub-Stellar Initial Mass Function of Simple and Composite Populations. p. 115
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA& A, 41, 57
  * Lada et al. (1984) Lada C. J., Margulis M., Dearborn D., 1984, ApJ, 285, 141
  * Leigh et al. (2015) Leigh N. W. C., Giersz M., Marks M., Webb J. J., Hypki A., Heinke C. O., Kroupa P., Sills A., 2015, MNRAS, 446, 226
  * Leinert et al. (1993) Leinert C., Zinnecker H., Weitzel N., Christou J., Ridgway S. T., Jameson R., Haas M., Lenzen R., 1993, A&A, 278, 129
  * Liu (2019) Liu C., 2019, MNRAS, 490, 550
  * Marcy & Butler (2000) Marcy G. W., Butler R. P., 2000, PASP, 112, 137
  * Marks et al. (2015) Marks M., Janson M., Kroupa P., Leigh N., Thies I., 2015, MNRAS, 452, 1014
  * Marks & Kroupa (2011) Marks M., Kroupa P., 2011, MNRAS, 417, 1702
  * Marks & Kroupa (2012) Marks M., Kroupa P., 2012, A&A, 543, A8
  * Marks et al. (2012) Marks M., Kroupa P., Dabringhausen J., Pawlowski M. S., 2012, MNRAS, 422, 2246
  * Marks et al. (2011) Marks M., Kroupa P., Oh S., 2011, MNRAS, 417, 1684
  * Marks et al. (2014) Marks M., Leigh N., Giersz M., Pfalzner S., Pflamm-Altenburg J., Oh S., 2014, MNRAS, 441, 3503
  * Marks et al. (2017) Marks M., Martín E. L., Béjar V. J. S., Lodieu N., Kroupa P., Manjavacas E., Thies I., Rebolo López R., Velasco S., 2017, A& A, 605, A11
  * McCarthy & Zuckerman (2004) McCarthy C., Zuckerman B., 2004, AJ, 127, 2871
  * Megeath et al. (2016) Megeath S. T., Gutermuth R., Muzerolle J., Kryukova E., Hora J. L., Allen L. E., Flaherty K., Hartmann L., Myers P. C., Pipher J. L., Stauffer J., Young E. T., Fazio G. G., 2016, AJ, 151, 5
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Offner et al. (2014) Offner S. S. R., Clark P. C., Hennebelle P., Bastian N., Bate M. R., Hopkins P. F., Moraux E., Whitworth A. P., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI The Origin and Universality of the Stellar Initial Mass Function. p. 53
  * Oh & Kroupa (2018) Oh S., Kroupa P., 2018, MNRAS, 481, 153
  * Oh et al. (2015) Oh S., Kroupa P., Pflamm-Altenburg J., 2015, ApJ, 805, 92
  * Parker et al. (2011) Parker R. J., Goodwin S. P., Allison R. J., 2011, MNRAS, 418, 2565
  * Petr et al. (1998) Petr M. G., Coudé du Foresto V., Beckwith S. V. W., Richichi A., McCaughrean M. J., 1998, ApJ, 500, 825
  * Pflamm-Altenburg et al. (2007) Pflamm-Altenburg J., Weidner C., Kroupa P., 2007, ApJ, 671, 1550
  * Prantzos & Aubert (1995) Prantzos N., Aubert O., 1995, A&A, 302, 69
  * Press et al. (1992) Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in C. The art of scientific computing, 2nd ed.. Cambridge University Press, Cambridge
  * Raghavan et al. (2010) Raghavan D., McAlister H. A., Henry T. J., Latham D. W., Marcy G. W., Mason B. D., Gies D. R., White R. J., ten Brummelaar T. A., 2010, ApJS, 190, 1
  * Rastegaev (2010) Rastegaev D. A., 2010, AJ, 140, 2013
  * Ratzka et al. (2005) Ratzka T., Köhler R., Leinert C., 2005, A&A, 437, 611
  * Reipurth et al. (2007) Reipurth B., Guimarães M. M., Connelley M. S., Bally J., 2007, AJ, 134, 2272
  * Sana et al. (2009) Sana H., Gosset E., Evans C. J., 2009, MNRAS, 400, 1479
  * Sana et al. (2008) Sana H., Gosset E., Nazé Y., Rauw G., Linder N., 2008, MNRAS, 386, 447
  * Scally et al. (1999) Scally A., Clarke C., McCaughrean M. J., 1999, MNRAS, 306, 253
  * Schulz et al. (2015) Schulz C., Pflamm-Altenburg J., Kroupa P., 2015, A&A, 582, A93
  * Stamatellos & Whitworth (2009) Stamatellos D., Whitworth A. P., 2009, MNRAS, 392, 413
  * Thies & Kroupa (2007) Thies I., Kroupa P., 2007, ApJ, 671, 767
  * Thies et al. (2010) Thies I., Kroupa P., Goodwin S. P., Stamatellos D., Whitworth A. P., 2010, ApJ, 717, 577
  * Thies et al. (2015) Thies I., Pflamm-Altenburg J., Kroupa P., Marks M., 2015, ApJ, 800, 72
  * Tokovinin et al. (2020) Tokovinin A., Petr-Gotzens M., Briceno C., 2020, arXiv e-prints, p. arXiv:2010.08543
  * Tout (1991) Tout C. A., 1991, MNRAS, 250, 701
  * Wang et al. (2019) Wang L., Kroupa P., Jerabkova T., 2019, MNRAS, 484, 1843
  * Wang et al. (2016) Wang L., Spurzem R., Aarseth S., Giersz M., Askar A., Berczik P., Naab T., Schadow R., Kouwenhoven M. B. N., 2016, MNRAS, 458, 1450
  * Ward-Duong et al. (2015) Ward-Duong K., Patience J., De Rosa R. J., Bulger J., Rajan A., Goodwin S. P., Parker R. J., McCarthy D. W., Kulesa C., 2015, MNRAS, 449, 2618
  * Warren & Hesser (1977) Warren W. H. J., Hesser J. E., 1977, ApJS, 34, 115
  * Weidner & Kroupa (2005) Weidner C., Kroupa P., 2005, ApJ, 625, 754
  * Weidner et al. (2010) Weidner C., Kroupa P., Bonnell I. A. D., 2010, MNRAS, 401, 275
  * Weidner et al. (2004) Weidner C., Kroupa P., Larsen S. S., 2004, MNRAS, 350, 1503
  * Wielen et al. (1983) Wielen R., Jahreiß H., Krüger R., 1983, in Philip A. G. D., Upgren A. R., eds, IAU Colloq. 76: Nearby Stars and the Stellar Luminosity Function The Determination of the Luminosity Function of Nearby Stars. pp 163–170
  * Yan et al. (2020) Yan Z., Jerabkova T., Kroupa P., 2020, A& A, 637, A68
  * Zahn & Bouchet (1989) Zahn J. P., Bouchet L., 1989, A& A, 223, 112
  * Zhang et al. (2018) Zhang Z.-Y., Romano D., Ivison R. J., Papadopoulos P. P., Matteucci F., 2018, Nature, 558, 260