text:
\section{Introduction}
\section{Introduction}
\label{sec:intro}
Species estimation is an important problem in numerous scientific
disciplines.  Initially used to estimate ecological
diversity~\cite{C84,CL92,BF93,CCG12}, it was subsequently applied to
assess vocabulary size~\cite{ET76,TE87},
database attribute variation~\cite{HNSS95},
and
password innovation~\cite{FH07}.
Recently it has found a number of
bio-science applications including estimation of
bacterial and microbial diversity~\cite{KPD99,PBG01,HHJTB01,GCB07},
immune
receptor diversity~\cite{R09},
and unseen genetic variations~\cite{ICL09}.

All approaches to the problem incorporate a statistical model,
with the most popular being the \emph{extrapolation model}
introduced by Fisher, Corbet, and Williams~\cite{FCW43} in 1943.
It assumes that $n$ independent samples $X^n\stackrel{\mathrm{def}}{=} X_1,\ldots, X_n$ were collected from
an unknown distribution $p$, and calls for estimating
\[
U
\stackrel{\mathrm{def}}{=}
\UXnm
\stackrel{\mathrm{def}}{=}
\left|\sets{X_{n+1}^{n+m}}\backslash\sets{X^n}\right|,
\]
the number of hitherto unseen symbols that would be observed
if $m$ additional samples $X_{n+1}^{n+m} \stackrel{\mathrm{def}}{=} X_{n+1},\ldots, X_{n+m}$,
were collected from the same distribution.

In 1956,
Good and Toulmin~\cite{GT56} predicted $\aUXnm$ by
a fascinating estimator
that has since intrigued statisticians and a
broad range of scientists alike~\cite{K86}.
For example, in the Stanford University Statistics Department
brochure~\cite{stanfordF}, published in the early 90's
and slightly abbreviated here,
Bradley Efron credited the problem and its elegant solution with
kindling his interest in statistics.
As we shall soon see, Efron, along with Ronald Thisted,
went on to make significant contributions to this problem.

\begin{tcolorbox}
In the early 1940's, naturalist Corbet had spent two years trapping butterflies in Malaya. At the
end of that time he constructed a table (see below) to show how
many times he had trapped various butterfly species. For example,
118 species were so rare that Corbet had trapped only
one specimen of each, 74 species had been trapped twice each, etc.\\[2mm]
\centerline{
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|} \hline
Frequency & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14
& 15\\ \hline
Species  & 118& 74& 44& 24& 29& 22& 20& 19& 20& 15& 12& 14& 6& 12&
6\\ \hline
\end{tabular}
}\\[1mm]
\indent Corbet returned to England with his table, and asked R.A. Fisher, the
greatest of all statisticians, how many new species he would see if he
returned to Malaya for another two years of trapping. This question
seems impossible to answer, since it refers to a column of Corbet's
table that doesn't exist, the ``0'' column. Fisher provided an
interesting answer that was later improved on [by Good and Toulmin].
The number of new species you can expect to see in two years of
additional trapping is
\[
118-74+44-24+ \ldots -12+6 = 75.
\]
\end{tcolorbox}

This example evaluates the Good-Toulmin estimator for the special case
where the original and future samples are of equal size, namely $m=n$.
To describe the estimator's general form we need only
a modicum of nomenclature.

The \emph{prevalence} $\prev i\stackrel{\mathrm{def}}{=}\prev i(X^n)$
of an integer $i\ge0$ in $X^n$ is the number of
symbols appearing $i$ times in $X^n$.
For example, for $X^7$=\emph{bananas}, $\prev1=2$ and
$\prev2=\prev3=1$, and in Corbet's table, $\prev1=118$ and $\prev2=74$.
Let $t \stackrel{\mathrm{def}}{=} \frac{m}{n}$ be the ratio of the number of future
and past samples so that $m=tn$.
Good and Toulmin
estimated $U$
by the
surprisingly simple formula
\begin{equation}
\label{eq:GT}
\aUGTXnm
\stackrel{\mathrm{def}}{=}
\Usss{GT}{\Xn}{t}
\stackrel{\mathrm{def}}{=}
-\sum_{i=1}^\infty \Paren{-t}^i\prev{i}
\end{equation}
They showed that for all $t\leq 1$, $\aUGTxnm$ is nearly unbiased,
and that while $\aUXnm$ can be as high as $nt$,\footnote{For $a,b >0$, denote $a \lesssim b$ or $b \gtrsim a$ if $\frac{a}{b} \leq c $
 for some universal constant $c$. Denote $a \asymp b$ if both $a \lesssim b$ and $a \gtrsim b$.}
\[
\mathbb{E}
(\aUGTXnm - \aUXnm)^2 \lesssim nt^2,
\]
hence in expectation,
$\aUGTxnm$ approximates $\aUXnm$ to within just $\sqrt{n}t$.
\arxiv{Figure~\ref{fig:GT_basic}}{(Fig. 1)} shows that for the ubiquitous
Zipf distribution, $\aUGTxnm$ indeed approximates $\aUXnm$ well for
all $t<1$.
\begin{figure}[ht]
\centering
\arxiv{\subfigure[]{\includegraphics[scale=0.5]{FGT_zipf-crop.pdf}}}{\subfigure[]{\includegraphics[scale=0.5]{FGT_zipf-crop.pdf}}}
\arxiv{\subfigure[]{\includegraphics[scale=0.5]{FGT_zipf_2-crop.pdf}}}{\subfigure[]{\includegraphics[scale=0.5]{FGT_zipf_2-crop.pdf}}}
\caption{GT estimate as a function of $t$ for two realizations random samples of
size $n=5000$ generated by a Zipf distribution $p_i \propto 1/(i+10)$
for $1\le i\le 10000$.}
\label{fig:GT_basic}\end{figure}
Naturally, we would like to estimate $\aUXnm$ for as large a $t$ as
possible.
However, as $t>1$ increases, $\aUGTxnm$ grows as
$(-t)^i\prev{i}$ for
the largest $i$ such that $\prev i>0$.  Hence whenever any symbol
appears more than once, $\aUGTxnm$ grows super-linearly in $t$,
eventually far exceeding $\aUXnm$ that grows at most linearly in
$t$.
\arxiv{Figure~\ref{fig:GT_basic}}{(Fig. 1)} also shows that for the same Zipf distribution,
for $t>1$ indeed  $\aUGTxnm$ does not approximate $\aUXnm$ at all.

To predict $\aUXnm$ for $t>1$, Good and Toulmin~\cite{GT56}
suggested using the \emph{Euler transform}~\cite{AS64}
that converts an alternating series into another series with
the same sum, and heuristically often converges faster.
Interestingly, Efron and Thisted~\cite{ET76} showed that when the
Euler transform of $\UGT$ is truncated after $k$ terms,
it can be expressed as another simple linear estimator,
\begin{equation*}
U^{\scriptscriptstyle\rm ET}
\stackrel{\mathrm{def}}{=}
\sum_{i=1}^n \coeffs{\rm ET}_i\cdot \prev{i}
\label{eq:UET}
\end{equation*}
where
\begin{equation*}
\label{eq:coeffET}
\coeffs{\rm ET}_i
\stackrel{\mathrm{def}}{=}
-\Paren{-t}^{i}\cdot\prob{\mathrm{Bin}\Big(k, \frac{1}{1+t}\Big) \geq i},
\end{equation*}
and
\[
\prob{\mathrm{Bin}\Paren{k,\frac{1}{1+t}} \geq i}
=
\begin{cases}
\sum_{j=i}^{k}
\binom{k}{j}\frac{t^{k-j}}{(1+t)^{k}}
& i \leq k,\\
0 & i > k,
\end{cases}
\]
is the \emph{binomial} tail probability that decays with $i$,
thereby moderating the rapid growth of $(-t)^i$.

Over the years, $U^{\scriptscriptstyle\rm ET}$ has been used by numerous researchers in a variety
of scenarios and a multitude of applications.
Yet despite its wide-spread use and robust empirical results,
no provable guarantees have been established for its performance
or that of any related estimator when $t>1$.
The lack of theoretical understanding, has also precluded
clear guidelines for choosing the parameter $k$ in $U^{\scriptscriptstyle\rm ET}$.

	\label{sec:new}
\section{Approach and results}
We construct a family of estimators that
\emph{provably} predict $\aUXnm$ optimally not just for constant $t>1$,
but all the way up to $t \propto \log n$.
This shows that per each observed sample, we can infer properties
of $\log n$ yet unseen samples.
The proof technique is general and provides a disciplined guideline
for choosing  the parameter $k$ for $U^{\scriptscriptstyle\rm ET}$ and, in addition, a modification
that outperforms $U^{\scriptscriptstyle\rm ET}$.

\subsection{Smoothed Good-Toulmin (SGT) estimator}
	\label{sec:SGT}

To obtain a new class of estimators, we too start with $\UGT$, but
unlike $U^{\scriptscriptstyle\rm ET}$ that was derived from $\UGT$ via analytical
considerations aimed at improving the convergence rate,
we take a probabilistic view that controls the bias and
variance of $\UGT$ and balances the two to obtain a more
efficient estimator.

Note that what renders $\UGT$
inaccurate when $t > 1$ is not its bias but mainly its high variance
due to the exponential growth of the coefficients $(-t)^i$ in \prettyref{eq:GT}; in fact $\UGT$
is the unique unbiased estimator for all $t$ and $n$ in the closely related Poisson
sampling model (see \prettyref{sec:def}).
Therefore it is tempting to truncate the series \prettyref{eq:GT} at the
$\ell^\mathrm{th}$ term and use the partial sum as an estimator:
\begin{equation}
\label{eq:UGT-ell}
U^{\ell}
\stackrel{\mathrm{def}}{=}
-\sum_{i=1}^\ell \Paren{-t}^i\prev{i}
\end{equation}
However, for $t>1$, it can be shown that for certain distributions most of the symbols typically appear
$\ell$ times and hence the last term in \prettyref{eq:UGT-ell} dominates, resulting in a large bias
and inaccurate estimates regardless of the choice of $\ell$ (see \prettyref{sec:taylor} for a rigorous justification).

To resolve this problem, we truncate the Good-Toulmin estimator at a
random location, denoted by an independent random nonnegative integer $L$, and
average over the distribution of $L$, which yields the following estimator:
\begin{equation}
U^{\scriptscriptstyle\rm L} = \mathbb{E}_L \left[- \sum_{i=1}^L
  \Paren{-t}^i\prev{i}\ignore{(x^n)}\right].
	\label{eq:UL0}
\end{equation}
 The key insight is that since the bias of $U^\ell$
typically alternates signs as $\ell$ grows, averaging over
different cutoff locations takes advantage of the cancellation and
dramatically reduces the bias.
Furthermore, the estimator \prettyref{eq:UL0} can be expressed simply as a linear combination of prevalences:
\begin{equation}
U^{\scriptscriptstyle\rm L} =  \mathbb{E}_L \left[-
  \sum_{i\geq 1} \Paren{-t}^i\prev{i} \mathbbm{1}_{i \leq L} \ignore{(x^n)}
  \right] = - \sum_{i \geq 1} \Paren{-t}^i \prob{L \ge i} \prev{i}.
	\label{eq:UL}
\end{equation}
We shall refer to estimators of the form \prettyref{eq:UL} \emph{Smoothed Good-Toulmin (SGT)} estimators and the distribution of $L$ the \emph{smoothing distribution}.

Choosing different smoothing distributions results a variety of
linear estimators, where the tail probability $\prob{L \ge i}$ compensates the exponential growth of $\Paren{-t}^i$ thereby stabilizing the variance.  Surprisingly, though the motivation and approach are quite different,
SGT estimators include $U^{\scriptscriptstyle\rm ET}$ in \prettyref{eq:UET} as a special case which corresponds to
the binomial smoothing $L\sim\mathrm{Bin}(k,  \frac{1}{1+t})$.
This provides an intuitive probabilistic interpretation of $U^{\scriptscriptstyle\rm ET}$, which was
originally derived via Euler's transform and analytic considerations.
As we show in the next section, this interpretation leads to the first theoretical
guarantee for $U^{\scriptscriptstyle\rm ET}$ as well as improved estimators that are provably optimal.

\subsection{Main results}
\label{sec:main}

Since $\aUXnm$ takes in values between $0$ and $nt$,
 we measure the performance of an estimator $U^{\scriptscriptstyle\rm E}$ by the worst-case \emph{normalized mean-square error (NMSE)},
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm E})
\stackrel{\mathrm{def}}{=}
\max_p \mathbb{E}_{p}\Paren{\frac{U^{\scriptscriptstyle\rm E} - \aUXnm}{nt}}^2.
\]
Observe that this criterion conservatively evaluates the performance of the
estimator for the worst possible
distribution.
The trivial estimator that always predicts $nt/2$ new elements has NMSE equal to $1/4$,
and we would like to construct estimators with vanishing NMSE,
which can estimate $\aUXnm$ up to an error that diminishes with $n$,
regardless of the data-generating distribution; in particular, we are interested in the largest $t$ for which this is possible.

Relating the bias and variance of $U^{\scriptscriptstyle\rm L}$ to the expectation of $t^L$
and another functional we obtain the following performance guarantee for SGT estimators with appropriately chosen smoothing distributions.
\begin{Theorem}
\label{thm:res_main}
For Poisson or binomially distributed $L$ with the parameters given in
Table\arxiv{~\ref{tab:values}}{~1}, for all $t \geq 1$ and $n\in \mathbb{N}$,
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm L})
\lesssim
\frac{1}{n^{1/t}}.
\]
\end{Theorem}
\begin{table}[ht]
\begin{center}
\begin{tabular}{|c|l|l|}
\hline
Smoothing distribution &\hfill Parameters\hfill${}$ &\hfill  $\cE_{n,t}(U^{\scriptscriptstyle\rm L})\lesssim$\hfill${}$\\ \hline
Poisson $(r)$ & $r=\frac{1}{2t}\log_e \frac{n(t+1)^2}{t-1}$ & $n^{-1/t}$\\[0.4cm]
Binomial $(k, q)$ & $k =  \left \lceil \frac{1}{2} \log_2  \frac{nt^2}{t-1} \right \rceil $, $q =\frac{1}{t+1} $ & $n^{-\log_2(1+1/t)}$  \\[0.4cm]
Binomial $(k, q)$ & $k =  \left \lceil \frac{1}{2} \log_3  \frac{nt^2}{t-1} \right \rceil $, $q = \frac{2}{t+2}$ & $n^{-\log_3(1+2/t)}$ \\[0.4cm]   \hline
\end{tabular}
\caption{NMSE of SGT estimators for three smoothing distributions.
Since for any $t \geq 1$, $\log_3(1+2/t) \geq \log_2(1+1/t) \geq 1/t$,
binomial smoothing with $q = 2/(2+t)$ yields the best convergence rate.}
\label{tab:values}
\end{center}
\end{table}

\prettyref{thm:res_main} provides a principled way for choosing the
parameter $k$ for $U^{\scriptscriptstyle\rm ET}$ and the first provable guarantee for
its performance, shown in \prettyref{tab:values}.
Furthermore, the result shows that a modification of $U^{\scriptscriptstyle\rm ET}$
with $q= \frac{2}{t+2}$ enjoys even faster convergence rate
and, as experimentally demonstrated in \prettyref{sec:experiments}, outperforms
the original version of Efron-Thisted as well as other state-of-the-art estimators.

Furthermore,
SGT estimators are essentially optimal as witnessed by the following matching minimax lower bound.
\begin{Theorem}
\label{thm:res_lb}
There exist universal constant $c,c'$ such that for any $t\geq c$, any $n\in \mathbb{N}$,  and any estimator $U^{\scriptscriptstyle\rm E}$
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm E}) \gtrsim \frac{1}{n^{c'/t}}.
\]
\end{Theorem}

Theorems \ref{thm:res_main} and \ref{thm:res_lb}
determine the limit of predictability up to a
constant multiple.
\begin{Corollary}
For any $\delta > 0 $,
\[
\lim_{n \to \infty} \frac{\max\sets{t: \cE_{n,t}(U^{\scriptscriptstyle\rm E})<\delta\textrm{\,\, for some $U^{\scriptscriptstyle\rm E}$}}}{\log n}
\asymp
\frac{1}{\log \frac{1}{\delta}}
\]
\end{Corollary}

The rest of the paper is organized as follows: In
Section~\ref{sec:def}, we describe the four statistical models
commonly used across various scientific disciplines, namely, the
multinomial, Poisson, hypergeometric, and Bernoulli product models.
Among the four models Poisson is the simplest to analyze and hence in
Sections~\ref{sec:prelim} and~\ref{sec:upper}, we first prove \prettyref{thm:res_main}
for the Poisson model and in Section~\ref{sec:extensions} we prove
similar results for the other three statistical models. In
Section~\ref{sec:lower}, we prove the lower bound for the multinomial and Poisson
models. Finally, in Section~\ref{sec:experiments} we demonstrate the
efficiency and practicality of our estimators on a variety of
synthetic and data sets.
\section{Statistical models}
\label{sec:def}
The extrapolation paradigm has been applied to several statistical
models.
In all of them, an initial sample of size related to $n$
is collected, resulting in a set $S_{\text{old}}$ of observed elements.
We consider collecting a new sample of size related to $m$,
that would result in a yet unknown set $S_{\text{new}}$ of observed elements,
and we would like to estimate
\[
|S_{\text{new}}\backslashS_{\text{old}}|,
\]
the number of unseen symbols that will appear in the new sample.
For example, for the observed sample \texttt{bananas} and future
sample \texttt{sonatas},
$S_{\text{old}}=\sets{\tt a,b,n,s}$,
$S_{\text{new}}=\sets{\tt a,n,o,s,t}$, and $|S_{\text{new}}\backslashS_{\text{old}}|=|\sets{\tt o,t}|=2$.

Four statistical models have been commonly used in the
literature (cf.~survey \cite{BF93} and \cite{CCG12}), and our
results apply to all of them. The first three statistical models are
also referred as the~\emph{abundance models} and the last one is often
referred to as the~\emph{incidence model} in ecology~\cite{CCG12}.
\begin{description}
\item[Multinomial:]
This is Good and Toulmin's original model where the samples are independently and identically distributed (i.i.d.),
and the initial and new samples consist of exactly $n$ and $m$
elements respectively. Formally,
$X^{n+m}=X_1,\ldots, X_{n+m}$ are generated independently according
to an unknown discrete distribution of finite or even infinite support,
$S_{\text{old}}=\sets{X^n}$, and $S_{\text{new}}=\sets{X_{n+1}^{n+m}}$.

\item[Hypergeometric:]
This model corresponds to a sampling-without-replacement variant of the multinomial model.
Specifically, $X^{n+m}$ are drawn uniformly without replacement from an unknown
collection of symbols that may contain repetitions, for example, an urn
with some white and black balls. Again, $S_{\text{old}}=\sets{X^n}$
and $S_{\text{new}}=\sets{X_{n+1}^{n+m}}$.

\item[Poisson:]
As in the multinomial model, the samples are also \emph{i.i.d.\xspace},
but the sample sizes, instead of
being fixed, are Poisson distributed.
Formally, $N\sim\mathrm{poi}(n)$, $M\sim\mathrm{poi}(m)$, $X^{N+M}$ are generated
independently according to an unknown discrete distribution,
$S_{\text{old}}=\sets{X^{N}}$, and $S_{\text{new}}=\sets{X_{N+1}^{N+M}}$.

\item[Bernoulli-product:]
In this model we observe signals
from a collection of independent processes over subset of an unknown set
$\mathcal{X}$. Every $x\in\mathcal{X}$ is associated with an unknown probability $0\le
p_x\le 1$, where the probabilities do not necessarily sum to 1.  Each
sample $X_i$ is a \emph{subset} of $\mathcal{X}$ where symbol $x\in\mathcal{X}$
appears with probability $p_x$ and is absent with probability $1-p_x$,
independently of all other symbols. $S_{\text{old}}=\cup_{i=1}^n X_{i}$ and
$S_{\text{new}}=\cup_{i=n+1}^{n+m}X_i$.
\end{description}

For theoretical analysis in Sections \ref{sec:prelim} and \ref{sec:upper} we use the Poisson sampling model as the leading example
 due to its
simplicity. Later in \prettyref{sec:extensions}, we show that very similar
results continue to hold for the other three models.

We close this section by discussing two problems that are closely related to the
extrapolation model, namely,
\emph{support size estimation} and \emph{missing mass estimation}, which correspond to $m=\infty$ and $m=1$ respectively.
Indeed, the probability that the next sample is new is precisely the expected value of $U$ for $m=1$,
which is the goal in the basic Good-Turing problem~\cite{G53,R68, MS00,OS15}.
On the other hand, any estimator $\aUEXnm$ for $U$ can be converted to a (not necessarily good) support size
estimator by adding the number of observed symbols.
Estimating the support size of an underlying distribution has been
studied by both ecologists~\cite{C84, CL92, BF93} and
theoreticians~\cite{RRSS09, VV11, VV13, WY14};
\ignore{ In ecology, several lower
  bounds on the underlying support size are known.  Most ecology
  papers study the support estimation in this context~\cite{C84, CL92,
    BF93}. While no upper bounds on the support can be
  rigorously proven, provable lower bounds can be obtained~\cite{}.}%
  however, to make the problem non-trivial,
all statistical models impose a lower bound on the minimum non-zero
probability of each symbol, which is assumed to be known to the statistician.  We
discuss these estimators and their differences to our results in Section~\ref{sec:support}.

\section{Preliminaries and the Poisson model}
\label{sec:prelim}
Throughout the paper,
we use standard asymptotic notation,
e.g., for any positive sequences $\{a_n\}$ and $\{b_n\}$, denote $a_n=\Theta(b_n)$ or $a_n  \asymp b_n$
if $1/c\le a_n/ b_n \le c$ for some universal constant $c>0$.
Let $\mathbbm{1}_{A}$ denote the indicator random variable of an event $A$.
Let $\mathrm{Bin}(n,p)$ denote the binomial distribution with $n$ trials and success probability $p$ and let $\mathrm{poi}(\lambda)$ denote the Poisson distribution with mean $\lambda$.
All logarithms are with respect to the natural base unless otherwise specified.

Let $p$ be a probability distribution over a discrete set $\mathcal{X}$, namely $p_x \geq 0$ for all $x\in \mathcal{X}$
 and $\sum_{x \in \mathcal{X}} p_x = 1$.
Recall that the sample sizes are Poisson distributed: $N \sim \mathrm{poi}(n)$, $M \sim \mathrm{poi}(m)$, and  $
t = \frac{m}{n}$. We abbreviate the number of unseen symbols by
\[
U \stackrel{\mathrm{def}}{=} \UXNM,
\]
and we denote an estimator by $U^{\scriptscriptstyle\rm E}\stackrel{\mathrm{def}}{=}\UEXNm$.

Let $\mul{x}$ and $\mul{x}'$ denote the multiplicity of a symbol $x$
in the current samples and future samples, respectively.
Let $\lambda_x\stackrel{\mathrm{def}}{=} np_x$. Then a symbol $x$
appears $\mul{x}\sim\mathrm{poi}(np_x)=\mathrm{poi}(\lambda_x)$ times, and for any $i\geq 0$,
\[
\mathbb{E}[\mathbbm{1}_{\mul{x} = i}] = e^{-\lambda_x} \frac{\lambda_x^i}{i!}.
\]
Hence
\[
\mathbb{E}[\prev{i}] = \mathbb{E}\left[\sum_{x} \mathbbm{1}_{\mul{x} = i} \right] =
\sum_{x} e^{-\lambda_x} \frac{\lambda_x^i}{i!}.
\]
A helpful property of Poisson sampling is that the multiplicities of
different symbols are independent of each other.
Therefore, for any function $f(x,i)$,
\[
\mathrm{Var}\left(\sum_x f(x,N_x)\right) = \sum_x \mathrm{Var}(f(x,N_x)).
\]
Many of our derivations rely on these three equations.
For example,
\begin{align*}
\mathbb{E}[U]
 = \sum_{x} \mathbb{E}[\mathbbm{1}_{\mul{x} = 0}] \cdot \mathbb{E}[\mathbbm{1}_{\mul{x}' > 0}]
= \sum_{x} e^{-\lambda_x}\cdot (1-e^{-t\lambda_x}),
\end{align*}
and
\begin{align*}
\mathrm{Var}(U)
&=
\mathrm{Var} \Paren{\sum_x \mathbbm{1}_{\mul{x} = 0} \cdot \mathbbm{1}_{\mul{x}' > 0}}
=
\sum_x \mathrm{Var} \Paren{\mathbbm{1}_{\mul{x} = 0} \cdot \mathbbm{1}_{\mul{x}' > 0}}  \\
&
\le
\sum_x \mathbb{E} \left[\mathbbm{1}_{\mul{x} = 0}  \cdot  \mathbbm{1}_{\mul{x}' > 0}\right]
=
\mathbb{E}\left[ U\right].
\end{align*}
Note that these equations imply that the standard deviation of $U$ is
at most $\sqrt{\mathbb{E}[U]} \ll \mathbb{E}[U]$, hence $U$ highly
concentrates around its expectation, and estimating $U$ and $\mathbb{E}[U]$ are
essentially the same.

\subsection{The Good-Toulmin estimator}
Before proceeding with general estimators, we prove a few properties of $\UGT$.  Under the Poisson model, $\UGT$ is in fact the \emph{unique} unbiased estimator for $U$.
\begin{Lemma}[\cite{ET76}]
For any distribution,
\[
\mathbb{E}[U] = \mathbb{E}[\UGT].
\]
\end{Lemma}
\begin{proof}
\begin{align*}
\mathbb{E} [U]
&=
\mathbb{E} \left[\sum_{x} \mathbbm{1}_{\mul{x} = 0} \cdot \mathbbm{1}_{\mul{x} > 0}\right]
=
\sum_{x} e^{-\lambda_x} \cdot  \left(1 - e^{-t\lambda_x} \right) \\
&=
-\sum_x e^{-\lambda_x} \cdot \sum^\infty_{i=1} \frac{(-t\lambda_x)^i}{i!}
=
-\sum^\infty_{i=1} (-t)^i\cdot \sum_x e^{-\lambda_x} \frac{\lambda_x^i}{i!}  \\
& =
-\sum^\infty_{i=1} (-t)^i\cdot \mathbb{E}[ \prev{i}] = \mathbb{E}[\UGT]. \qedhere
\end{align*}
\end{proof}
Even though $\UGT$ is unbiased for all $t$, for $t>1$ it has
high variance and hence does not estimate $U$ well even for the simplest distributions
\begin{Lemma}
\label{lem:FGT}
For any $t >1$,
\begin{align*}
\lim_{n \to \infty} \cE_{n,t}(U^{\scriptscriptstyle \rm GT}) = \infty.
\end{align*}
\end{Lemma}
\begin{proof}
Let $p$ be the uniform distribution over two symbols $a$ and $b$, namely,
$p_a = p_b = 1/2$. First consider even $n$.
Since $(\UGT - U)^2$ is always nonnegative,
\[
\mathbb{E}[(\UGT-U)^2] \geq \Pr(\mul{a} =\mul{b} = n/2) (2(-t)^{n/2})^2 =
\left(e^{-n/2} \frac{(n/2)^{n/2}}{(n/2)!} \right)^2 4 t^n \geq \frac{4t^n}{e^2n},
\]
where we used the fact that $k! \leq (\frac{k}{e})^k \sqrt{k} e$.
Hence for $t > 1$,
\[
\lim_{n \to \infty} \frac{\mathbb{E}[(\UGT-U)^2] }{(nt)^2}
\geq \lim_{n \to \infty}  \frac{4t^n}{e^2n (nt)^2} = \infty.
\]
The case of odd $n$ can be shown similarly by considering the event $\mul{a} = \lfloor n/2\rfloor, \mul{b} = \lceil n/2\rceil$.
\ignore{
In this case, according to Eq.\prettyref{eq:Upoisson}
we have $U = \prev{0} (1-e^{-nt/2})$ and hence
\begin{align*}
\UGT - U
= -\sum^{\infty}_{i=1} \prev{i} (-t)^{i} - \prev{0} (1-e^{-nt/2})
 = \sum_{x \in \{a,b\}} \left(-\sum^{\infty}_{i=1} \mathbbm{1}_{\mul{x} = i} \cdot (-t)^{i} -  \mathbbm{1}_{\mul{x} = 0}\cdot (1-e^{-nt/2})\right).
\end{align*}
Since the GT estimator is unbiased, we have
\begin{align*}
\mathbb{E}[(\UGT - U )^2]
& = \mathrm{Var}(\UGT - U) \\
& = \mathrm{Var} \left( \sum_{x \in \{a,b\}} -\left(\sum^{\infty}_{i=1} \mathbbm{1}_{\mul{x} = i} (-t)^{i} -  \mathbbm{1}_{\mul{x} = 0} (1-e^{-nt/2}) \right)\right) \\
& =  \sum_{x \in \{a,b\}} \mathrm{Var} \left(-\sum^{\infty}_{i=1} \mathbbm{1}_{\mul{x} = i} (-t)^{i} -  \mathbbm{1}_{\mul{x} = 0} (1-e^{-nt/2}) \right) \\
& \stackrel{(a)}{=} 2 \mathrm{Var} \left(-\sum^{\infty}_{i=1} \mathbbm{1}_{\mul{a} = i} (-t)^{i} -  \mathbbm{1}_{\mul{a} = 0} (1-e^{-nt/2}) \right) \\
& = 2 \mathbb{E} \left[ \left(-\sum^{\infty}_{i=1} \mathbbm{1}_{\mul{a} = i} (-t)^{i} -  \mathbbm{1}_{\mul{a} = 0} (1-e^{-nt/2}) \right)^2\right] \\
& \stackrel{(b)}{=} 2 \mathbb{E} \left[\sum^{\infty}_{i=1} \mathbbm{1}_{\mul{a} = i} (-t)^{2i} +  \mathbbm{1}_{\mul{a} = 0} (1-e^{-nt/2})^2 \right] \\
& = 2 \sum^\infty_{i=1} e^{-n/2} \frac{(n/2)^i t^{2i}}{i!} + 2e^{-n/2}(1-e^{-nt/2}) \\
& = 2 e^{\frac{n(t^2-1)}{2}} - 2e^{-n(t+1)/2},
\end{align*}
where (a) follows from symmetry as both symbols have identical probability, and
(b) follows from the fact that for any indicator random variable $\mathbbm{1}^2 =  \mathbbm{1}$. Hence for any $t > 1$
\[
\lim_{n \to \infty} \frac{\mathbb{E}[(\UGT-U)^2] }{(nt)^2} = \infty.
 \qedhere
\]
}%
\end{proof}
\subsection{General linear estimators}
Following \cite{ET76}, we consider general linear estimators of the form
\begin{equation}
U^{\scriptscriptstyle\rm h} = \sum^\infty_{i=1} \prev{i} \cdot h_i,
	\label{eq:Uh}
\end{equation}
which can be identified with
a formal power series $h(y)=
\sum^\infty_{i=1}\frac{h_i y^i}{i!}$. For example, $\UGT$ in \prettyref{eq:GT}
corresponds to the function $h(y) = 1-e^{-yt}$.
The next lemma bounds the bias and variance of any linear estimator $U^{\scriptscriptstyle\rm h}$ using
properties of the function $h$. In \prettyref{sec:smooth} we apply this result to
the SGT estimator
whose coefficients are of the specific form:
\begin{equation*}
h_i = -\Paren{-t}^{i}\cdot\prob{L \geq i}.
	\label{eq:hL}
\end{equation*}
Let $\prev{+} \stackrel{\mathrm{def}}{=} \sum^\infty_{i=1} \prev{i}$ denote the number of observed symbols.
\begin{Lemma}
\label{lem:general_bounds}
The bias of $U^{\scriptscriptstyle\rm h}$ is
\begin{equation*}
\label{eq:Upoisson}
\mathbb{E}[U^{\scriptscriptstyle\rm h} - U] = \sum_{x} e^{-\lambda_x} \left( h(\lambda_x) - (1-e^{-t\lambda_x}) \right),
\end{equation*}
and the variance satisfies
\[
\mathrm{Var}(U^{\scriptscriptstyle\rm h} - U) \leq \mathbb{E}[\prev{+}]  \cdot \sup_{i \geq 1} h^2_i + \mathbb{E}[U].
\]
\end{Lemma}
\begin{proof}
Note that
\begin{align}
U^{\scriptscriptstyle\rm h} - U
& = \sum^\infty_{i=1} \prev{i} h_i -  \sum_x \mathbbm{1}_{\mul{x} = 0} \cdot \mathbbm{1}_{\mul{x}' > 0}\nonumber \\
& = \sum^\infty_{i=1}
  \sum_{x} \mathbbm{1}_{\mul{x}=i} \cdot h_i - \sum_x \mathbbm{1}_{\mul{x} = 0}
  \cdot  \mathbbm{1}_{\mul{x}' > 0}
\nonumber \\
& = \sum_x \left(\sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i -
  \mathbbm{1}_{\mul{x} = 0} \cdot  \mathbbm{1}_{\mul{x}' > 0} \right).
\nonumber
\end{align}
For every symbol $x$,
\begin{align}
\mathbb{E} \left[\sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i -
  \mathbbm{1}_{\mul{x} = 0} \cdot  \mathbbm{1}_{\mul{x}' > 0} \right]
& =
 \sum^\infty_{i=1} e^{-\lambda_x}
\frac{\lambda^i_x}{i!} \cdot h_i - e^{-\lambda_x} \cdot (1 -
e^{-t\lambda_x})  \nonumber \\
& = e^{-\lambda_x} \left(
\sum^\infty_{i=1} \frac{\lambda^i_x h_i}{i!} - (1-e^{-t\lambda_x})
\right) \nonumber \\
& =  e^{-\lambda_x} \left( h(\lambda_x) -
(1-e^{-t\lambda_x}) \right),
\nonumber
\end{align}
from which \prettyref{eq:Upoisson} follows.
\ignore{
The bias of such an estimator is
\begin{align*}
\mathbb{E}[U^{\scriptscriptstyle\rm h} - U] & = \mathbb{E}\left[\sum^\infty_{i=1} \prev{i} h_i
  \right] - \mathbb{E}\left[ \sum_x \mathbbm{1}_{\mul{x} = 0} \cdot \left(1 -
  e^{-t\lambda_x} \right)\right] \\
& = \mathbb{E}\left[\sum^\infty_{i=1}
  \sum_{x} \mathbbm{1}_{\mul{x}=i} \cdot h_i - \sum_x \mathbbm{1}_{\mul{x} = 0}
  \cdot \left(1 - e^{-t\lambda_x} \right)\right] \\
& = \sum_{x}
\mathbb{E}\left[\sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i -
  \mathbbm{1}_{\mul{x} = 0} \cdot \left(1 - e^{-t\lambda_x} \right)\right]
\\ & = \sum_{x} \left( \sum^\infty_{i=1} e^{-\lambda_x}
\frac{\lambda^i_x}{i!} \cdot h_i - e^{-\lambda_x} \cdot \left(1 -
e^{-t\lambda_x} \right)\right) \\
& = \sum_{x} e^{-\lambda_x} \left(
\sum^\infty_{i=1} \frac{\lambda^i_x h_i}{i!} - (1-e^{-t\lambda_x})
\right) \\
& = \sum_{x} e^{-\lambda_x} \left( h(\lambda_x) -
(1-e^{-t\lambda_x}) \right).
\end{align*}
}%
\ignore{
For every symbol $x$,
\begin{align*}
\mathrm{Var} \left(
\sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i - \mathbbm{1}_{\mul{x} = 0}
\cdot \left(1 - e^{-t\lambda_x} \right) \right)
& \leq \mathbb{E}
\left[ \left( \sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i -
  \mathbbm{1}_{\mul{x} = 0} \cdot \left(1 - e^{-t\lambda_x} \right)
  \right)^2 \right] \\
& \stackrel{(a)}{=}  \mathbb{E} \left[ \sum^\infty_{i=1}
  \mathbbm{1}_{\mul{x}=i}h_i^2 \right] + \mathbb{E}[\mathbbm{1}_{\mul{x} = 0}] \cdot
\left(1 - e^{-t\lambda_x} \right)^2 \\
& \leq \sum^\infty_{i=1}  \mathbb{E}
[  \mathbbm{1}_{\mul{x}=i}]h_i^2 +
\mathbb{E}[\mathbbm{1}_{\mul{x} = 0}] \cdot \left(1 - e^{-t\lambda_x} \right).
\end{align*}
For every $i \neq j$, $\mathbb{E}[ \mathbbm{1}_{\mul{x}=i}  \mathbbm{1}_{\mul{x}=j}] = 0$ and hence $(a)$.
}%
For the variance, observe that for every symbol $x$,
\begin{align*}
\mathrm{Var} \left(
\sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i - \mathbbm{1}_{\mul{x} = 0}
\cdot \mathbbm{1}_{\mul{x}'>0} \right)
& \leq \mathbb{E}
\left[ \left( \sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i -
  \mathbbm{1}_{\mul{x} = 0} \cdot \mathbbm{1}_{\mul{x}' >0}
  \right)^2 \right] \\
& \stackrel{(a)}{=}  \mathbb{E} \left[ \sum^\infty_{i=1}
  \mathbbm{1}_{\mul{x}=i}h_i^2 \right] + \mathbb{E}[\mathbbm{1}_{\mul{x} = 0}] \cdot
 \mathbb{E}[\mathbbm{1}_{\mul{x}' >0}] \\
& = \sum^\infty_{i=1}  \mathbb{E}
[  \mathbbm{1}_{\mul{x}=i}]\cdot h_i^2 +
\mathbb{E}[\mathbbm{1}_{\mul{x} = 0}] \cdot  \mathbb{E}[\mathbbm{1}_{\mul{x}' >0}],
\end{align*}
where $(a)$ follows as for every $i \neq j$, $\mathbb{E}[ \mathbbm{1}_{\mul{x}=i}  \mathbbm{1}_{\mul{x}=j}] = 0$.
Since the variance of a sum of independent random variables is the sum of variances,
\begin{align*}
\mathrm{Var}(U^{\scriptscriptstyle\rm h} - U)
& \leq  \sum_{x}
\sum^\infty_{i=1}\mathbb{E}[  \mathbbm{1}_{\mul{x}=i} ]h_i^2 +
\sum_x \mathbb{E}[\mathbbm{1}_{\mul{x}=0}] \cdot  \mathbb{E}[\mathbbm{1}_{\mul{x}' >0}]\\
& =
\sum^\infty_{i=1} \mathbb{E}[\prev{i}] \cdot h_i^2 + \mathbb{E}[U] \\
& \leq
\mathbb{E}[\prev{+}] \cdot \sup_{i \geq 1} h^2_i + \mathbb{E}[U].
\qedhere
\end{align*}
\ignore{\begin{align*}
\mathrm{Var}(U^{\scriptscriptstyle\rm h} - U) & = \mathrm{Var} \left( \sum^\infty_{i=1} \sum_{x}
\mathbbm{1}_{\mul{x}=i} \cdot h_i - \sum_x \mathbbm{1}_{\mul{x} = 0} \cdot
(1 - e^{-t\lambda_x})\right) \\
 & = \sum_x \mathrm{Var} \left(
\sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i - \mathbbm{1}_{\mul{x} = 0}
\cdot \left(1 - e^{-t\lambda_x} \right) \right) \\
& \leq \sum_x \mathbb{E}
\left[ \left( \sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \cdot h_i -
  \mathbbm{1}_{\mul{x} = 0} \cdot \left(1 - e^{-t\lambda_x} \right)
  \right)^2 \right] \\
& = \sum_x \left( \mathbb{E} \left[ \sum^\infty_{i=1}
  \mathbbm{1}_{\mul{x}=i} \right]h_i^2 + \mathbb{E}[\mathbbm{1}_{\mul{x} = 0}] \cdot
\left(1 - e^{-t\lambda_x} \right)^2 \right)\\
& \leq \sum_x \left( \mathbb{E}
\left[ \sum^\infty_{i=1} \mathbbm{1}_{\mul{x}=i} \right]h_i^2 +
\mathbb{E}[\mathbbm{1}_{\mul{x} = 0}] \cdot t\lambda_x \right)\\
& =
\sum^\infty_{i=1} \mathbb{E}[\prev{i}] \cdot h_i^2 + \mathbb{E}[\prev{1}] t \\
& =
&\leq
\mathbb{E}[\prev{+}] \cdot \left(t + \sup_{i \geq 1} \frac{h^2_i}{i} \right).
\qedhere
\end{align*}
}
\end{proof}
\prettyref{lem:general_bounds} enables us to reduce the estimation
problem to a task on approximating functions. Specifically, in view of \eqref{eq:Upoisson}, the goal is to approximate $1-e^{-yt}$ by a
function $h(y)$ whose derivatives at zero all have small magnitude.

\subsection{Estimation via polynomial approximation and support size estimation}
\label{sec:support}
Approximation-theoretic techniques for estimating norms and other
properties such as support size and entropy have been successfully used in the
statistics literature. For example, estimating the $L_p$ norms in Gaussian models
\cite{LNS99,CL11} and  estimating entropy~\cite{WY14,JVHW15} and
support size~\cite{WY14b} of discrete distributions.
Among the aforementioned
problems, support size estimation is closest to ours. Hence, we now
discuss the difference between the approximation technique we use and the
those used for support size estimation.

The support size of a discrete distribution $p$ is
\begin{equation}
S(p) = \sum_x \mathbbm{1}_{p_x > 0}.
	\label{eq:Sp}
\end{equation}
At the first glance, estimating $S(p)$ may
appear similar to species estimation problem as one can convert a
support size estimator $\hat{\supp}$ to $\hat U$ by
\[
\hat{U} = \hat{\supp} - \sum^\infty_{i=1} \prev{i}.
\]
However, without any assumption on the distribution it is impossible to estimate the support size.
For example, regardless how many samples are collected, there could be infinitely many symbols
with arbitrarily small probabilities that will never be observed.  A common assumption is therefore
that the minimum non-zero probability of the underlying distribution $p$, denoted by $p^+_{\min}$,
is at least $1/k$, for some known $k$. Under this assumption~\cite{VV11} used a linear programming estimator similar to the one in \cite{ET76}, to estimate the support size within an additive error of $k \epsilon$ with constant probability using
$\Omega(\frac{k}{\log k} \frac{1}{\epsilon^2})$ samples. Based on {best polynomial approximations} recently \cite{WY14b} showed that the minimax risk
of support size estimation satisfies
\[
\min_{\hat{\supp}} \max_{p: p^+_{\min} \geq 1/k} \mathbb{E}_p
    [(\hat{\supp}-S(p))^2] = k^2 \exp \left(- \Theta\pth{\max \sth{
    \sqrt{\frac{k\log k}{n}}, \frac{k}{n}, 1}} \right)
\]
and that the optimal sample complexity of for estimating $S(p)$ within an additive error of $k \epsilon$ with constant probability is in fact $\Theta(\frac{k}{\log k} \log^2 \frac{1}{\epsilon})$.
Note that the assumption $p^+_{\min} \geq 1/k$ is crucial for this result to hold for otherwise estimation is impossible; in contrast, as we
show later, for species estimation no such assumptions are necessary. The intuition is that if there exist a large number of very improbable symbols, most likely they will not appear in the new samples anyway.

To estimate the support size, in view of \prettyref{eq:Sp} and the assumption $p^+_{\min} \geq 1/k$, the technique of \cite{WY14b} is to approximate the indicator function $y\mapsto\mathbbm{1}_{y \geq 1/k}$ in the
range $\{0\} \cup [1/k, \log k/n]$ using Chebyshev polynomials.
Since by assumption no $p_x$ lies in $(0,\frac{1}{k})$, the approximation error in this interval is irrelevant. For example, in \arxiv{Figure~\ref{fig:fsa}}{(\prettyref{fig:fsa})}, the red curve is a useful
  approximation for the support size, even though it behaves badly over
  $(0,1/k)$. To estimate the average number of unseen symbols $U$, in view of \prettyref{eq:Upoisson}, we need to approximate $y\mapsto 1-e^{-yt}$ over the entire
  $[0,\infty)$ as in, \emph{e.g., } \arxiv{Figure~\ref{fig:fsb}}{(\prettyref{fig:fsb})}.
\begin{figure}[ht]
\centering
\subfigure[]{\label{fig:fsa}
\arxiv{\includegraphics[scale = 0.4]{suppfigure.pdf}}{\includegraphics[scale = 0.4]{suppfigure.pdf}}}
\subfigure[]{\label{fig:fsb}
\arxiv{\includegraphics[scale = 0.4]{specifigure.pdf}}{\includegraphics[scale = 0.4]{specifigure.pdf}}}
\caption{(a) a good approximation for support size; (b) a good
  approximation for species estimation.}
\end{figure}
Concurrent to this work,~\cite{VV15} proposed a linear programming algorithm to
estimate $U$.  However, their NMSE is $O(\frac{t}{\log n})$ compared
to the optimal result $O(n^{-1/t})$ in \prettyref{thm:res_main}, thus exponentially weaker for $t =
o(\log n)$. Furthermore, the computational cost far exceeds those of
our linear estimators
\section{Results for the Poisson model}
\label{sec:upper}
In this section, we provide the performance guarantee for SGT estimators under the Poisson sampling model.  We
first show that the truncated GT estimators incurs a high bias.  We
then introduce the class of smoothed GT estimators
obtained by averaging several truncated GT estimators and bound their
mean squared error in Theorem~\ref{thm:general_bound} for an arbitrary smoothing distribution. We then apply
this result to obtain NMSE bounds for Poisson and Binomial smoothing
in Corollaries~\ref{cor:poi_tail} and~\ref{cor:bin_tail}
respectively, which imply the main result (Theorem~\ref{thm:res_main})
announced in \prettyref{sec:main} for the Poisson model.

\subsection{Why truncated Good-Toulmin does not work}
\label{sec:taylor}
Before we discuss the SGT estimator, we first show that the naive approach of truncating the GT estimator described in \prettyref{sec:SGT}
leads to bad performance when $t>1$.
Recall from \prettyref{lem:general_bounds}
that designing a good linear estimator boils to approximating $1-e^{-yt}$ by an analytic function $h(y)=\sum_{i\geq 1} \frac{h_i}{i!} y^i$ such
that all its derivatives at zero are small, namely, $\sup_{i\geq 1} |h_i|$ is small. The GT estimator
corresponds to the perfect approximation
\[
h^{\scriptscriptstyle\rm GT}(y) = 1-e^{-yt};
\]
however, $\sup_{i\geq 1} |h_i| = \max(t,t^\infty)$, which is infinity if $t > 1$ and leads to large variance.
To avoid this situation, a natural approach is to use use the $\ell$-term Taylor expansion of $1-e^{-yt}$ at $0$, namely,
\begin{equation}
h^\ell(y) = - \sum^{\ell}_{i=1} \frac{(-yt)^i}{i!},
	\label{eq:h-ell}
	\end{equation}
	which corresponds to the estimator $U^\ell$ defined in \prettyref{eq:UGT-ell}.
Then $\sup_{i \geq 1} |h_i| = t^{\ell}$ and, by
Lemma~\ref{lem:general_bounds}, the variance is at most $n (t^{\ell}
+ t)$. Hence if ${\ell} \leq \log_t m$, the variance is at most
$n(m+t)$.  However, note that the $\ell$-term Taylor approximation is a degree-$\ell$ polynomial which eventually diverges and
deviates from $1-e^{-yt}$ as $y$ increases, thereby incurring a large bias. \arxiv{Figure~\ref{fig:taylor}}{(Fig. S2a)} illustrates this phenomenon by plotting the function $1-e^{-yt}$ and its Taylor expansion
with $5,10,$ and $20$ terms.
\begin{figure}[ht]
\centering
\subfigure[]{\label{fig:taylor}
\arxiv{\includegraphics[scale=0.4]{Taylor.pdf}}{\includegraphics[scale=0.4]{Taylor.pdf}}}
\subfigure[.]{\label{fig:average_taylor}
\arxiv{\includegraphics[scale=0.4]{Average_Taylor.pdf}}{\includegraphics[scale=0.4]{Average_Taylor.pdf}}}
\caption{(a) Taylor approximation for $t = 2$, (b) Averages of $10$ and $11$ term Taylor approximation $t = 2$. }
\label{fig:taylor_all}
\end{figure}
Indeed, the next result (proved in \prettyref{app:ell}) rigorously shows that the NMSE of truncated GT estimator never vanishes:
\begin{Lemma}
\label{lmm:UGT-ell}
There exist a constant $c > 0$ such that for any $\ell \geq 0$, any $t>1$ and any $n\in\mathbb{N}$,
	\[
	L^{E^\ell}(n,t) \geq \frac{c(t-1)^5}{ t^4}.
	\]
\end{Lemma}

\subsection{Smoothing by random truncation}
\label{sec:smooth}
As we saw in the previous section, the $\ell$-term Taylor
approximation, where all the coefficients after the $\ell^\mathrm{th}$ term
are set to zero results in large bias. Instead, one can choose a weighted average of several Taylor
series approximations, whose biases cancel each other leading to significant bias reduction.
For example, in \arxiv{Figure~\ref{fig:average_taylor}}{(Fig. S2b)},
 we plot
\[
w  h^{10} + (1-w) h^{11}
\]
for various values of $w\in[0,1]$. Notice that the weight $w=0.6$ leads to better approximation of $1-e^{-yt}$ than both $h^{10}$ and $h^{11}$.

A natural generalization of the above argument entails taking the weighted average of various Taylor approximations
with respect to a given probability distribution over $\ZZ_{+}
\stackrel{\mathrm{def}}{=} \{0,1,2,\ldots\}$. For a $\ZZ_{+}$-valued random variable $L$, consider the power series
\[
h^{\scriptscriptstyle\rm L}(y) = \sum^{\infty}_{\ell=0} \Pr(L = \ell) \cdot h^{\ell}(y),
\]
where $h^\ell$ is defined in	\prettyref{eq:h-ell}.
Rearranging terms, we have
\[
h^{\scriptscriptstyle\rm L}(y) = \sum^\infty_{\ell=0} \Pr(L =\ell) \sum^{\ell}_{i=1} \frac{-(-yt)^i}{i!}
= - \sum^\infty_{i=1}  \frac{(-yt)^i}{i!} \Pr(L \geq i).
\]
Thus, the linear estimator with coefficients
\ignore{
A natural way to reduce the bias
would be to use a $L$-term Taylor series approximation, where
$L$ is distributed according to some known distribution over $\ZZ_{+}
\stackrel{\mathrm{def}}{=} \{0,1,2,\ldots\}$.  Let $U'$ be the estimator whose coefficients
are given by the formal power series
\[
h' = -\sum^{L}_{i=1} \frac{(-yt)^i}{i!}.
\]
However, by Jensen's inequality,
\[
\mathbb{E}[(U' - U)^2] \geq \mathbb{E}[(\mathbb{E}_L[U'] - U)^2].
\]
Thus, instead of using $U'$, we can use $\mathbb{E}_L[U']$ which has a
smaller mean squared loss. Thus, we propose to use estimators of the
form $U^{\scriptscriptstyle\rm L} \stackrel{\mathrm{def}}{=} \mathbb{E}_L[U']$, where $L$ is distributed
according to some known distribution.
\[
h^{\scriptscriptstyle\rm L} = \sum^\infty_{\ell=0} \Pr(L =\ell) \sum^{\ell}_{i=1} \frac{-(-yt)^i}{i!}
= - \sum^\infty_{i=1}  \frac{(-yt)^i}{i!} \Pr(L \geq i).
\]
Thus, randomizing the cutoff point $L$ amounts to attenuating the coefficients of the GT estimator by the tail probability of $L$, resulting in a linear estimator with the following coefficients:
}%
\begin{equation}
	h^{\scriptscriptstyle\rm L}_i = -(-t)^i \prob{L \geq i},
	\label{eq:hrandi}
\end{equation}
is precisely the SGT estimator $U^{\scriptscriptstyle\rm L}$ defined in \prettyref{eq:UL}. Special cases of smoothing distributions include:
\begin{itemize}
\item $L=\infty$: This corresponds to the original Good-Toulmin estimator \prettyref{eq:GT} without smoothing;
	\item $L=\ell$ deterministically: This leads to the estimator $U^\ell$ in \prettyref{eq:UGT-ell} corresponding to the $\ell$-term Taylor approximation;
\item $L\sim \mathrm{Bin}(k, 1/(1+t))$: This recovers the Efron-Thisted estimator \prettyref{eq:UET}, where $k$ is a tuning parameter to be chosen.
\end{itemize}
\begin{figure}[ht]
\centering
\subfigure[]{\label{fig:errors}\includegraphics[scale=0.5]{diff_est_approx-crop.pdf}}
\subfigure[]{\label{fig:coeff}\includegraphics[scale=0.5]{Coeff-crop.pdf}}
\caption{Comparisons of approximations of $h^{\scriptscriptstyle\rm L}(\cdot)$ with  $\mathbb{E}[L] = 2$ and $t=2$. (a) $e^{-y}(1-e^{-yt} - h^{\scriptscriptstyle\rm L}(y))$ as a function of $y$. (b) Coefficients $h^{\scriptscriptstyle\rm L}_i$ as a function of index $i$.}
\end{figure}
We study the performance of linear
estimators corresponding to the Poisson smoothing and the Binomial smoothing. To this
end, we first systematically upper bound the bias and variance for any
probability smoothing $L$.  We plot the error that corresponds to each
smoothing in
\arxiv{Figure~\ref{fig:errors}}{(Fig. S3a)}. Notice that the Poisson and binomial
smoothings have significantly small error compared to the Taylor series
approximation. The coefficients of the resulting estimator is plotted
in
\arxiv{Figure~\ref{fig:coeff}}{(Fig. S3b)}. It is easy so see that the maximum absolute
value of the coefficient is higher for the Taylor series approximation
compared to the Poisson or binomial smoothings.
\begin{Lemma}
For a random variable $L$ over $\ZZ_{+}$ and $t \geq 1$,
\[
\mathrm{Var}(U^{\scriptscriptstyle\rm L}-U) \leq \mathbb{E}[\prev{+}] \cdot  {\mathbb{E}}^2[t^{L}] + \mathbb{E}[U].
\]
\end{Lemma}
\begin{proof}
By Lemma~\ref{lem:general_bounds}, to bound the variance it suffices to bound the highest coefficient in $h^{\scriptscriptstyle\rm L}$.
\begin{equation}
\label{eq:coeff_rand}
|h^{\scriptscriptstyle\rm L}_i| \leq   t^i  \Pr(L \geq i ) = t^i \sum^\infty_{j=i}
\Pr(L =j) \leq   \sum^\infty_{j=i}
\Pr(L =j)  t^j  \leq \mathbb{E}[t^L].
\end{equation}
The above bound together with Lemma~\ref{lem:general_bounds} yields the result.
\end{proof}
To bound the bias, we need few definitions. Let
\begin{equation}
  g(y) \stackrel{\mathrm{def}}{=} - \sum_{i=1}^\infty \frac{\prob{L \geq i}}{i!} (-y)^i.
	\label{eq:EGF}
\end{equation}
Under this definition, $h^{\scriptscriptstyle\rm L}(y) = g(yt)$.
We use the following auxiliary lemma to bound the bias.
\begin{Lemma}
	\label{lem:Delta}
	For any random variable $L$ over $\ZZ_{+}$,
\begin{equation*}
	g(y) - (1-e^{-y}) = -e^{-y} \int_0^y \mathbb{E} \left[{\frac{(-s)^L}{L!}}\right] e^s \text{d} s.
\end{equation*}
\end{Lemma}
\begin{proof}
Subtracting
\prettyref{eq:EGF} from the Taylor series expansion of $1-e^{-y}$
,
	\begin{align*}
g(y) - (1-e^{-y})
= & ~ \sum_{i=1}^\infty \frac{\prob{L < i}}{i!} (-y)^i	\\
= & ~ \sum_{i=1}^\infty \sum_{j=0}^{i-1} \frac{(-y)^i}{i!} \prob{L = j}	\\
= & ~ \sum_{j=0}^\infty \pth{\sum_{i=j+1}^{\infty} \frac{(-y)^i}{i!} } \prob{L = j}.
\end{align*}
Note that $\sum_{i=j+1}^{\infty} \frac{z^i}{i!}$ can be expressed (via incomplete Gamma function) as
\[
\sum_{i=j+1}^{\infty} \frac{z^i}{i!} = \frac{e^{z}}{j!} \int_0^z \tau^j e^{-\tau} \text{d} \tau.
\]
Thus by Fubini's theorem,
	\begin{align*}
g(y) - (1-e^{-y})
= & ~ \sum_{j=0}^\infty \frac{e^{-y}}{j!} \int_0^{-y} \tau^j e^{-\tau} \text{d} \tau \prob{L = j}	\\
= & ~ e^{-y} \int_0^{-y} e^{-\tau}\text{d}  \tau \pth{\sum_{j=0}^\infty \frac{\tau^j}{j!}    \prob{L = j}}	\\
= & ~ -e^{-y} \int_0^{y} e^{s} \text{d}  s \pth{\sum_{j=0}^\infty \frac{(-s)^j}{j!}    \prob{L = j}}	\\
= & ~ -e^{-y} \int_0^y \mathbb{E} \left[{\frac{(-s)^L}{L!}}\right] e^s \text{d} s. \qedhere
\end{align*}
\end{proof}
To bound the bias, we need one more definition. For a random variable $L$ over  $\ZZ_{+}$, let
\[
\xi_{\rand}(t) \stackrel{\mathrm{def}}{=}  \max_{0 \leq s < \infty} \left  \lvert \mathbb{E} \left[{\frac{(-s)^L}{L!}}\right] \right \rvert e^{-s/t},
\]
\begin{Lemma}
\label{lem:rand_bias}
For a random variable $L$ over $\ZZ_{+}$,
\[
|\mathbb{E}[U^{\scriptscriptstyle\rm L} - U]| \leq (\mathbb{E}[\prev{+}] + \mathbb{E}[U])\cdot \xi_{\rand}(t).
\]
\end{Lemma}
\begin{proof}
By Lemma~\ref{lem:Delta},
\begin{align*}
|g(y) - (1-e^{-y})|
& \leq e^{-y} \int_0^y \left  \lvert \mathbb{E} \left[{\frac{(-s)^L}{L!}}\right] \right \rvert e^s \text{d} s \\
& \leq \max_{s \leq y}  \left \lvert \mathbb{E} \left[{\frac{(-s)^L}{L!}}\right]  \right \rvert e^{-y} \int_0^y e^s \text{d} s \\
&=    \max_{s \leq y} \left  \lvert \mathbb{E} \left[{\frac{(-s)^L}{L!}}\right] \right \rvert (1-e^{-y})   .
\end{align*}
For a symbol $x$,
\begin{align*}
e^{-\lambda_x} \left( h^{\scriptscriptstyle\rm L}(\lambda_x) - (1-e^{-\lambda_x t}) \right)
= e^{-\lambda_x} \left( g(\lambda_xt) - (1-e^{-\lambda_x t}) \right).
\end{align*}
Hence,
\begin{align*}
\lvert e^{-\lambda_x} \left( h^{\scriptscriptstyle\rm L}(\lambda_x) - 1-e^{-\lambda_x t} \right) \rvert
& \leq  (1-e^{-\lambda_xt}) \max_{0 \leq y \leq \infty} e^{-y} \max_{0 \leq s \leq yt} \left  \lvert \mathbb{E} \left[{\frac{(-s)^L}{L!}}\right] \right \rvert \\
& \leq  (1-e^{-\lambda_xt}) \max_{0 \leq s \leq \infty} \left  \lvert \mathbb{E} \left[{\frac{(-s)^L}{L!}}\right] \right \rvert e^{-s/t}.
\end{align*}
The lemma follows by summing over all the symbols and substituting
$\sum_x 1-e^{-\lambda_xt} \leq \sum_x 1-e^{-\lambda_x(t+1)} =  \mathbb{E}[\prev{+}] + \mathbb{E}[U]$.
\end{proof}
The above two lemmas yield our main result.
\begin{Theorem}
\label{thm:general_bound}
For any random variable $L$ over $\ZZ_{+}$ and $t \geq 1$,
\[
\mathbb{E}[(U^{\scriptscriptstyle\rm L} - U)^2] \leq \mathbb{E}[\prev{+}] \cdot {\mathbb{E}}^2[t^{L}] + \mathbb{E}[U]
+  (\mathbb{E}[\prev{+}] + \mathbb{E}[U] )^2  \xi_{\rand}(t)^2.
\]
\end{Theorem}
We have therefore reduced the problem of computing mean-squared loss,
to that of computing expectation of certain function of the random
variable. We now apply the above theorem for Binomial and Poisson
smoothings. Notice that the above bound is distribution dependent and can
be used to obtain stronger results for certain distributions. However,
in the rest of the paper, we concentrate on obtaining minimax
guarantees.

\subsection{Poisson smoothing}

\begin{Corollary}
\label{cor:poi_tail}
For $t \geq 1$, $L \sim \mathrm{poi}(r)$ with $r=  \frac{1}{2t} \log  \left(\frac{n(t+1)^2}{t-1} \right)$,
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm L})
\leq \frac{c_t}{n^{1/t}},
\]
where $0 \leq c_t \leq 3$ and $\lim_{t \to\infty} c_t = 1$.
\end{Corollary}
\begin{proof}
For $L \sim \mathrm{poi}(r)$,
\begin{equation}
\label{eq:poivar}
\mathbb{E}[t^L] = e^{-r} \sum^\infty_{\ell=0} \frac{(r t)^{\ell}}{\ell!} = e^{r(t-1)}.
\end{equation}
Furthermore,
\[
\mathbb{E} \left[{\frac{(-s)^L}{L!}}\right] = e^{-r} \sum_{j=0}^\infty \frac{(-s r)^j}{(j!)^2} = e^{-r} J_0(2\sqrt{sr}),
\]
where $J_0$ is the Bessel function of first order which takes values in $[-1,1]$ cf.~\cite[9.1.60]{AS64}. Therefore
\begin{equation}
\label{eq:poibias}
\xi_{\rand}(t) \leq e^{-r}.
\end{equation}
Equations~\eqref{eq:poivar} and~\eqref{eq:poibias} together with Theorem~\ref{thm:general_bound} yields
\[
\mathbb{E}[(U^{\scriptscriptstyle\rm L} - U)^2]  \leq \mathbb{E}[\prev{+}] \cdot e^{2r(t-1)} + \mathbb{E}[U]
+  (\mathbb{E}[\prev{+}]  + \mathbb{E}[U])^2  \cdot  e^{-2r}.
\]
Since $\mathbb{E}[\prev{+}]  \leq n$ and $\mathbb{E}[U] \leq nt$,
\[
\mathbb{E}[(U^{\scriptscriptstyle\rm L} - U)^2]  \leq n e^{2r(t-1)} + nt +  (n+nt)^2 e^{-2r}.
\]
Choosing $r = \frac{1}{2t} \log \frac{n(t+1)^2}{t-1}$ yields
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm L})  \leq \frac{1}{(nt)^{1/t}} \cdot \left(\frac{t(t-1)}{(t+1)^2} \right)^{\frac{1-t}{t}} + \frac{1}{nt},
\]
and the lemma with $c_t \stackrel{\mathrm{def}}{=}  \frac{1}{t^{1/t}} \cdot \left(\frac{t(t-1)}{(t+1)^2} \right)^{\frac{1-t}{t}} + \frac{1}{t}$.
\end{proof}
\subsection{Binomial smoothing}
We now prove the results when $L \sim \mathrm{Bin}(k,q)$. Our analysis holds
for all $q \in [0,2/(2+t)]$ and in this range, the performance of the estimator improves as $q$ increases, and hence the NMSE bounds are strongest for $q=2/(2+t)$.
Therefore, we consider binomial smoothing for two cases: the Efron-Thisted suggested value $q = 1/(1+t)$ and the optimized value $q = 2/(2+t)$.
\begin{Corollary}
\label{cor:bin_tail}
For $t \geq 1$ and $L \sim \mathrm{Bin}(k,q)$, if $k = \left \lceil \frac{1}{2} \log_2  \frac{nt^2}{t-1} \right \rceil \text{ and } q = \frac{1}{t+1}$,
then
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm L}) \leq \frac{c_t}{n^{\log_2(1+1/t)}},
\]
where $c_t$ satisfies $0 \leq c_t \leq 6$ and $\lim_{t \to\infty} c_t = 1$; if $k = \left \lceil \frac{1}{2} \log_3  \frac{nt^2}{t-1} \right \rceil \text{ and } q = \frac{2}{t+2}$,
then
\[
	\cE_{n,t}(U^{\scriptscriptstyle\rm L}) \leq \frac{c'_t}{(nt)^{\log_3(1+2/t)}},
\]
where $c'_t$ satisfies $0 \leq c'_t \leq 6$ and $\lim_{t \to\infty} c'_t = 1$.
\ignore{
\begin{align*}
k = \left \lceil \frac{1}{2} \log_2  \frac{nt^2}{t-1} \right \rceil \text{ and } q = \frac{1}{t+1}
\Rightarrow  & ~ \cE_{n,t}(U^{\scriptscriptstyle\rm L}) \leq \frac{c_t}{n^{\log_2(1+1/t)}},	\\
k = \left \lceil \frac{1}{2} \log_3  \frac{nt^2}{t-1} \right \rceil \text{ and } q = \frac{2}{t+2} \Rightarrow   & ~ 	\cE_{n,t}(U^{\scriptscriptstyle\rm L}) \leq \frac{c'_t}{(nt)^{\log_3(1+2/t)}},
\end{align*}
where $c_t$ satisfies $c_t \leq 6$ for $t \geq 1$ and $\lim_{t \to \infty} c_t = 1$, and
$c'_t $ satisfies $c'_t \leq 11$ for all $t \geq 1$ and $\lim_{t \to \infty} c'_t = 1$.
}
\end{Corollary}
\begin{proof}
If $L\sim \mathrm{Bin}(k,q)$,
\[
\mathbb{E}[t^L] = \sum^k_{\ell=0} {k \choose \ell} (tq)^{\ell}(1-q)^{k - \ell} = (1+q(t-1))^{k}.
\]
Furthermore,
\[
\mathbb{E} \left[{\frac{(-s)^L}{L!}}  \right]
=  \sum_{j=0}^k \frac{(-s)^j}{j!} \binom{k}{j} (q)^j (1-q)^{k-j} = (1-q)^{k}L_k\pth{\frac{qs}{1-q}},
\]
where
\begin{equation}
L_{k}(y) = 	\sum_{j=0}^k \frac{(-y)^j}{j!}  \binom{k}{j}
	\label{eq:lag}
\end{equation}
 is the Laguerre polynomial of degree $k$.
If $\frac{tq}{2(1-q)} \leq 1$, for any $s \geq 0$,
\[
e^{-\frac{s}{t}} \left \lvert \mathbb{E} \left[{\frac{(-s)^L}{L!}}  \right] \right \rvert
\leq (1-q)^{k} e^{-\frac{s}{t}} e^{\frac{qs}{2(1-q)}}
\leq (1-q)^{k},
\]
where the second inequality follows from the fact cf.~\cite[22.14.12]{AS64} that
for all $y \geq 0$ and all $k\geq 0$,
\begin{equation}
|L_k(y)| \leq e^{y/2}.
	\label{eq:lagbound}
\end{equation}
Hence for $q\leq 2/(t+2)$,
\[
\mathbb{E}[(U^{\scriptscriptstyle\rm L} - U)^2]  \leq \mathbb{E}[\prev{+}] \cdot (1+q(t-1))^{2k} + \mathbb{E}[U]
+  (\mathbb{E}[\prev{+}] +\mathbb{E}[U])^2  \cdot  (1-q)^{2k}.
\]
Since $\mathbb{E}[U] \leq nt$ and $\mathbb{E}[\prev{+}] \leq n$,
\begin{equation}
\label{eq:gen_bin}
\mathbb{E}[(U^{\scriptscriptstyle\rm L} - U)^2]  \leq n\cdot (1+q(t-1))^{2k} + nt
+   (nt+n)^2  \cdot  (1-q)^{2k}.
\end{equation}
Substituting the Efron-Thisted suggested $q = \frac{1}{t+1}$ results in
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm L}) \leq \left(\frac{2^{2k}}{nt^2} + \frac{(t+1)^2}{t^2}\right)\left( \frac{t}{t+1}\right)^{2k} + \frac{1}{nt}.
\]
Choosing $k = \left \lceil \frac{1}{2}\log_2 \frac{nt^2}{t-1} \right \rceil $
yields the first result with $c_t \stackrel{\mathrm{def}}{=}  \left(\frac{4}{t-1} + \left(\frac{t+1}{t}\right)^2 \right)  \cdot \left(\frac{t-1}{t^2} \right)^{\log_2(1+1/t)} + \frac{1}{t}$.
For the second result, substituting $q = \frac{2}{t+2}$ in \eqref{eq:gen_bin}
results in
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm L}) \leq \left(\frac{3^{2k}}{nt^2} + \frac{(t+1)^2}{t^2}\right)\left( \frac{t}{t+2}\right)^{2k} + \frac{1}{nt}.
\]
Choosing $k =\left \lceil  \frac{1}{2} \log_3 \frac{nt^2}{t-1} \right \rceil $ yields the result with $c'_t \stackrel{\mathrm{def}}{=} \left(\frac{9}{t-1} + \frac{(t+1)^2}{t^2} \right) \cdot \left(\frac{t-1}{t^2} \right)^{\log_3(1+2/t)} + \frac{1}{t}$ .
\end{proof}
\ignore{
\subsection{Results on the smoothed Good-Toulmin estimators}
Combining the bounds on bias and variance leads to the following upper bound on the normalized mean square loss:
\begin{Theorem}
\label{thm:rupper}
For any strictly positive $n, t = m/n$ and $r$,
the proposed estimator satisfies
\[
L^{\scriptscriptstyle\rm J}_{n}(m) \leq e^{-2r}  + \frac{1}{nt^2}\max\sth{t^{2r},t^2,\frac{e^{2r(t-1)}}{\max\{1, \frac{8e}{27} rt\}}} + \frac{1}{nt}.
\]
\end{Theorem}
Substituting a suitable value of $r$ results in
\begin{Theorem}
\label{thm:main}
Let $t = m/n \geq 1$ and $nt \geq 7$.
Let
\begin{equation}
r = \frac{\log (nt) + \log \log(nt)}{2t}.
	\label{eq:bestr}
\end{equation}
Then
\[
L^{\scriptscriptstyle\rm J}_{n}(m) \leq \frac{1}{nt} + \pth{1 + \max \sth{\frac{5}{2t} , \frac{(n t \log(nt))^{\log(et)/t}}{nt^2}} } \frac{1}{(n t \log(nt))^{1/t}}
\]
The above bound is at most $3/n^{1/t}$ and hence the predictability horizon is at least $M^{\scriptscriptstyle\rm J}_n(\delta) \geq \frac{n \log n}{\log (3/\delta)}$.
\end{Theorem}
\begin{proof}
From Lemmas~\ref{lem:variance} and \ref{lem:bias}, the mean square loss is at most
\begin{align*}
\mathbb{E}_p[(U^{\scriptscriptstyle\rm J}-U)^2]
\leq
 & ~ n t+ e^{-2r} (nt)^2 +  n \max \sth{\frac{e^{2r(t-1)}}{\max\{1, \frac{8e}{27} rt\}}, t^{2r}, t^2} \\
= & ~ n t+ e^{-2r} (nt)^2 + n \max \sth{\frac{27e^{2r(t-1)}}{8e rt}, t^{2r}},
\end{align*}
where the last equality follows from the assumption that $t \geq 1$ and $nt\geq 7$ so that $rt  > 27/(8e)$.
Therefore the normalized mean square loss satisfies
\[
L^{\scriptscriptstyle\rm J}_{n}(m)
\leq   \frac{1}{n t} + e^{-2r} + \frac{1}{nt^2} \max\sth{\frac{27e^{2r(t-1)}}{8e rt} ,  t^{2r}}.
\]
With the choice of $r$ in \prettyref{eq:bestr}, $e^{-2r} = \frac{1}{(n t \log(nt))^{1/t}}$, which is the leading term.
As for the last two terms,
\[
\frac{1}{nt^2} \frac{27e^{2r(t-1)}}{8e rt} \frac{1}{e^{-2r} } = \frac{27}{4et} \frac{\log(nt)}{\log(nt)+\log\log(nt)} \leq \frac{5}{2t}
\]
 and
\[
\frac{1}{nt^2} \frac{t^{2r}}{e^{-2r}}  =  \frac{(n t \log(nt))^{\log(et)/t}}{nt^2},
\]
 completing the proof.
\ignore{
With the choice of $r$ in \prettyref{eq:bestr}, $e^{-2r} = \frac{1}{(n t \log(nt))^{1/t}}$, which is the leading term.
As the last two terms,
$\frac{1}{nt^2} \frac{27e^{2r(t-1)}}{8e rt} \frac{1}{e^{-2r} } = \frac{27}{8e} \frac{\log(nt)}{\log(nt)+\log\log(nt)} \leq \frac{2}{t}$ and
$\frac{1}{nt^2} \frac{t^{2r}}{e^{-2r}}  =  \frac{(n t \log(nt))^{\log(et)/t}}{nt^2}$, completing the proof.
}
\end{proof}
} In terms of the exponent, the  result is strongest for $L \sim
\mathrm{Bin}(k, 2/(t+2))$. Hence, we state the following asymptotic
result, which is a direct consequence of \prettyref{cor:bin_tail}:
\begin{Corollary}
\label{cor:asym}
For  $L \sim \mathrm{Bin}(k,q)$, $q =
\frac{2}{t+2}$,$k = \lceil \log_3 (\frac{nt^2}{t-1})\rceil$, and any fixed $\delta$, the maximum $t$ till which $U^{\scriptscriptstyle\rm L}$ incurs a NMSE of $\delta$ is
\[
\lim_{n \to \infty} \frac{\max\sets{t: \cE_{n,t}(U^{\scriptscriptstyle\rm L})<\delta}}{\log n}
\geq \frac{2}{\log
  3\cdot \log \frac{1}{\delta}}.
\]
\end{Corollary}
\begin{proof}
By Corollary~\ref{cor:bin_tail}, if $t \to\infty$, then
\begin{equation*}
\cE_{n,t}(U^{\scriptscriptstyle\rm L}) \leq (1+o(1)) n^{- \frac{2+o(1)}{t\log 3}}.
	\label{eq:convrate}
\end{equation*}
where $o(1)=o_t(1)$ is uniform in $n$.
Consequently, if $t = (\alpha + o(1)) \log n$ and $n\to\infty$, then
\[
\limsup_{n \to \infty} \cE_{n,t}(U^{\scriptscriptstyle\rm L}) \leq e^{-\frac{2}{\alpha\log3}}.
\]
Thus for any fixed $\delta$, the maximum $t$ till which $U^{\scriptscriptstyle\rm L}$ incurs a NMSE of $\delta$ is
\[
\lim_{n \to \infty} \frac{\max\sets{t: \cE_{n,t}(U^{\scriptscriptstyle\rm L})<\delta}}{\log n}
\geq \frac{2}{\log
  3\cdot \log \frac{1}{\delta}}. \qedhere
\]
\end{proof}
Corollaries~\ref{cor:poi_tail} and~\ref{cor:bin_tail} imply Theorem~\ref{thm:res_main}
for the Poisson model.
\section{Extensions to other models}
\label{sec:extensions}
Our results so far have been developed for the Poisson
model. Next we extend them to the multinomial model (fixed sample size), the Bernoulli-product model,
and the hypergeometric model (sampling without replacement)
\cite{BF93}, for which upper bounds of NMSE for general smoothing distributions that are
analogous to \prettyref{thm:general_bound} are presented in Theorem~\ref{thm:multinomial}, \ref{thm:bernoulli} and \ref{thm:hyper}, respectively.
Using these results, we obtain the NMSE for Poisson and Binomial
smoothings similar to Corollaries~\ref{cor:poi_tail}
and~\ref{cor:bin_tail}.  We remark that up to multiplicative
constants, the NMSE under multinomial and Bernoulli-product model are
similar to those of Poisson model; however, the NMSE under
hypergeometric model is slightly larger.

\subsection{The multinomial model}
The multinomial model corresponds to the setting described in \prettyref{sec:intro}, where upon observing $n$ i.i.d. samples, the objective is
to estimate the expected number of new symbols $\UXnm$ that would be observed
if we took $m$ more samples.
We can write the expected
number of new symbols as
\begin{equation*}
\label{eq:binomial}
\UXnm = \sum_{x} \mathbbm{1}_{\mul{x} = 0} \cdot \mathbbm{1}_{\mul{x}' > 0}.
\end{equation*}
As before we abbreviate \[
U \stackrel{\mathrm{def}}{=} \UXnm
\]
and similarly $U^{\scriptscriptstyle\rm E}\stackrel{\mathrm{def}}{=}\UEXnm$ for any estimator $E$.
The difficulty in handling multinomial distributions is that, unlike the Poisson model, the number of
occurrences of symbols are correlated; in particular, they sum up to
$n$. This dependence renders the analysis cumbersome.  In the multinomial
setting each symbol is distributed according to $\mathrm{Bin}(n,p_x)$ and hence
\[
\mathbb{E}[\mathbbm{1}_{\mul{x} = i}] = {n \choose i} p^i_x (1-p_x)^{n-i}.
\]
As an immediate consequence,
\[
\mathbb{E}[\prev{i}] = \mathbb{E}\left[\sum_{x} \mathbbm{1}_{\mul{x} = i} \right] =
\sum_{x} {n \choose i} p^i_x (1-p_x)^{n-i}.
\]
We now bound the bias and variance of an arbitrary linear estimator $U^{\scriptscriptstyle\rm h}$.
 We
first show that the bias
$\mathbb{E}[U^{\scriptscriptstyle\rm h} - U]$ under the multinomial model is close to that
under the Poisson model, which is $\sum_{x} e^{-\lambda_x} (
h(\lambda_x) - (1-e^{-t\lambda_x}) )$ as given in \prettyref{eq:Upoisson}.
\begin{Lemma}
\label{lem:bias_multi}
The bias of $U^{\scriptscriptstyle\rm h} = \sum^\infty_{i=1} \prev{i} h_i$ satisfies
\[
\left\lvert \mathbb{E}[U^{\scriptscriptstyle\rm h} - U] - \sum_{x} e^{-\lambda_x} \left( h(\lambda_x) - (1-e^{-t\lambda_x}) \right) \right \rvert \leq 2\sup_i |h_i| + 2.
\]
\end{Lemma}
\begin{proof}
First we recall a result on Poisson approximation: For $X\sim \mathrm{Bin}(n,p)$ and $Y\sim \mathrm{poi}(np)$,
\begin{equation}
|\mathbb{E}[f(X)] - \mathbb{E}[f(Y)]| \leq 2p \sup_{i} |f(i)|,
	\label{eq:poi-approx}
\end{equation}
which follows from the total variation bound $d_{\rm TV}(\mathrm{Bin}(n,p),\mathrm{poi}(np)) \leq p$ \cite[Theorem 1]{BH84} and the fact that $d_{\rm TV}(\mu,\nu) = \frac{1}{2}\sup_{\|f\|_\infty\leq1} \int f d\mu-\int f d\nu$. In particular, taking $f(x) = \mathbbm{1}_{x=0}$ gives
\begin{equation*}
0 \leq  e^{-np} - (1-p)^n \leq 2 p.
\end{equation*}
Note that the linear estimator can be expressed as $U^{\scriptscriptstyle\rm h} = \sum_x h_{\mul{x}}$.
Under the multinomial model,
\begin{equation*}
\mathbb{E}[U^{\scriptscriptstyle\rm h} - U] = \sum_x \mathbb{E}_{\mul{x}\sim \mathrm{Bin}(n,p_x)}[h_{\mul{x}}] - \sum_x (1-p_x)^n (1-(1-p_x)^m).
	\label{eq:bias-hsub}
\end{equation*}
Under the Poisson model,
\[
\sum_{x} e^{-\lambda_x} \left( h(\lambda_x) - (1-e^{-t\lambda_x}) \right) = \sum_x \mathbb{E}_{\mul{x}\sim \mathrm{poi}(np_x)}[h_{\mul{x}}] - \sum_x e^{-np_x} (1-e^{-mp_x}).
\]
Then
\[
\left|\sum_x \mathbb{E}_{\mul{x}\sim \mathrm{Bin}(n,p_x)}[h_{\mul{x}}] - \sum_x \mathbb{E}_{\mul{x}\sim \mathrm{poi}(np_x)}[h_{\mul{x}}] \right| \overset{\prettyref{eq:poi-approx}}{\leq} 2 \sup_i |h_i| \sum_x p_x = 2 \sup_i |h_i|.
\]
Furthermore,
\begin{align*}
	& ~ \sum_x (1-p_x)^n (1-(1-p_x)^m) - \sum_x e^{-np_x} (1-e^{-mp_x}) \\
\leq & ~ \sum_x e^{-n p_x} ( e^{-mp_x} - (1-p_x)^m) \overset{\prettyref{eq:poi-approx0}}{\leq} \sum_x e^{-n p_x} 2 p_x \leq 2.
\end{align*}
Similarly, $ \sum_x (1-p_x)^n (1-(1-p_x)^m) - \sum_x e^{-np_x} (1-e^{-mp_x}) \geq -2$. Assembling the above proves the lemma.
\end{proof}

The next result bounds the variance.
\begin{Lemma}
\label{lem:variance_multi}
For any linear estimator $U^{\scriptscriptstyle\rm h}$,
\[
\mathrm{Var}(U^{\scriptscriptstyle\rm h} - U) \leq 8 n \max\sth{\sup_{i\geq 1} h^2_i,1} + 8m.
\]
\end{Lemma}
\begin{proof}
Recognizing that
$U^{\scriptscriptstyle\rm h} - U$ is a function of $n+m$ independent random variables, namely,  $X_1,\ldots,X_{n+m}$ drawn i.i.d. from $p$, we apply Steele's variance inequality~\cite{S86} to bound its variance.
Similar to \prettyref{eq:bias-hsub},
\[
U^{\scriptscriptstyle\rm h} - U = \sum_{x}  h_{\mul{x}} + \mathbbm{1}_{\mul{x}=0} \mathbbm{1}_{\mul{x}'> 0}
\]
Changing the value of any one of the first $n$ samples changes the multiplicities of
two symbols, and hence the value of $U^{\scriptscriptstyle\rm h} - U$ can change by at most $4\max(\max_{i \geq 1} |h_i|,1)$.
Similarly, changing any one of the last $m$ samples changes the value of $U^{\scriptscriptstyle\rm h}-U$ by at most four.
Applying Steele's inequality gives the lemma.
\end{proof}
Lemmas~\ref{lem:bias_multi} and~\ref{lem:variance_multi} are analogous to Lemma~\ref{lem:general_bounds}. Together with \eqref{eq:coeff_rand} and Lemma~\ref{lem:rand_bias}, we obtain the main result for the multinomial model.
\begin{Theorem}
\label{thm:multinomial}
For $t \geq 1$ and any random variable $L$ over $\ZZ_{+}$,
\[
\mathbb{E}[(U^{\scriptscriptstyle\rm L}-U)^2] \leq 8 n\, {\mathbb{E}}^2[t^L] + 8m + \left((n(t+1) \xi_{\rand}(t)
+ 2 \mathbb{E}[t^L] + 2\right)^2.
\]
\end{Theorem}
Similar to Corollaries~\ref{cor:poi_tail} and~\ref{cor:bin_tail}, one can compute the NMSE for Binomial and Poisson smoothings. We remark that up to multiplicative constants the results are identical to those for the Poisson model.

\subsection{Bernoulli-product model}
\label{sec:bernoulli}
Consider the following species assemblage model. There are $k$
distinct species and each one can be found in one of $n$ independent
sampling units. Thus every species can be present in multiple sampling
units simultaneously and each sampling unit can capture multiple
species.  For example species $x$ can be found in sampling units $1,3$
and $5$ and species $y$ can be found in units $2,3$, and $4$.  Given
the data collected from $n$ sampling units, the objective
is to estimate the expected number of new species that would be
observed if we placed $m$ more units.

The aforementioned problem is typically modeled as by the
\emph{Bernoulli-product model}. Since, in this model each sample only
has presence-absence data, it is often referred to as incidence
model~\cite{CCG12}.  For notational simplicity, we use the same
notation as the other three models.  In Bernoulli-product model, for a
symbol $x$, $\mul{x}$ denotes the number of sampling units in which
$x$ appears and $\prev{i}$ denotes the number of symbols that appeared
in $i$ sampling units.
Given a set of distinct symbols (potentially infinite), each symbol
$x$ is observed in each sampling unit independently with probability
$p_x$ and the observations from each sampling unit are independent of
each other.
To distinguish from the multinomial and Poisson sampling models where
each sample can be only one symbol, we refer to samples here as
sampling units.  Given the results of $n$ sampling units, the goal is
to estimate the expected number of new symbols that would appear in
the next $m$ sampling units. Let $p_{_S} = \sum_{x} p_x$.  Note that
$p_{_S}$ is also the expected number of symbols that we observe for
each sampling unit and need not sum to $1$. For example, in the
species application, probability of catching bumble bee can be $0.5$
and honey bee be $0.7$.

This model is significantly different from the multinomial model in
two ways. Firstly, here given $n$ sampling units the number of
occurrences of symbols are independent of each other. Secondly, $p_{_S}
\stackrel{\mathrm{def}}{=} \sum_x p_x$ need not be $1$.  In the Bernoulli-product model, the
probability observing each symbol at a particular sample is $p_x$ and
hence in $n$ samples, the number of occurrences is distributed
$\mathrm{Bin}(n,p_x)$.
Therefore the probability that $x$ is be observed in $i$ sampling units is
\[
\mathbb{E}[\mathbbm{1}_{\mul{x} = i}] = {n \choose i} p^i_x (1-p_x)^{n-i},
\]
and an immediate consequence on the number of distinct symbols that
appear $i$ sampling units is
\[
\mathbb{E}[\prev{i}] = \mathbb{E}\left[\sum_{x} \mathbbm{1}_{\mul{x} = i} \right] =
\sum_{x} {n \choose i} p^i_x (1-p_x)^{n-i}.
\]
Furthermore, the expected total number of symbols is $np_{_S}$ and hence
\begin{equation*}
\sum^n_{i=1} \mathbb{E}[\prev{i}] i = n p_{_S}.
\end{equation*}
Under the Bernoulli-product model the objective is to estimate the
number of new symbols that we observe in $m$ more sampling
units and is
\[
\UXnm = \sum_{x} \mathbbm{1}_{\mul{x} = 0} \cdot \mathbbm{1}_{\mul{x}' > 0}.
\]
As before, we abbreviate \[
U \stackrel{\mathrm{def}}{=} \UXnm
\]
and similarly $U^{\scriptscriptstyle\rm E}\stackrel{\mathrm{def}}{=}\UEXnm$ for any estimator $E$.
Since the probabilities need not add up to $1$, we redefine our definition of $\cE_{n,t}(U^{\scriptscriptstyle\rm E})$ as
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm E})
\stackrel{\mathrm{def}}{=}
\max \mathbb{E}_p\Paren{\frac{U-U^{\scriptscriptstyle\rm E}}{n t p_{_S}}}^2.
\]
Under this model, the SGT estimator satisfy similar results to that of Corollaries~\ref{cor:poi_tail} and~\ref{cor:bin_tail}, up to multiplicative constants.
The main ingredient is to bound the bias and variance (like Lemma~\ref{lem:general_bounds}).
We note that since the marginal of $N_x$ is $\mathrm{Bin}(n,p_x)$ under both the multinomial and the Bernoulli-product model, the bias bound follows entirely analogously as in \prettyref{lem:bias_multi}. The proof of variance bound is very similar to that of \prettyref{lem:general_bounds} and hence is omitted.
\begin{Lemma}
\label{lem:bias_bernoulli}
The bias of the linear estimator $U^{\scriptscriptstyle\rm h}$ is
\[
\left\lvert \mathbb{E}[U^{\scriptscriptstyle\rm h} - U] - \sum_{x} e^{-\lambda_x} \left( h(\lambda_x) - (1-e^{-t\lambda_x}) \right) \right \rvert \leq 2 p_{_S} \pth{\sup_i |h_i| + 1},
\]
and the variance
\[
\mathrm{Var}(U^{\scriptscriptstyle\rm h} - U) \leq n p_{_S}  \cdot \left(t + \sup_{i \geq 1} {h^2_i} \right).
\]
\end{Lemma}
The above lemma together with \eqref{eq:coeff_rand} and Lemma~\ref{lem:rand_bias} yields  the main result for the Bernoulli-product model.
\begin{Theorem}
\label{thm:bernoulli}
For any random variable $L$ over $\ZZ_{+}$ and $t \geq 1$,
\[
\mathbb{E}[(U^{\scriptscriptstyle\rm L}-U)^2] \leq  n p_{_S}  \cdot \left(t +{\mathbb{E}}^{2}[t^{L}] \right) + (n(t+1)p_{_S} \xi_{\rand}(t) + 2 p_{_S}(\mathbb{E}[t^{L}] + 1))^2.
\]
\end{Theorem}
Similar to Corollaries~\ref{cor:poi_tail} and~\ref{cor:bin_tail}, one can compute the normalized mean squared loss for Binomial and Poisson smoothings. We remark that up to multiplicative constants the results would be similar to that for the Poisson model.

\subsection{The hypergeometric model}
\label{sec:hyper}
The hypergeometric model considers the population estimation problem with samples drawn without replacement.
Given $n$ samples drawn uniformly at random, without replacement from a set $\{y_1,\ldots, y_{R}\}$ of $R$ symbols,
the objective is to estimate the number of new symbols that would be
observed if we had access to $m$ more random samples without replacement, where $n+m \leq R$.
Unlike the Poisson,
multinomial, and Bernoulli-product models we have considered so far, where the samples are
independently and identically distributed, in the hypergeometric model the samples are \emph{dependent} hence a modified analysis is needed.

Let $r_x \stackrel{\mathrm{def}}{=} \sum^R_{i=1} \mathbbm{1}_{y_i = x}$ be the number of occurrences of
symbol $x$ in the $R$ symbols, which satisfies $\sum _x r_x=R$. Denote by $N_x$ the number of times $x$ appears in the $n$ samples drawn without replacements, which is
distributed according to the hypergeometric distribution $\mathrm{Hyp}(R,r_x,n)$ with the following probability mass function:\footnote{We adopt the convention that $\binom{n}{k}=0$ for all $k<0$ and $k>n$ throughout.}
\[
\mathbb{P}(\mul{x} = i) = \frac{{r_x \choose i}{R-r_x \choose n-i}}{{R \choose n}}.
\]
We also denote the joint distribution of $\{N_x\}$, which is multivariate hypergeometric, by $\mathrm{Hyp}(\{r_x\},n)$.
Consequently,
\[
\mathbb{E}[\prev{i}] = \sum_{x} \mathbb{P}(\mul{x} = i) = \sum_x \frac{{r_x \choose i}{R-r_x \choose n-i}}{{R \choose n}}.
\]
Furthermore, conditioned on $\mul{x} = 0$, $\mul{x}'$ is distributed as $\mathrm{Hyp}(R-n,r_x,m)$
and hence
\begin{equation}
\label{eq:Umean-without_replacement}
\mathbb{E}[U] = \sum_{x} \mathbb{E}[\mathbbm{1}_{\mul{x}=0}] \cdot  \mathbb{E}[\mathbbm{1}_{\mul{x}' > 0} | \mathbbm{1}_{\mul{x} = 0}]
= \sum_{x} \frac{{R-r_x \choose n}}{{R \choose n}} \cdot \left( 1- \frac{{R-n-r_x \choose m}}{{R-n \choose m}} \right).
\end{equation}
As before, we abbreviate \[
U \stackrel{\mathrm{def}}{=} \UXnm
\]
which we want to estimate and similarly for any estimator $U^{\scriptscriptstyle\rm E}\stackrel{\mathrm{def}}{=}\UEXnm$.
We now bound the variance and bias of a linear estimator $U^{\scriptscriptstyle\rm h}$ under the hypergeometric model.
\begin{Lemma}
\label{lem:hyper_var}
For any linear estimator $U^{\scriptscriptstyle\rm h}$,
\[
\mathrm{Var}(U^{\scriptscriptstyle\rm h} - U) \leq 12 n \sup_{i} h^2_i + 6n + 3m.
\]
\end{Lemma}
\begin{proof}
We first note that for a random variable $Y$ that lies in the
interval $[a,b]$,
\[
\mathrm{Var}(Y) \leq \frac{(a-b)^2}{4}.
\]

For notational convenience define $h_0=0$.
Then $U^{\scriptscriptstyle\rm h} =\sum_x h_{N_x}$.
Let $Z = \sum\mathbbm{1}_{N_x=0}$ and $Z' = \sum\mathbbm{1}_{N_x=N_x'=0}$ denote the number of unobserved symbols in the first $n$ samples and the total $n+m$ samples, respectively. Then $U=Z-'Z$. Since the collection of random variables
$\mathbbm{1}_{\mul{x}=0}$ indexed by $x$ are negatively correlated, we have
\[
\mathrm{Var}\bigl(Z) \leq \sum_{x} \mathrm{Var}(\mathbbm{1}_{\mul{x}=0}\bigr) = \sum_{x} \mathbb{E}[ \mathbbm{1}_{\mul{x}=0}(1 - \mathbbm{1}_{\mul{x}=0})] \leq \sum_x \mathbb{E}\qth{\mathbbm{1}_{\mul{x}>0}} \leq n.
\]
Analogously, $\mathrm{Var}(Z') \leq n+m$ and hence
\[
\mathrm{Var}(U^{\scriptscriptstyle\rm h}-U) = \mathrm{Var}(U^{\scriptscriptstyle\rm h}-Z+Z') \leq 3 \mathrm{Var}(U^{\scriptscriptstyle\rm h}) + 3 \mathrm{Var}(Z') + 3  \mathrm{Var}(Z) \leq 3 \mathrm{Var}(U^{\scriptscriptstyle\rm h}) + 6 n + 3m.
\]
Thus it remains to show
\begin{equation}
\mathrm{Var}(U^{\scriptscriptstyle\rm h}) \leq 4 n \sup_{i} h^2_i.
	\label{eq:varUh}
\end{equation}
 By induction on $n$, we show that for any $n \in \mathbb{N}$,
any set of nonnegative integers $\{r_x\}$ and any function $(x,k)\mapsto f(x,k)$ with $k\in\mathbb{Z}_+$ satisfying $f(x,0)= 0$,
\begin{equation}\label{eq:var-induction}
\mathrm{Var}\pth{\sum_x f(x,N_x) } \leq 4 n \|f\|_\infty^2,
\end{equation}
where $\{N_x\} \sim \mathrm{Hyp}(\{r_x\},n)$ and $\|f\|_\infty = \sup_{x,k} |f(x,k)|$. Then the desired Equation~\prettyref{eq:varUh} follows from \prettyref{eq:var-induction} with $f(x,k) = h_k$.

 We first prove \prettyref{eq:var-induction} for $n=1$, in which case
 exactly one of $\mul{x}$'s is one and the rest are zero.
 Hence,
$|\sum_x f(x,N_x)| \leq \|f\|_\infty$ and $\mathrm{Var}(\sum_x f(x,N_x)) \leq \|f\|_\infty^2$.

Next assume the induction hypothesis holds for $n-1$.
Let $X_1$ denote the first sample and let $\tilde N_x$ denote the number of occurrences of symbol $x$ in samples $X_2,\ldots,X_n$.
Then $N_x = \tilde N_x +\mathbbm{1}_{X_1=x}$. Furthermore,  conditioned on $X_1= y$,  $\{\tilde N_{x}\} \sim \mathrm{Hyp}(\{\tilde r_x\}, n-1)$, where $\tilde r_x=r_x-\mathbbm{1}_{x=y}$.
By the
law of total variance, we have
\begin{equation}
\mathrm{Var}\pth{\sum_x f(x,N_x)}=
\mathbb{E} \qth{V(X_1)} + \mathrm{Var} \pth{g(X_1)}.
	\label{eq:totalvar}
\end{equation}
where
\[
V(y) \stackrel{\mathrm{def}}{=} \mathrm{Var}\pth{\sum_x f(x,N_x) \Bigg\vert X_1=y}, \quad g(y) \stackrel{\mathrm{def}}{=} \mathbb{E}\qth{\sum_x f(x,N_x) \Bigg\vert
  X_1=y}
\]
For the first term in \prettyref{eq:totalvar}, note that
\[
V(y) = \mathrm{Var}\pth{\sum_x f(x,\tilde N_x+\mathbbm{1}_{x=y}) \Bigg\vert X_1=y}  =  \mathrm{Var}\pth{\sum_x f_y(x,\tilde N_x) \Bigg\vert X_1=y}.
\]
where we defined $f_y(x,k) \stackrel{\mathrm{def}}{=} f(x,k+\mathbbm{1}_{x=y})$. Hence, by the induction hypothesis,
$V(y) \leq 4(n-1)\|f_y\|_\infty^2 \leq 4(n-1)\|f\|_\infty^2$ and
$\mathbb{E}  \qth{ V(X_1)} \leq  4(n-1) ||f||^2_\infty$.

For the second term in \prettyref{eq:totalvar}, observe that for any $y\neq z$
\[
g(y)= \mathbb{E}[f(y,\tilde N_y+1)|X_1=y] + \mathbb{E}[f(z,\tilde N_{z})|X_1=y] + \mathbb{E}\qth{\sum_{x \neq y,z} f(x, \tilde N_{x})\Bigg|
  X_1=y},
\]
and
\[
g(z)= \mathbb{E}[f(z,\tilde N_z+1)|X_1=z] + \mathbb{E}[f(y,\tilde N_{y})|X_1=z] + \mathbb{E}\qth{\sum_{x \neq y,z} f(x, \tilde N_{x})\Bigg|
  X_1=z},
\]
Observe that $\{N_x\}_{x \neq y,z}$ have the same joint distribution conditioned on either $X_1=y$ or $X_1=z$ and hence
$\mathbb{E}[\sum_{x \neq y,z} f(x, \tilde N_{x})|
  X_1=y] = \mathbb{E}[\sum_{x \neq y,z} f(x, \tilde N_{x})|
  X_1=z]$.
  Therefore $|g(y) - g(z)| \leq 4 \|f\|_\infty$ for any $y\neq z$.
  This implies that the function $g$ takes values in an interval of length at most $4 \|f\|_\infty$.
  Therefore $\mathrm{Var}(g(X_1)) \leq \frac{1}{4} (4 \|f\|_\infty)^2 = 4 \|f\|_\infty^2$.
  This completes the proof of \prettyref{eq:var-induction} and hence the lemma.
  \end{proof}

Let
\begin{align*}
B(h,r_x) \stackrel{\mathrm{def}}{=} \sum^{r_x}_{i=1} {r_x \choose i} \left(\frac{n}{R} \right)^i
\left(1-\frac{n}{R} \right)^{r_x-i} h_i - \left(1-\frac{n}{R}
\right)^{r_x} \left(1 - \left(1-\frac{m}{R-n} \right)^{r_x}\right).
\end{align*}
To bound the bias, we first prove an auxiliary result.
\begin{Lemma}
\label{lem:bin-hyper}
For any linear estimator $U^{\scriptscriptstyle\rm h}$,
\begin{align*}
 \left \lvert
\mathbb{E}[U^{\scriptscriptstyle\rm h} - U] -  \sum_{x}  B(h,r_x)
 \right \rvert  \leq 4\max\pth{\sup_{i} |h_i|,1} + \frac{2R}{R-n}.
\end{align*}
\end{Lemma}
\begin{proof}
Recall that $N_x\sim\mathrm{Hyp}(R,r_x,n)$. Let $\tilde{N}_x$ be a random variable distributed as $\mathrm{Bin}(r_x,n/R)$.
Since $\mathrm{Hyp}(R,r_x,n)$ coincides with $\mathrm{Hyp}(R,n,r_x)$, we have
\[
d_{\rm TV}(\mathrm{Bin}(r_x,n/R),\mathrm{Hyp}(R,r_x,n))
= d_{\rm TV}(\mathrm{Bin}(r_x,n/R),\mathrm{Hyp}(R,n,r_x))
\leq \frac{2r_x}{R},
\]
where the last inequality follows from \cite[Theorem 4]{DF80}.
 Since $d_{\rm TV}(\mu,\nu) = \frac{1}{2}\sup_{\|f\|_\infty\leq1} \int f d\mu-\int f d\nu = \sup_E \mu(E)-\nu(E)$, we have
\begin{equation}
\left|\mathbb{E}[f(N_x)] - \mathbb{E}[f(\tilde{N}_x)]\right| \leq \frac{4r_x}{R} \sup_{i} |f(i)|,
	\label{eq:hyper-approx}
\end{equation}
and
\begin{equation}
\left \lvert \frac{{R-n-r_x \choose m}}{{R-n \choose m}} - \left(1-\frac{m}{R-n}\right)^{r_x} \right \rvert
\leq d_{\rm TV}(\mathrm{Bin}(r_x,m/(R-n)),\mathrm{Hyp}(R-n,m,r_x)) \leq \frac{2r_x}{R-n}.
	\label{eq:poi-approx0}
\end{equation}
Define $f_x(i) = h_i - \mathbbm{1}_{i=0} \left(1 - \left(1-\frac{m}{R-n}\right)^{r_x} \right)$. In view of \prettyref{eq:Umean-without_replacement} and the fact that $\sum r_x=R$, we have
\[
\left \lvert \mathbb{E}[U^{\scriptscriptstyle\rm h} - U] - \sum_{x} \mathbb{E}[f_x(N_x)]\right \rvert \leq \frac{2R}{R-n}.
\]
Applying \prettyref{eq:hyper-approx} yields
\[
\sum_{x} \left \lvert  \mathbb{E}[f_x(\tilde{N}_x) ]
- \mathbb{E}\left[f_x(\mul{x})\right]\right \rvert \leq 4 \sup_{i} |f_x(i)| \leq
4\max\pth{\sup_{i} |h_i|,1}.
\]
The above equation together with \eqref{eq:poi-approx0} results in the lemma since $B(h,r_x)=\mathbb{E}[f_x(\tilde N_x)]$.
\end{proof}
Note that to upper bound the bias, we need to bound $\sum_x B(h,r_x)$.
It is easy to verify for the GT coefficients $h^{\scriptscriptstyle\rm GT}_i = -\left(-t \right)^{i}$ with $t=m/n$, $B(h^{\scriptscriptstyle\rm GT},r_x) = 0$. Therefore, if we choose $h=h^{\scriptscriptstyle\rm L}$ based on the tail of random variable $L$ with $h^{\scriptscriptstyle\rm L}_i = h^{\scriptscriptstyle\rm GT}_i \prob{L \geq i}$ as defined in \prettyref{eq:hrandi}, we have
\begin{align}
B(h^{\scriptscriptstyle\rm L},r_x)
= & ~ \sum^{r_x}_{i=1} {r_x \choose i} \left(\frac{n}{R} \right)^i
\left(1-\frac{n}{R} \right)^{r_x-i} (-t)^i \mathbb{P}(L < i)	\nonumber \\
= & ~ \left(1-\frac{n}{R} \right)^{r_x}
\sum^{r_x}_{i=1} {r_x \choose i} \left(- \frac{m}{R-n} \right)^i \mathbb{P}(L < i).	\label{eq:Btail}
\end{align}

Similar to \prettyref{lem:Delta}, our strategy is to find an integral presentation of the bias. This is done in the following lemma.
\begin{Lemma}
\label{lem:hyper-equal}
For any $y \geq 0$ and any $k \in \mathbb{N}$,
\begin{equation}
\sum^{k}_{i=1} {k \choose i} (-y)^i \mathbb{P}(L < i) =  - k(1-y)^{k}   \int^{y}_0 \mathbb{E} \left[{k-1 \choose L} (-s)^L \right] (1-s)^{-k-1} d s.
	\label{eq:hyper-equal}
\end{equation}
\end{Lemma}
\begin{remark}
\label{rmk:hyper-equal}
	For the special case of $y=1$, \prettyref{eq:hyper-equal} is understood in the limiting sense: Letting $\delta=1-y$ and $\beta=\frac{1-s}{\delta}$, we can rewrite the right-hand side as
	\[- k \int^{1/\delta}_1 \mathbb{E} \left[{k-1 \choose L} (\beta\delta-1)^L \right] k \beta^{-k-1} d \beta.
\]
	For all $|\delta| \leq 1$ and hence $0 \leq 1 - \beta \delta \leq 2$,
	 we have \[
\left|\mathbb{E} \left[{k-1 \choose L} (\beta\delta-1)^L \right]\right| = \Big|\mathbb{E} \left[{k-1 \choose L} (\beta\delta-1)^L \mathbbm{1}_{L < k} \right]\Big| \leq 4^k.
\] 	By dominated convergence theorem, as $\delta\to0$, the right-hand side converges to $- \mathbb{E}\qth{\binom{k-1}{L} (-1)^L}$ and coincides with the left-hand side, which can be easily obtained by applying $\binom{k}{i}=\binom{k-1}{i} + \binom{k-1}{i-1}$.
\end{remark}

\begin{proof}
Denote the left-hand side of \prettyref{eq:hyper-equal} by $F(y)$. Using $i\binom{k}{i}=k\binom{k-1}{i-1}$, we have
\begin{align}
F'(y)
= & ~ \sum^{k}_{i=1} {k \choose i} (-i)(-y)^{i-1} \mathbb{P}(L < i) 	= - k \sum^{k}_{i=1} {k-1 \choose i-1} (-y)^{i-1} \mathbb{P}(L < i) 	\nonumber \\
= & ~ - k \sum^{k}_{i=1} {k-1 \choose i-1} (-y)^{i-1} \mathbb{P}(L < i-1) 	- k \sum^{k}_{i=1} {k-1 \choose i-1} (-y)^{i-1} \mathbb{P}(L = i-1). \label{eq:Fprime}
\end{align}
The second term is simply $- k \mathbb{E} \left[{k-1 \choose L} (-y)^L \right] \stackrel{\mathrm{def}}{=} G(y)$. For the first term, since $L \geq 0$ almost surely and  $\binom{k}{i}=\binom{k-1}{i} + \binom{k-1}{i-1}$, we have
\begin{align}
k \sum^{k}_{i=1} {k-1 \choose i-1} (-y)^{i-1} \mathbb{P}(L < i-1)
= & ~ k \sum^{k}_{i=1} {k-1 \choose i} (-y)^{i} \mathbb{P}(L < i)
	\nonumber \\
= & ~ k \sum^{k}_{i=1} {k \choose i} (-y)^{i} \mathbb{P}(L < i) - k \sum^{k}_{i=1} {k-1 \choose i-1} (-y)^{i} \mathbb{P}(L < i) \nonumber	\\
= & ~ k F(y) - y F'(y). \label{eq:Fprime2}
\end{align}
Combining \prettyref{eq:Fprime} and \prettyref{eq:Fprime2} yields the following ordinary differential equation:
\[
F'(y)(1-y) + kF(y) = G(y), \quad F(0)=0,
\]
whose solution is readily obtained as $F(y) = (1-y)^{k} \int_0^y (1-s)^{-k-1} G(s) d s$, \emph{i.e., } the desired Equation~\prettyref{eq:hyper-equal}.
\end{proof}
Combining \prettyref{lem:bin-hyper}--\ref{lem:hyper-equal}
yields the following bias bound:
\begin{Lemma}
\label{lem:hyper_bias}
For any random variable $L$ over $\ZZ_{+}$ and $t=m/n \geq 1$,
\[
|\mathbb{E}[U^{\scriptscriptstyle\rm L} - U] | \leq nt \cdot \max_{0 \leq s \leq 1} \left \lvert  \mathbb{E} \left[{r_x-1 \choose L} (-s)^L \right] \right \rvert
+  4 \mathbb{E}[t^L] + \frac{2R}{R-n}.
\]
\end{Lemma}
\begin{proof}
Recall the coefficient bound \prettyref{eq:coeff_rand} that $\sup_i |h_i| \leq \mathbb{E}[t^L]$.
By Lemma~\ref{lem:bin-hyper} and the assumption that $t \geq 1$,
\[
\left \lvert
\mathbb{E}[U^{\scriptscriptstyle\rm h} - U] - \sum_{x} B(h^{\scriptscriptstyle\rm L},r_x) \right \rvert \leq  4\mathbb{E}[t^L] + \frac{2R}{R-n}.
\]
Thus it suffices to bound $ \sum_{x} B(h^{\scriptscriptstyle\rm L},r_x) $.
For every $x$, using \prettyref{eq:Btail} and applying \prettyref{lem:hyper-equal} with $y = \frac{m}{R-n}$ and $k=r_x$, we obtain
\begin{align*}
B(h^{\scriptscriptstyle\rm L},r_x)
& = -
\left(1-\frac{n+m}{R} \right)^{r_x}\int^{\frac{m}{R-n}}_0  \mathbb{E} \left[{r_x-1 \choose L} (-s)^L \right]  r_x (1-s)^{-r_x-1} d s.
\end{align*}
Since $0 \leq \frac{m}{R-n} \leq1$, letting
$K = \max_{0 \leq s \leq 1} \big|  \mathbb{E} \big[{r_x-1 \choose L} (-s )^L  \big] \big|$,
we have
\begin{align*}
 |B(h^{\scriptscriptstyle\rm L},r_x)|
& \leq \left(1-\frac{n+m}{R} \right)^{r_x} K \int^{\frac{m}{R-n}}_0 r_x (1-s)^{-r_x-1} d s. \\
& = K
\left(\left(1-\frac{n}{R} \right)^{r_x} - \left(1-\frac{n+m}{R} \right)^{r_x}  \right)
\leq  K
 \left(1-\frac{n}{R} \right)^{r_x-1} \frac{mr_x}{R}, \label{eq:claim}
\end{align*}
where the last inequality follows from the convexity of $x \mapsto (1-x)^{r_x}$.
Summing over all symbols $x$ results in the lemma.
\end{proof}
Combining \prettyref{lem:hyper_bias} and \prettyref{lem:hyper_var} gives the following NMSE bound:
\begin{Theorem}
\label{thm:hyper}
Under the assumption of \prettyref{lem:hyper_bias},
\[
\mathbb{E}[(U^{\scriptscriptstyle\rm L} - U)^2] \leq 12(n+1){\mathbb{E}}^2[t^L] + 6n + 3m +  \frac{12R^2}{(R-n)^2} + 3m^2 \max_{1 \geq\alpha > 0}  \left \lvert  \mathbb{E} \left[{r_x-1 \choose L} (-\alpha)^L \right] \right \rvert^2.
\]
\end{Theorem}
As before, we can choose various smoothing distribution and obtain upper bounds on the mean squared error.
\begin{Corollary}
\label{cor:hyper}
If $L\sim\mathrm{poi}(r)$ and $R-n \geq m \geq n$, then
\[
\mathbb{E}[(U^{\scriptscriptstyle\rm L} - U)^2] \le 12(n+1)e^{2r(t-1)} +  3m^2 e^{-r} + 9m + 48.
\]
Furthermore, if $r = \frac{1}{2t-1} \cdot \log (nt^2)$,
\[
\cE_{n,t}(U^{\scriptscriptstyle\rm L})
 \leq \frac{27}{(nt^2)^{\frac{1}{2t-1}}} + \frac{9nt + 48}{(nt)^2}.
\]
\end{Corollary}
\begin{proof}
For $L\sim\mathrm{poi}(r)$, $\mathbb{E}[t^L] = e^{r(t-1)}$ and
\[
\max_{0 \leq \alpha \leq 1}  \left \lvert  \mathbb{E} \left[{r_x-1 \choose L} (-\alpha)^L \right] \right \rvert
= e^{-r} \max_{0 \leq \alpha \leq 1} |L_{r_x-1} \left(\alpha r \right)| \leq e^{-r/2},
\]
where $L_{r_x-1}$ is the Laguerre polynomial of degree $r_x-1$ defined in \prettyref{eq:lag} and the last equality follows the bound \prettyref{eq:lagbound}.
Furthermore, $R/(R-n) = 1 + n/(R-n) \leq 1 + n/m \leq 2$ and $n \leq m$, and hence the first part of the lemma.
The second part follows by substituting the value of $r$.
\end{proof}
\section{Lower bounds}
\label{sec:lower}
Under the multinomial model (i.i.d.~sampling), we lower bound the risk
$\cE_{n,t}(U^{\scriptscriptstyle\rm E})$ for any estimator $U^{\scriptscriptstyle\rm E}$ using the support size estimation
lower bound in~\cite{WY14b}. Since the lower bound in~\cite{WY14b}
also holds for the Poisson model, so does our lower bound.

Recall that for a discrete distribution $p$, $S(p) = \sum_x \mathbbm{1}_{p_x > 0}$ denotes its support size.
It is shown that given $n$ i.i.d.~samples drawn from a distribution $p$ whose minimum non-zero mass $p^+_{\min}$ is at least $1/k$, the minimax mean-square error for estimating $S(p)$ satisfies
\begin{equation}
\min_{\hat{S}} \max_{p : p^+_{\min} \geq 1/k} \mathbb{E}[(\hat{S} - S(p))^2] \geq c' k^2 \cdot \exp \left(- c \max\left(\sqrt{\frac{n\log k}{k}}, \frac{n}{k} \right) \right).
	\label{eq:supplb}
\end{equation}
where $c,c'$ are universal positive constants with $c > 1$. We prove Theorem~\ref{thm:res_lb} under the multinomial model with $c$ being the universal constant from \prettyref{eq:supplb}.

Suppose we have an estimator $\hat{U}$ for $U$ that can accurately predict the number of new symbols arising in the next $m$ samples,
we can then produce an estimator for the support size by adding the number of symbols observed, $\prev{+}$, in the current $n$ samples, namely,
\begin{equation}
\hat{S} = \hat{U} +  \prev{+}.
	\label{eq:hatS}
\end{equation}
Note that $U = \sum_x \mathbbm{1}_{N_x=0} \mathbbm{1}_{N'_x>0}$.
When $m=\infty$,
$U$ is the total number of unseen symbols and we have $S(p) = U + \prev{+}$. Consequently, if $\hat{U}$ can foresee too far into the future (\emph{i.e., } for too large an $m$), then \prettyref{eq:hatS} will constitute a support size estimator that is too good to be true.

Combining \prettyref{thm:res_lb} with the positive result (\prettyref{cor:poi_tail} or \ref{cor:bin_tail})
yields the following characterization of the minimax risk:
\begin{Corollary}
\label{cor:rate}
For all $t \geq c$, we have
\[
\inf_{U^{\scriptscriptstyle\rm E}} \cE_{n,t}(U^{\scriptscriptstyle\rm E}) = \exp\pth{ - \Theta\pth{\max\sth{\frac{\log n}{t}, 1}}}
\]
Consequently, as $n\to\infty$, the minimax risk $\inf_{U^{\scriptscriptstyle\rm E}} \cE_{n,t}(U^{\scriptscriptstyle\rm E}) \to 0$ if and only if $t = o(\log n)$.
\end{Corollary}

\begin{proof}[Proof of \prettyref{thm:res_lb}]
Recall that $m=nt$.
Let $\hat U$ be an arbitrary estimator for $U$.
For the support size estimator $\hat S = \hat{U} + \prev{+}$ defined in \prettyref{eq:hatS}, it must obey the lower bound \prettyref{eq:supplb}. Hence there exists some $p$ satisfying $p_{\min}^+ \geq 1/k$,
such that
\begin{equation}
\mathbb{E}[(S(p)- {\hat{S}})^2] \geq c' k^2 \cdot \exp \left(- c\max\left(\sqrt{\frac{n\log k}{k}}, \frac{n}{k} \right) \right).
	\label{eq:supplb1}
\end{equation}
Let $S=S(p)$ denote the support size, which is at most $k$.
Let $\tilde{U} \stackrel{\mathrm{def}}{=} \mathbb{E}_{X^{n+m}_{n+1}}[U]$ be the expectation of $U$ over the unseen samples $X^{n+m}_{n+1}$ conditioned on the available samples $X_1^n$. Then
$\tilde U = \sum_x \mathbbm{1}_{\mul{x} = 0} \left(1 - (1-p_x)^{nt} \right)$.
Since the estimator $\hat{U}$ is independent of $X^{n+m}_{n+1}$, by convexity,
\begin{equation}
\mathbb{E}_{X^{n+m}_1}
[(U-\hat{U})^2] \geq
\mathbb{E}_{X^{n}_{1}}[(\mathbb{E}_{X^{n+m}_{n+1}}[U-\hat{U}])^2] = \mathbb{E}[(\tilde{U}-\hat{U})^2].
\label{eq:UUtilde}
\end{equation}
Notice that with probability one,
\begin{equation}
|S - \tilde{U} - \prev{+} |\leq  Se^{-nt/k}	 \leq k e^{-nt/k},
	\label{eq:Snewt}
\end{equation}
which follows from
\[
\tilde{U} +  \prev{+}
=  \sum_{x: p_x > 0} \mathbbm{1}_{\mul{x} = 0} \left(1 - (1-p_x)^{nt} \right) +  \mathbbm{1}_{\mul{x} > 0}  \leq  S,
\]
and, on the other hand,
\begin{align*}
\tilde{U} +  \prev{+}
&=  \sum_{x: p_x \geq 1/k} \mathbbm{1}_{\mul{x} = 0} \left(1 - (1-p_x)^{nt} \right) +  \mathbbm{1}_{\mul{x} > 0}  \\
& \geq  \sum_x \mathbbm{1}_{\mul{x} = 0} \left(1 - (1-1/k)^{nt} \right) +  \mathbbm{1}_{\mul{x} > 0}  \geq S(1 -  (1-1/k)^{nt}) \geq S(1-e^{-nt/k}).
\end{align*}
Expanding the left hand side of \prettyref{eq:supplb1},
\begin{align*}
\mathbb{E}[(S - {\hat{S}})^2]
& = \mathbb{E}\qth{\pth{  S - \tilde{U} - \prev{+} + \tilde{U} - \hat{U} }^2} \leq  2\mathbb{E}[(S - \tilde{U} - \prev{+} )^2] + 2\mathbb{E}[ (\tilde{U} - \hat{U}))^2] \\
& \overset{\prettyref{eq:Snewt}}{\leq} 2 k^2 e^{-2nt/k}  + 2\mathbb{E}[ (\tilde{U} - \hat{U}))^2] \overset{\prettyref{eq:UUtilde}}{\leq} 2 k^2 e^{-2nt/k}  + 2\mathbb{E}[ (U - \hat{U}))^2] \\
\end{align*}
Let
\[
k = \min\sth{\frac{nt^2}{c^2 \log \frac{nt^2}{c^2}}, \frac{nt}{\log \frac{4}{c'}}},
\]
which ensures that
\begin{equation}
c' k^2 \cdot \exp \pth{- c \max\sth{ \sqrt{\frac{n\log k}{k}}, \frac{n}{k} } } \geq 4 k^2 e^{-2nt/k}.
	\label{eq:lb1}
\end{equation}
 Then
\[
\mathbb{E}[ (U - \hat{U})^2] \geq k^2 e^{-2nt/k},
\]
establishes the following lower bound  with
$\alpha \stackrel{\mathrm{def}}{=} \frac{c'^2}{4 \log^2 (4/c')}$
and $\beta \stackrel{\mathrm{def}}{=} c^2$:
\begin{equation*}
\min_{E} \cE_{n,t}(U^{\scriptscriptstyle\rm E})
\geq \min\sth{\alpha, \frac{4t^2}{\beta^2 \log^2 \frac{nt^2}{\beta}} \pth{\frac{\beta}{nt^2}}^{2\beta/t}}.
\end{equation*}
To verify \prettyref{eq:lb1}, since $t \geq c$ by assumption, we have
$\exp(\frac{2tn}{k} - \frac{cn}{k}) \geq \exp(\frac{nt}{k}) \geq \frac{4}{c'}$.  Similarly,
since $k \log k \leq \frac{nt^2}{c^2}$ by definition, we have
$\frac{2nt}{k} \geq 2 c' \sqrt{\frac{n \log k}{k}}$ and hence $\exp\big(\frac{2tn}{k} - c\sqrt{\frac{n \log k}{k}}\big) \geq \exp(\frac{nt}{k}) \geq \frac{4}{c'}$, completing the proof of \prettyref{eq:lb1}.

Thus we have shown that there exist universal
positive constants $\alpha,\beta$ such that
\begin{equation*}
\min_{E} \cE_{n,t}(U^{\scriptscriptstyle\rm E})
\geq \min\sth{\alpha, \frac{4t^2}{\beta^2 \log^2 \frac{nt^2}{\beta}} \pth{\frac{\beta}{nt^2}}^{2\beta/t}}.
\end{equation*}
Let $y = \left( \frac{nt^2}{\beta}\right)^{2\beta/t}$, then
\[
\min_{E} \cE_{n,t}(U^{\scriptscriptstyle\rm E}) \geq \min\sth{\alpha, 16 \frac{1}{y\log^2 y}}.
\]
Since $y > 1$, $y^3 \geq y \log^2 y$ and hence for some constants $c_1, c_2 > 0$,
\[
\min_{E} \cE_{n,t}(U^{\scriptscriptstyle\rm E}) \geq \min\sth{\alpha, 16 \frac{1}{y^3}}
\geq \min \sth{\alpha, \left( \frac{\beta}{nt^2}\right)^{6\beta/t}}
\geq c_1 \min \sth{1,  \left( \frac{1}{n}\right)^{c_2/t}}
\geq \frac{c_1}{n^{c_2/t}}. \qedhere
\]
\end{proof}
\section{Experiments}
\label{sec:experiments}
We demonstrate the efficacy of our estimators by comparing their
performance with that of several state-of-the-art support-size estimators
currently used by ecologists:
Chao-Lee estimator~\cite{C84, CL92},
Abundance Coverage Estimator (ACE)~\cite{C05},
and the jackknife estimator~\cite{SV84},
combined with the Shen-Chao-Lin unseen-species estimator~\cite{SCL03}.
We consider various natural synthetic distributions and established datasets.
Starting with the former, \arxiv{Figure~\ref{fig:synthetic}}{(Fig.\ 2)} shows the \emph{species discovery curve},
the prediction of $U$ as a function of $t$ of several
predictors for various distributions.
The true value is shown in black, and the other estimators are
color coded, with the solid line representing their mean estimate,
and the shaded area corresponding to one standard deviation.
Note that the Chao-Lee and ACE estimators are designed specifically
for uniform distributions, hence in \arxiv{Figure~\ref{fig:uniform}}{(Fig.\ 2a)} they coincide with the
true value, but for all other distributions, our proposed smoothed Good-Toulmin
estimators outperform the existing ones.
\begin{figure}
\centering
\subfigure[Uniform]{\label{fig:uniform}\arxiv{\includegraphics[width=60mm]{syn_uniform.pdf}}{\includegraphics[width=60mm]{syn_uniform.pdf}}}
\subfigure[Two steps: $\frac{1}{2k} \times \frac{k}{2} \cup \frac{3}{2k} \times \frac{k}{2}$  ]{\label{fig:step}\arxiv{{\includegraphics[width=60mm]{syn_step.pdf}}}{\includegraphics[width=60mm]{syn_step.pdf}}}
\subfigure[Zipf-$1$: $p_i \propto \frac{1}{i}$]{\label{fig:zipf1}\arxiv{\includegraphics[width=60mm]{syn_zipf1.pdf}}{\includegraphics[width=60mm]{syn_zipf1.pdf}}}
\subfigure[Zipf-$1.5$: $p_i \propto \frac{1}{i^{1.5}}$]{\label{fig:zipf32}\arxiv{\includegraphics[width=60mm]{syn_zipf32.pdf}}{\includegraphics[width=60mm]{syn_zipf32.pdf}}}
\subfigure[Dirichlet-$1$ prior]{\label{fig:dir1}\arxiv{\includegraphics[width=60mm]{syn_dir1.pdf}}{\includegraphics[width=60mm]{syn_dir1.pdf}}}
\subfigure[Dirichlet-$1/2$ prior]{\label{fig:dir12}\arxiv{\includegraphics[width=60mm]{syn_dir12.pdf}}{\includegraphics[width=60mm]{syn_dir12.pdf}}}
\begin{tabular}{| c                                                                                | l l | l l |}
\hline
True value  & Previous & &  Proposed & \\ \hline
 {\color{black} \rule[0.08cm]{0.05\textwidth}{1pt}}    & Chao-Lee  & {\color{myred} \rule[0.08cm]{0.05\textwidth}{1pt}} & Poisson smoothing  & {\color{myblue}\rule[0.08cm]{0.05\textwidth}{1pt}}  \\
       & ACE   & {\color{mycyan} \rule[0.08cm]{0.05\textwidth}{1pt}} & Binomial smoothing $q=1/(t+1)$ &  {\color{myorange} \rule[0.08cm]{0.05\textwidth}{1pt}}\\
       & Jackknife  & {\color{myyellow} \rule[0.08cm]{0.05\textwidth}{1pt}}  & Binomial smoothing $q = 2/(t+2)$ & {\color{mygreen} \rule[0.08cm]{0.05\textwidth}{1pt}} \\
\hline
\end{tabular}
\caption{Comparisons of the estimated number of unseen species
as a function of $t$. All experiments have distribution support size $10^6$,
$n=5 \cdot 10^5$, and are averaged over $100$ iterations.}
\label{fig:synthetic}
\end{figure}
Of the proposed estimators, the binomial-smoothing estimator
with parameter $q = \frac{2}{2+t}$ has a stronger theoretical
guarantee and performs slightly better than the others.
Hence when considering real data we plot only its performance
and compare it with the other state-of-the art estimators.
We test the estimators on three real datasets taken from various
scientific applications where the samples size $n$ ranges from few
hundreds to a million. For all these date sets, our estimator
outperforms the existing procedures

\begin{figure}
\centering
\subfigure[Hamlet random]{\label{fig:hamlet_random}\arxiv{\includegraphics[scale=0.35]{Hamlet_data-crop.pdf}}{\includegraphics[scale=0.35]{Hamlet_data-crop.pdf}}}
\subfigure[Hamlet consecutive]{\label{fig:hamlet_consecutive}\arxiv{\includegraphics[scale=0.35]{Hamlet_consecutive-crop.pdf}}{\includegraphics[scale=0.35]{Hamlet_consecutive-crop.pdf}}}
\subfigure[SLOTUs]{\label{fig:haplo}\arxiv{\includegraphics[scale=0.35]{Haplo_data-crop.pdf}}{\includegraphics[scale=0.35]{Haplo_data-crop.pdf}}}
\subfigure[Last names]{\label{fig:name}\arxiv{\includegraphics[scale=0.35]{Name_data-crop.pdf}}{\includegraphics[scale=0.35]{Name_data-crop.pdf}}}
\begin{tabular}{| c                                                                                | l l | l l |}
\hline
True value  & Previous & &  Proposed & \\ \hline
 {\color{black} \rule[0.08cm]{0.05\textwidth}{1pt}}
 & Chao-Lee
& {\color{myred} \rule[0.08cm]{0.05\textwidth}{1pt}}
& Binomial smoothing $q = 2/(t+2)$ & {\color{mygreen} \rule[0.08cm]{0.05\textwidth}{1pt}} \\
       & ACE   & {\color{mycyan} \rule[0.08cm]{0.05\textwidth}{1pt}}   & & \\
       & Jackknife  & {\color{myyellow} \rule[0.08cm]{0.05\textwidth}{1pt}}  &  & \\
& Empirical & {\color{mybrown} \rule[0.08cm]{0.05\textwidth}{1pt}}   & & \\
\hline
\end{tabular}
\caption{Estimates for number of: $(a)$ distinct words in Hamlet with
random sampling $(b)$ distinct words in Hamlet with consecutive
sampling $(c)$ SLOTUs on human skin $(d)$ last names.}
\label{fig:hamlet}
\end{figure}

\arxiv{Figure~\ref{fig:hamlet_random}}{(Fig.\ 3a)} shows the first real-data experiment, predicting
vocabulary size based on partial text.
Shakespeare's play \emph{Hamlet} consists of $n_{\text{total}}=31999$ words,
of which $4804$ are distinct.
We randomly select $n$ of the $n_{\text{total}}$ words without replacement,
predict the number of unseen words in $n_{\text{total}}-n$ new ones,
and add it to those observed.
The results shown are averaged over $100$ trials.
Observe that the new estimator outperforms existing ones and
that as little as $20\%$ of the data
already yields an accurate estimate of the total number of
distinct words.
\arxiv{Figure~\ref{fig:hamlet_consecutive}}{(Fig.\ 3b)} repeats the experiment but instead of random sampling,
uses the first $n$ consecutive words, with similar conclusions.

\arxiv{Figure~\ref{fig:haplo}}{(Fig.\ 3c)} estimates the number of bacterial species on the human skin.
\cite{GCB07} considered forearm skin biota of six subjects.
They identified $n_{\text{total}}=1221$ clones consisting
of $182$ different species-level operational taxonomic units (SLOTUs).
As before, we select $n$ out of the $n_{\text{total}}$ clones without replacement and
predict the number of distinct SLOTUs found.
Again the estimates are more accurate than those of existing
estimators and are reasonably accurate already with $20\%$ of the data.

Finally, \arxiv{Figure~\ref{fig:name}}{(Fig.\ 3d)} considers the 2000 United States
 Census~\cite{Census00}, which lists all U.S. last names
corresponding to at least 100 individuals.
With these many repetitions, even just a small fraction of the data
will cover all names, hence we
first subsampled the data $n_{\text{total}}=10^6$ and obtained a list of
100328 distinct last names.
As before we estimate for this number
using $n$ randomly chosen names, again with similar conclusions.
\section*{Acknowledgments}
\label{sec:acknowledgements}
This work was partially completed while the authors were visiting the
Simons Institute for the Theory of Computing at UC Berkeley, whose
support is gratefully acknowledged.  We thank Dimitris Achlioptas,
David
Tse, Chi-Hong Tseng, and Jinye Zhang for helpful discussions and comments.

redpajama_set_name: RedPajamaArXiv