11affiliationtext: Department of Statistical Science, Duke University, Durham,
North Carolina 27708, U.S.A.22affiliationtext: Department of Economics,
Management and Statistics, University of Milano–Bicocca, 20126 Milano,
Italy33affiliationtext: Organismal and Evolutionary Biology Research
Programme, University of Helsinki, Helsinki, Finland

# Bayesian nonparametric modelling of sequential discoveries

Alessandro Zito  Tommaso Rigon  Otso Ovaskainen  David Dunson

###### Abstract

We aim at modelling the appearance of distinct tags in a sequence of labelled
objects. Common examples of this type of data include words in a corpus or
distinct species in a sample. These sequential discoveries are often
summarised via accumulation curves, which count the number of distinct
entities observed in an increasingly large set of objects. We propose a novel
Bayesian nonparametric method for species sampling modelling by directly
specifying the probability of a new discovery, therefore allowing for flexible
specifications. The asymptotic behavior and finite sample properties of such
an approach are extensively studied. Interestingly, our enlarged class of
sequential processes includes highly tractable special cases. We present a
subclass of models characterized by appealing theoretical and computational
properties. Moreover, due to strong connections with logistic regression
models, the latter subclass can naturally account for covariates. We finally
test our proposal on both synthetic and real data, with special emphasis on a
large fungal biodiversity study in Finland.

##  1 Introduction

Our goal is to develop a flexible procedure for modelling the appearance of
previously unobserved objects in a sequence. The sequential recording of
distinct entities can be represented through an accumulation curve, namely the
cumulative number of distinct entities $K_{n}$ within a collection of $n$
objects. These entities can be of various nature, including biological species
(Good, 1953; Good and Toulmin, 1956), words (Efron and Thisted, 1976; Thisted
and Efron, 1987), genes (Ionita-Laza et al., 2009), bacteria (Hughes et al.,
2001; Gao et al., 2007) and cell types (Camerlenghi et al., 2020). The
analysis of accumulation curves has a rich history in statistics, as testified
by the early contributions of Fisher et al. (1943), Good (1953), and Good and
Toulmin (1956). We refer to Bunge and Fitzpatrick (1993) for a historical
account on the topic. Several nonparametric approaches have been developed in
more recent years, often aiming at i) predicting the number of unseen entities
(e.g. Shen et al., 2003), or ii) estimating the probability of a new discovery
(e.g. Chao and Shen, 2004; Mao, 2004; Favaro et al., 2012).

Our work builds on Bayesian nonparametric methods, whose development has been
spurred by the seminal paper of Ferguson (1973) on the Dirichlet process. In
our motivating application, we aim to assess how many of the species present
in a sample are missed when a given number of dna barcode sequences are
obtained through high-throughput sequencing. Let $(X_{n})_{n\geq 1}$ be a
sequence of objects, such as fungal dna sequences in a single biological soil
or air sample (Abrego et al., 2020), taking values in $\mathds{X}$, which is
the space of fungal species in our case. Among the first $n$ observed objects
$X_{1},\dots,X_{n}$, there will be $K_{n}\leq n$ distinct entities, or
species, representing the $n$th value of the accumulation curve. The values
$(X_{n})_{n\geq 1}$ are randomly generated in a sequential manner, so that the
tag $X_{n+1}$ is either new or equal to one of the previously observed
objects. For instance, in the Dirichlet process case, the sequential
allocation mechanism for any $n\geq 1$ proceeds as follows:

$(X_{n+1}\mid X_{1},\dots,X_{n})=\begin{cases}\textrm{``new''},&\textrm{with
probability}\quad\alpha/(\alpha+n),\\\ X_{i},&\textrm{with probability}\quad
1/(\alpha+n),\quad(i=1,\ldots,n),\end{cases}$ (1)

where $\alpha>0$ controls the rate of new discoveries; see also Blackwell and
MacQueen (1973).

The predictive scheme in (1) is restrictive in depending on a single parameter
and in inducing a logarithmic growth for the accumulation curve
$(K_{n})_{n\geq 1}$. These limitations motivated the development of more
general random processes that allow for polynomial growth rates. These include
the two parameter Poisson–Dirichlet process of Perman et al. (1992), often
called the Pitman–Yor process when the number of species is assumed to be
infinite or the Dirichlet-multinomial process in the finite case (Pitman and
Yor, 1997), and the general classes of Gibbs-type priors (Gnedin and Pitman,
2005) and species sampling models (Pitman, 1996). The derivation of Bayesian
nonparametric estimators for accumulation curves, under general Gibbs-type
priors and the Pitman–Yor process, is due to Lijoi et al. (2007) and Favaro et
al. (2009), respectively.

Unfortunately, tractable generalizations of (1) such as the Pitman–Yor process
are too restrictive for many real-world scenarios. This is evident from Figure
1, which shows in- and out-of-sample performance in estimating the number of
distinct fungi species in a given number of fungal dna-barcode sequences.
‘Species’ are defined in this article based on genetic sequences being
sufficiently distinct, but the terminology used by ecologists is ‘operational
taxonomic units’ as determining species requires additional verification. The
Dirichlet process fails badly in sample, while the Pitman–Yor has good in-
sample fit but poor out-of-sample predictive accuracy. This is not surprising,
as the Pitman–Yor process depends on only two parameters and assumes that
$K_{n}\rightarrow\infty$ almost surely as $n\rightarrow\infty$. As there are
finitely many fungi species, $K_{n}$ should more realistically converge to a
finite constant. The Dirichlet-multinomial process allows finite
$\lim_{n\rightarrow\infty}K_{n}=K_{\infty}$ but the trajectory has similar
lack of fit as the Dirichlet process.

Potentially one could use a species sampling model that is more flexible than
the Pitman–Yor, while also allowing finite $K_{\infty}$; recent examples
include Camerlenghi et al. (2018); Lijoi et al. (2020). However, such
specifications involve cumbersome combinatorial structures in the sampling
mechanism, effectively preventing their application in the types of large
datasets that are now routinely collected in our motivating application areas.
For example, in fungi biodiversity studies, it is now common to sequence
millions of dna barcodes from 10,000s of species (e.g. Ovaskainen et al.,
2020).

Figure 1: Empirical and estimated accumulation curve in one air fungal DNA-
barcoding sample obtained from Finland. White dots indicate the observed
values. Left panel: the vertical line is the train-test set cutoff, set to
$1/3$ of the total number of genetic sequences. The parameters of the
Dirichlet, the Pitman–Yor and the Dirichlet-multinomial processes are based on
the training set. Right panel: the curves are instead estimated using the full
dataset

We address the above limitations through a novel modelling framework, which is
highly flexible, analytically tractable, and computationally efficient. The
key distinction compared to species sampling models, such as (1), is that we
directly specify a model for the accumulation curve $(K_{n})_{n\geq 1}$,
whereas the tags $(X_{n})_{n\geq 1}$ are regarded as nuisance parameters.
Specifically, we consider a collection of Bernoulli random variables
$(D_{n})_{n\geq 1}$ representing whether at the $(n+1)$th step a new entity
has been discovered or not, namely

$\textrm{pr}(D_{n+1}=1)=\textrm{pr}(X_{n+1}=\text{``new"}\mid
X_{1},\dots,X_{n}),\qquad n\geq 1,$

having set $D_{1}=1$. The accumulation curve is obtained by summing over these
binary indicators: $K_{n}=\sum_{i=1}^{n}D_{i},n\geq 1.$ Differently from
general species sampling models, in our framework, the Bernoulli indicators
$(D_{n})_{n\geq 1}$ are assumed to be _independent_ , albeit not identically
distributed. Hence, we aim at developing suitable formulations for the
probabilities $(\pi_{n})_{n\geq 1}$, with $\pi_{n}=\mathrm{pr}(D_{n}=1)$, for
any $n\geq 1$. It is natural to require these probabilities to be decreasing
over $n$, so that the discovery of a new entity is increasingly difficult the
more data we collect. Moreover, $\pi_{1}=\mathrm{pr}(D_{1}=1)=1$, since the
first entity of the sequence is necessarily new. Both requirements are
satisfied by the Dirichlet process, where $\pi_{n}=\alpha/(\alpha+n-1)$. We
propose a general strategy for the specification of $(\pi_{n})_{n\geq 1}$,
relying on the notion of survival functions, and study the impact of specific
choices on the asymptotic behavior of $K_{n}$.

A specific subclass of our framework is particularly appealing in terms of
analytic and computational simplicity, due to connections with logistic
regression. This subclass includes the Dirichlet process and naturally leads
to covariate-dependent extensions. Existing covariate-dependent species
sampling models are typically complex to implement; refer to Quintana et al.
(2020) for a recent overview. In contrast, our approach simply involves
implementing a logistic regression with certain constraints on the parameters.
We illustrate the flexibility and computational tractability through
application to fungi biodiversity data collected at different sampling sites
in Finland under different ecological conditions (Abrego et al., 2020).

##  2 A general modelling framework for accumulation curves

###  2.1 Background on species sampling models

In this Section we review key concepts about species sampling models that will
be used throughout the paper. For a broader overview, refer to Pitman (1996)
and De Blasi et al. (2015).

Let $(X_{n})_{n\geq 1}$ be a sequence of objects. Given the discrete nature of
the data, there will be ties among $X_{1},\dots,X_{n}$, comprising a total of
$K_{n}=k$ distinct entities $X_{1}^{*},\dots,X_{k}^{*}$, having frequencies
$n_{1},\dots,n_{k}$, with $\sum_{j=1}^{k}n_{j}=n$. Species sampling models
generalize the sequential allocation of the Dirichlet process in (1), so that
for any $n\geq 1$,

$(X_{n+1}\mid X_{1},\dots,X_{n})=\begin{cases}\textrm{``new''},&\textrm{with
probability}\quad p_{k+1}(n_{1},\dots,n_{k}),\\\ X_{j}^{*},&\textrm{with
probability}\quad p_{j}(n_{1},\dots,n_{k}),\quad(j=1,\ldots,k),\end{cases}$
(2)

for suitable probabilities $\sum_{j=1}^{k+1}p_{j}(n_{1},\dots,n_{k})=1$ and
with $X_{1}=\text{``new''}$. The discovery probabilities
$\textrm{pr}(X_{n+1}=\text{``new''}\mid X_{1},\dots,X_{n})$ depend on the
previous values only through $k$ and/or the frequencies $n_{1},\dots,n_{k}$.
Equation (2) only leads to a valid species sampling model if the resulting law
of $(X_{n})_{n\geq 1}$ is _exchangeable_ ; this is not automatic as discussed
in Lee et al. (2013). In the Pitman–Yor case,
$p_{k+1}(n_{1},\dots,n_{k})=(\alpha+k\sigma)/(\alpha+n)$ and
$p_{j}(n_{1},\dots,n_{k})=(n_{j}-\sigma)/(\alpha+n)$, for $j=1,\dots,k$ and
with $\sigma\in[0,1)$ and $\alpha>-\sigma$; the Dirichlet process is recovered
with $\sigma=0$. The Dirichlet-multinomial has the same sampling scheme of the
Pitman–Yor, having set $\sigma<0$ and $\alpha=H|\sigma|$, with
$H\in\mathds{N}$ representing the total number of species.

The sequential mechanism in (2) induces a law for the accumulation curve
$(K_{n})_{n\geq 1}$. For the remainder of the Section, we focus on the
Dirichlet process, since it provides a special case of our framework. In this
case, the distribution of $K_{n}$ is available in closed form (Antoniak,
1974),

$\mathrm{pr}(K_{n}=k)=\frac{\alpha^{k}}{(\alpha)_{n}}|s(n,k)|,\qquad\mbox{for
any $n\geq 1$},$ (3)

where $|s(n,k)|$ is the signless Stirling number of the first kind; see
Charalambides (2005). Moreover, the expectation of (3) is

$E(K_{n})=\sum_{i=1}^{n}\frac{\alpha}{\alpha+i-1},$ (4)

which provides the prior mean for the accumulation curve. One may also be
interested in the posterior distribution of $K_{m}^{(n)}$, the number of new
entities in a future sample of size $m$ conditioning on training data
$X_{1},\dots,X_{n}$. Under a Dirichlet process this distribution does not
depend on either the past values $X_{1},\dots,X_{n}$ or the observed number of
distinct values $K_{n}=k$; see Lijoi et al. (2007). Consequently, we obtain

$E(K_{m}^{(n)}\mid
X_{1},\dots,X_{n})=\sum_{i=1}^{m}\frac{\alpha}{\alpha+n+i-1}.$ (5)

The Dirichlet process is the only species sampling model for which such a
simplification occurs (Lijoi et al., 2007). For example, in the Pitman–Yor
process the posterior distribution of $K_{m}^{(n)}$ depends on the observed
number of distinct values $K_{n}=k$.

###  2.2 The model

In species sampling models, the distribution of the accumulation curve
$(K_{n})_{n\geq 1}$ is essentially a byproduct of the specification for the
values $(X_{n})_{n\geq 1}$. Instead, we propose a more direct formulation for
$(K_{n})_{n\geq 1}$ which avoids modelling of the sequence $(X_{n})_{n\geq
1}$.

Let $(D_{n})_{n\geq 1}$ be a collection of _independent_ binary indicators,
denoting the discoveries, with probabilities $(\pi_{n})_{n\geq 1}$. Moreover,
let $K_{n}=\sum_{i=1}^{n}D_{i}$ for any $n\geq 1$ be the accumulation curve.
We have $D_{n}=K_{n}-K_{n-1}$ for any $n\geq 2$ with $D_{1}=1$. Hence, the
discoveries $(D_{n})_{n\geq 1}$ and the accumulation curve $(K_{n})_{n\geq 1}$
carry the same information, having a one-to-one relationship. The resulting
distribution of $K_{n}$ is Poisson-binomial with parameters
$\pi_{1},\dots,\pi_{n}$, which we denote as
$K_{n}\sim\textsc{pb}(\pi_{1},\dots,\pi_{n})$. As previously discussed, the
probabilities $(\pi_{n})_{n\geq 1}$ must satisfy $\pi_{n}>\pi_{n+1}$ for every
$n\geq 1$ with $\pi_{1}=1$. In addition, we require
$\lim_{n\rightarrow\infty}\pi_{n}=0$, meaning that the probability of making a
new discovery should eventually approach zero. A general strategy for
constructing a set of probabilities satisfying these requirements is described
as follows.

######  Definition 1.

Let $T$ be a random variable on $(0,\infty)$ with strictly increasing
cumulative distribution function $F(t;\theta)$ indexed by
$\theta\in\Theta\subseteq\mathds{R}^{p}$. Moreover, let
$S(t;\theta)=1-F(t;\theta)$ be its survival function. The set of probabilities
$(\pi_{n})_{n\geq 1}$ are said to be directed by $S(t;\theta)$ if

$\pi_{n}=\textup{pr}(T_{n}>n-1)=S(n-1;\theta),\qquad\mbox{for any $n\geq 1$},$
(6)

where $(T_{n})_{n\geq 1}$ are independent and identically distributed random
variables following $F(t;\theta)$.

It is easy to check that a set of probabilities $(\pi_{n})_{n\geq 1}$ directed
by $S(t;\theta)$ satisfies the aforementioned requirements. Indeed, one has
that $\pi_{1}=S(0;\theta)=1$ for any $\theta\in\Theta$, since $T$ is supported
on $(0,\infty)$. Moreover, $\pi_{n}=S(n-1;\theta)>S(n;\theta)=\pi_{n+1}$,
because by assumption $S(t;\theta)$ is strictly decreasing. Furthermore, one
has that
$\lim_{n\rightarrow\infty}\pi_{n}=\lim_{n\rightarrow\infty}S(n-1;\theta)=0$,
as desired, since $S(t;\theta)$ is a survival function. Each binary random
variable $D_{n}$ may be represented as $D_{n}=\mathds{1}(T_{n}>n-1)$, with
$\mathds{1}(\cdot)$ denoting the indicator function.

If the probabilities $(\pi_{n})_{n\geq 1}$ are directed by $S(t;\theta)$, then
inferential statements about the parameter vector $\theta\in\Theta$ can be
based on the likelihood function $\mathscr{L}(\theta\mid D_{1},\dots,D_{n})$
or, equivalently, on $\mathscr{L}(\theta\mid K_{1},\dots,K_{n})$. The former
is readily available as

$\mathscr{L}(\theta\mid
D_{1},\dots,D_{n})\propto\prod_{i=2}^{n}S(i-1;\theta)^{D_{i}}\\{1-S(i-1;\theta)\\}^{1-D_{i}},$
(7)

having excluded the degenerate term $D_{1}=1$. The Dirichlet process
implicitly assumes $S(t;\alpha)=\alpha/(\alpha+t)$ with $\alpha>0$, which is
the survival function of a continuous random variable $T$, as will be shown in
Section 3. Since the Dirichlet process is both a species sampling model and a
member of our general framework, it is interesting to compare the information
contained in the data $X_{1},\dots,X_{n}$ with that carried by the discovery
indicators $D_{1},\dots,D_{n}$. This can be formally achieved by comparing the
likelihood functions $\mathscr{L}(\alpha\mid X_{1},\dots,X_{n})$ and
$\mathscr{L}(\alpha\mid D_{1},\dots,D_{n})$. The former may be obtained
following Antoniak (1974), whereas the latter coincides with (7), having set
$S(t;\alpha)=\alpha/(\alpha+t)$. Before stating the result, let us denote with
$(a)_{n}=a(a+1)\cdots(a+n-1)$ the Pochhammer symbol, for any $a>0$ and $n\geq
1$.

######  Theorem 1.

Let $(X_{n})_{n\geq 1}$ be a sequence of objects directed by a Dirichlet
process as in (1) and let $(D_{n})_{n\geq 1}$ be the associated discovery
indicators. Then for a sample $X_{1},\dots,X_{n}$ with $K_{n}=k$ distinct
values one has

$\mathscr{L}(\alpha\mid D_{1},\dots,D_{n})\propto\mathscr{L}(\alpha\mid
X_{1},\dots,X_{n})\propto\frac{\alpha^{k}}{(\alpha)_{n}}.$

Hence, it is equivalent to base inferences on the Dirichlet process parameter
$\alpha$ on the likelihood (7) for the discovery indicators instead of the
usual likelihood for $X_{1},\dots,X_{n}$. Broadly speaking, this occurs
because $K_{n}=\sum_{i=1}^{n}D_{i}$ is the minimal sufficient statistic for
$\alpha$ in the Dirichlet process; see also Lijoi et al. (2007) for similar
considerations. An implication is that the empirical Bayes estimate of
$\alpha$, obtained by maximizing $\alpha^{k}/(\alpha)_{n}$, coincides with the
maximizer of (7).

######  Remark 1.

The Dirichlet process in (1) is the only species sampling model in which the
discovery indicators $(D_{n})_{n\geq 1}$ are independent (Zabell, 1982; Lee et
al., 2013). Thus, our framework can be regarded as the result of some
generative mechanism such as (2), but the underlying sequence $(X_{n})_{n\geq
1}$ will not necessarily be exchangeable. The main implication is that
inference on the parameter $\theta\in\Theta$ may depend on the order of the
observations, which might be reasonable or not depending on the application.
However, the lack of exchangeability has mild practical implications, as
illustrated in the following.

###  2.3 Smoothing, prediction and posterior representations

In this Section, we present prior and posterior properties of $K_{n}$, which
may be useful for both smoothing and prediction. Supposing $(\pi_{n})_{n\geq
1}$ is directed by $S(t;\theta)$, then
$K_{n}\sim\textsc{pb}\\{1,S(1;\theta),\ldots,S(n-1;\theta)\\}$, with prior
mean and variance equal to

$E(K_{n})=\sum_{i=1}^{n}S(i-1;\theta),\quad\mathrm{var}(K_{n})=\sum_{i=1}^{n}S(i-1;\theta)\\{1-S(i-1;\theta)\\},\qquad
n\geq 1.$

These moment formulas may be useful in choosing the parametric form of
$S(t;\theta)$ and for prior elicitation for $\theta$. In ecology, in-sample
estimation of species accumulation curves is sometimes called rarefaction;
this amounts to smoothing of the $K_{1},\ldots,K_{n}$ values observed in the
training samples. In our framework, $E(K_{n})$ is the sum of the first $n$
discovery probabilities

$E(K_{n})=\sum_{i=1}^{n}\pi_{i}=\sum_{i=1}^{n}\textrm{pr}(D_{i}=1),\qquad
n\geq 1.$

This expectation does not depend on the ordering of the data, at least for any
fixed value of $\theta$.

Suppose we are given a sample of $D_{1},\dots,D_{n}$ discoveries displaying
$K_{n}=k$ distinct entities and that we are interested in predicting future
values of the accumulation curve $K_{n+1},\ldots,K_{n+m}$ or in predicting the
number of new entities within a future sample of size $m$,
$K_{m}^{(n)}=K_{n+m}-K_{n}=\sum_{i=n+1}^{n+m}D_{i}$. The posterior
distribution of $(K_{m}^{(n)}\mid D_{1},\dots,D_{n})$ is available in closed
form, as summarized in the next Proposition, immediately leading to the
distribution of $(K_{n+m}\mid D_{1},\dots,D_{n})=k+K_{m}^{(n)}$.

######  Proposition 1.

Let $(D_{n})_{n\geq 1}$ be a collection of independent discovery indicators
with probabilities $(\pi_{n})_{n\geq 1}$ directed by $S(t;\theta)$. Moreover,
let $K_{m}^{(n)}=K_{n+m}-K_{n}$. Then for any $n\geq 1$ and $m>n$

$(K_{m}^{(n)}\mid
D_{1},\dots,D_{n})\sim\textsc{pb}\\{S(n;\theta),\ldots,S(n+m-1;\theta)\\}.$

Hence, it follows that

$E(K_{m}^{(n)}\mid
D_{1},\dots,D_{n})=\sum_{i=n+1}^{n+m}\textup{pr}(D_{i}=1)=\sum_{j=1}^{m}S(j+n-1;\theta),$

implying that $E(K_{n+m}\mid D_{1},\dots,D_{n})=k+E(K_{m}^{(n)}\mid
D_{1},\dots,D_{n})$.

Within the ecological community, out-of-sample prediction through
$E(K_{n+m}\mid D_{1},\dots,D_{n})$ is sometimes called _extrapolation_ , which
again can be interpreted as a sum of discovery probabilities. Proposition 1
implies that the posterior distribution of $K_{m}^{(n)}$ is conjugate, being a
Poisson-binomial distribution with updated parameters. In addition, the
posterior law of $K_{n+m}$ only depends on $K_{n}=k$, meaning that
extrapolation of the accumulation curve also does not depend on the order of
$D_{1},\dots,D_{n}$ for any fixed value of $\theta$.

###  2.4 Asymptotic behavior of $K_{n}$

The limit of $K_{n}$ as $n\rightarrow\infty$ is often of inferential interest,
representing the random number of entities one would eventually discover.
Depending on the choice of $S(t;\theta)$, two scenarios can occur: i) the
number of distinct entities diverges, as in the Dirichlet process case, so
that $K_{n}\rightarrow\infty$ almost surely as $n\rightarrow\infty$. In this
regime, it is useful to study the growth rate of $K_{n}$. Alternatively, we
could find that ii) the number of distinct species converges to some non-
degenerate random variable $K_{n}\rightarrow K_{\infty}$, almost surely, as
$n\rightarrow\infty$. Within ecology the random variable $K_{\infty}$ is
called the _species richness_.

The asymptotic behaviour of $K_{n}$ is controlled by the structure of the
chosen survival function $S(t;\theta)$. Before stating our first result, let
us define
$E(T)=\int_{0}^{\infty}\mathrm{pr}(T>t)\mathrm{d}t=\int_{0}^{\infty}S(t;\theta)\mathrm{d}t,$
that is, the expectation of the latent variables in Definition 1.

######  Proposition 2.

Let $K_{n}\sim\textsc{pb}\\{1,S(1;\theta),\dots,S(n-1;\theta)\\}$. Then, there
exists a possibly infinite random variable $K_{\infty}$ such that
$\lim_{n\rightarrow\infty}K_{n}\to K_{\infty},$ almost surely, with
$E(K_{\infty})=\sum_{i=0}^{\infty}S(i)$. Moreover,

$E(T)\leq E(K_{\infty})\leq E(T)+1.$ (8)

Equation (8) provides lower and upper bounds for the asymptotic mean, which
can be used to summarize the species richness. Besides, the expected value of
$E(T)$ represents a simple tool to determine whether the accumulation curve
diverges or not, as the following Corollary clarifies.

######  Corollary 1.

Under the conditions of Proposition 2, $K_{\infty}=\infty$ almost surely if
and only if $E(T)=\infty$.

Let us consider the first asymptotic regime, corresponding to the
$K_{\infty}=\infty$ case. In this case, the rate of growth is controlled by
$S(t;\theta)$, as clarified in the following Theorem, which also presents a
central limit approximation.

######  Theorem 2.

Let $K_{n}\sim\textsc{pb}\\{1,S(1;\theta),\dots,S(n-1;\theta)\\}$ and suppose
$K_{\infty}=\infty$ almost surely. Then, as $n\rightarrow\infty$,
$K_{n}/b_{n}\to 1$ almost surely, for
$b_{n}=\int_{1}^{n}S(t-1;\theta)\mathrm{d}t$. In addition,

$\frac{K_{n}-E(K_{n})}{\mathrm{var}(K_{n})^{1/2}}\to N(0,1),\qquad
n\rightarrow\infty,$

in distribution.

Theorem 2 implies that the growth rate of $K_{n}$ corresponds to
$b_{n}=\int_{1}^{n}S(t-1;\theta)\mathrm{d}t$. In the Dirichlet process case,
$b_{n}=\alpha\log{(\alpha+n-1)}-\alpha\log{\alpha}$, corresponding to the
well-known growth rate $\alpha\log{n}$ (Korwar and Hollander, 1973). The
$N(0,1)$ limiting distribution allows one to assess uncertainty in $K_{n}$ for
large $n$.

Consider now the second asymptotic regime, namely the $K_{\infty}<\infty$
case. Although the distribution of $K_{\infty}$ is generally not available in
closed form, the first two moments are well defined.

######  Corollary 2.

Under the conditions of Proposition 2, if $K_{\infty}<\infty$ almost surely,
then $E(K_{\infty})=\sum_{i=1}^{\infty}S(i-1;\theta)<\infty$ and
$\textup{var}(K_{\infty})=\sum_{i=1}^{\infty}S(i-1;\theta)\\{1-S(i-1;\theta)\\}<\infty$.

Hence, a natural estimator for the species richness is $E(K_{\infty})$, which
may be numerically approximated; for instance by truncating the infinite
summation $E(K_{\infty})=\sum_{i=0}^{\infty}S(i)$. Alternatively, one could
exploit equation (8) and consider the arithmetic mean of the bounds, obtaining
the approximation $E(K_{\infty})\approx E(T)+1/2,$ which is often easier to
compute than $E(K_{\infty})$ and is highly accurate when the number of species
is not small. Despite the absence of a central limit theorem in this case,
there exist several approximations for Poisson-binomial distributions, which
may be used for $K_{\infty}$ and $K_{n}$ when $n$ is large (Hong, 2013).

Interestingly, Proposition 1 offers a natural estimator for the posterior
species richness as well, namely $E(K_{\infty}\mid D_{1},\ldots,D_{n})$.
Consider $E(K_{m+n}\mid D_{1},\ldots,D_{n})$ and let $m\to\infty$. Then, it is
straightforward to see that $E(K_{\infty}\mid
D_{1},\ldots,D_{n})=k+E(K_{\infty}^{(n)}\mid D_{1},\ldots,D_{n})$, where

$E(K_{\infty}^{(n)}\mid
D_{1},\ldots,D_{n})=\sum_{j=1}^{\infty}S(j+n-1;\theta).$

Hence, all the properties of $K_{\infty}$ can be naturally extended to the
posterior species richness.

##  3 Logistic models

###  3.1 The log-logistic distribution

The framework in the previous Section requires elicitation of $S(t;\theta)$.
In this Section, we focus on a class of survival functions, which lead to a
generalization of the Dirichlet process, enjoy appealing analytical and
computational properties and result in natural covariate-dependent extensions,
as described in Section 3.3. In particular, we first consider a two parameter
case

$S(t;\alpha,\sigma)=\frac{\alpha}{\alpha+t^{1-\sigma}},\qquad t\geq 0,$ (9)

where $\alpha>0$ and $\sigma<1$. The survival function $S(t;\alpha,\sigma)$
characterizes a two-parameter log-logistic distribution, and therefore we will
write $T\sim\textsc{ll}(\alpha,\sigma)$. Clearly, when $\sigma=0$,
$S(t;\alpha,0)$ reduces to the Dirichlet process case. The parameter $\sigma$
plays a similar role to the discount parameter of the Pitman–Yor process and
general Gibbs-type priors. For any $\sigma<0$, one has

$E(T)=\frac{\alpha^{1/(1-\sigma)}\pi}{(1-\sigma)\sin\\{\pi/(1-\sigma)\\}},$

implying that when $\sigma<0$ the limiting distribution $K_{\infty}<\infty$ is
non-degenerate, thanks to Corollary 1. Conversely, when $0\leq\sigma<1$, one
has that both $K_{\infty}=\infty$ and $E(T)=\infty$. The rate at which this
occurs is logarithmic in the Dirichlet process case in which $\sigma=0$. In
contrast, for $\sigma>0$, one can show that the growth of $K_{n}$ is
polynomial, so that in the notation of Theorem 2 one has
$b_{n}=\int_{1}^{n}S(t;\alpha,\sigma)\mathrm{d}t=\mathcal{O}(n^{\sigma})$.
These considerations reinforce the parallelism with Gibbs-type priors; see
Gnedin and Pitman (2005) and De Blasi et al. (2015) for details.

In the next Section, we describe a three-parameter extension of the log-
logistic distribution and derive combinatorial tools and distributional
properties that also apply to $S(t;\alpha,\sigma)$ in (9).

###  3.2 A three parameter log-logistic distribution

In this Section we extend the log-logistic specification by including an
additional parameter, denoted as $\phi$, which forces $K_{n}$ to converge to a
non-degenerate distribution. This allows us to restrict focus to the second
asymptotic regime. In particular, we let $\theta=(\alpha,\sigma,\phi)$ and

$S(t;\alpha,\sigma,\phi)=\frac{\alpha\phi^{t}}{\alpha\phi^{t}+t^{1-\sigma}},\qquad
t\geq 0,$ (10)

with $\alpha>0$, $\sigma<1$ and $0<\phi\leq 1$. The two parameter
specification is recovered when $\phi=1$. We call the distribution of
$S(t;\alpha,\sigma,\phi)$ a three-parameter log-logistic, written
$T\sim\textsc{ll}(\alpha,\sigma,\phi)$.

######  Proposition 3.

Let $K_{n}\sim\textsc{pb}\\{1,S(1;\theta),\dots,S(n-1;\theta)\\}$, with
$S(t;\theta)$ defined as in Equation (10). Then for any $0<\phi<1$ it holds
that $K_{n}\to K_{\infty}<\infty$ almost surely as $n\to\infty$.

Proposition 3 ensures that for $0<\phi<1$ the species richness is always
finite. For the remainder of the Section, we discuss some combinatorial
properties related to the law of $K_{n}$. While having their own theoretical
relevance, our results facilitate computation of the probability mass function
of $K_{n}$ and draw further parallels with Gibbs-type priors.

######  Definition 2.

Let $\alpha>0$, $\sigma<1$ and $0<\phi\leq 1$. Then for any $n\geq 1$ and
$0\leq k\leq n$ we define $\mathscr{C}_{n,k}(\sigma,\phi)$ as the coefficients
of the polynomial expansion
$\prod_{k=0}^{n-1}(\alpha+k^{1-\sigma}\phi^{-k})=\sum_{k=0}^{n}\alpha^{k}\>\mathscr{C}_{n,k}(\sigma,\phi),$
having set $\mathscr{C}_{0,0}(\sigma,\phi)=1$.

In the special case $\phi=1$ and $\sigma=0$ one recovers the definition of the
signless Stirling numbers of the first kind, namely
$\mathscr{C}_{n,k}(0,1)=|s(n,k)|$; see Charalambides (2005). In addition, the
coefficients $\mathscr{C}_{n,k}(\sigma,\phi)$ can be conveniently computed
through recursive formulas.

######  Theorem 3.

The coefficients $\mathscr{C}_{n,k}(\sigma,\phi)$ of Definition 2 satisfy the
triangular recurrence

$\mathscr{C}_{n+1,k}(\sigma,\phi)=\mathscr{C}_{n,k-1}(\sigma,\phi)+n^{1-\sigma}\phi^{-n}\mathscr{C}_{n,k}(\sigma,\phi),$

for any $n\geq 0$ and $1\leq k\leq n+1$, with initial conditions
$\mathscr{C}_{0,0}(\sigma,\phi)=1,$ $\mathscr{C}_{n,0}(\sigma,\phi)=0,$ $n\geq
1,$ $\mathscr{C}_{n,k}(\sigma,\phi)=0,$ $k>n.$ Moreover, for any $1\leq k\leq
n$ and $n\geq 2$, one has

$\mathscr{C}_{n,k}(\sigma,\phi)=\sum_{(i_{1},\ldots,i_{n-k})}\prod_{j=1}^{n-k}i_{j}^{1-\sigma}\phi^{-i_{j}},$

where the sum runs over the $(n-k)$-combinations of integers
$(i_{1},\dots,i_{n-k})$ in $\\{1,\dots,n-1\\}$.

We can now state the main theoretical result, namely the probability mass
function of $K_{n}$, which can be expressed in terms of the coefficients
$\mathscr{C}_{n,k}(\sigma,\phi)$.

######  Theorem 4.

Let
$K_{n}\sim\textsc{pb}\\{1,S(1;\alpha,\sigma,\phi),\ldots,S(n-1;\alpha,\sigma,\phi)\\}$
for every $n\geq 1$. Then,

$\mathrm{pr}(K_{n}=k)=\frac{\alpha^{k}}{\prod_{i=0}^{n-1}(\alpha+i^{1-\sigma}\phi^{-i})}\mathscr{C}_{n,k}(\sigma,\phi).$

Theorem 4 reduces to the distribution obtained by Antoniak (1974) and recalled
in equation (3) when $\sigma=0$ and $\phi=1$. Gibbs-type priors enjoy a
similar structure for the distribution of $K_{n}$, having replaced
$\mathscr{C}_{n,k}(\sigma,\phi)$ with the so-called generalized factorial
coefficients, which have similar properties; see Gnedin and Pitman (2005); De
Blasi et al. (2015) for further discussion.

###  3.3 Covariate-dependent models

Under the three parameter log-logistic specification, the discovery
probabilities are
$\pi_{n+1}=\mathrm{pr}(D_{n+1}=1)=\alpha\phi^{n}(\alpha\phi^{n}+n^{1-\sigma})^{-1}$
for $n\geq 1$ with $\pi_{1}=1$. An interesting and practically useful property
of our model is the following representation

$\log\frac{\pi_{n+1}}{1-\pi_{n+1}}=\log\alpha-(1-\sigma)\log
n+(\log\phi)n=\beta_{0}+\beta_{1}\log{n}+\beta_{2}n,\qquad n\geq 1,$ (11)

having set $\beta_{0}=\log\alpha$, $\beta_{1}=\sigma-1<0$ and
$\beta_{3}=\log\phi\leq 0$. Hence, equation (11) has the form of a logistic
regression for the binary indicators $D_{2},\ldots,D_{n}$, where the
regression coefficients $\beta_{2}$ and $\beta_{3}$ are constrained to be
negative. By letting $\beta_{1}=-1$ and $\beta_{2}=0$ one recovers the
discovery probability of the Dirichlet process. This representation has
computational advantages, which will be discussed in Section 4.

The logistic regression representation in (11) suggests natural extensions to
accommodate covariates. Suppose we are given a collection of $L$ accumulation
curves, namely $(K_{1n})_{n\geq 1},\dots,(K_{Ln})_{n\geq 1}$, representing for
example the sequential discoveries recorded at different geographical
locations. Each location is associated with a set of covariates
$z_{\ell}^{\mathrm{T}}=(z_{\ell 1},\ldots,z_{\ell p})\in\mathds{R}^{p}$ for
$\ell=1,\dots,L$. Let $(D_{\ell n})_{n\geq 1}$ be the sequence of discovery
indicators for the $\ell$th location, with probabilities $(\pi_{\ell
n})_{n\geq 1}$. The most flexible specification for $K_{\ell n}$ corresponds
to the case in which all the parameters are location-specific, so that for any
$n\geq 1$,

$\log\frac{\pi_{\ell n+1}}{1-\pi_{\ell n+1}}=\beta_{\ell 0}+\beta_{\ell
1}\log{n}+\beta_{\ell 2}n,\qquad(\ell=1,\dots,L).$

This specification can borrow information across locations via a hierarchical
model on $\beta_{l}=(\beta_{l0},\beta_{l1},\beta_{l2})^{\mathrm{T}}$ or by
fixing certain parameters. Alternatively, systematic variation across
locations can be modeled through including covariates $z_{\ell}$ via

$\log\frac{\pi_{\ell n+1}}{1-\pi_{\ell n+1}}=\beta_{\ell 0}+\beta_{\ell
1}\log{n}+\beta_{\ell
2}n=z_{\ell}^{\mathrm{T}}\gamma_{0}+(z_{\ell}^{\mathrm{T}}\gamma_{1})\log{n}+(z_{\ell}^{\mathrm{T}}\gamma_{2})n,$
(12)

for $\ell=1,\dots,L$, with $\gamma_{0},\gamma_{1},\gamma_{2}\in\mathds{R}^{p}$
being vectors of coefficients such that $z_{\ell}^{\mathrm{T}}\gamma_{2}<0$
and $z_{\ell}^{\mathrm{T}}\gamma_{2}\leq 0$. This specification is still in
the form of a logistic regression and therefore inference on the parameters
$\gamma_{0},\gamma_{1}$ and $\gamma_{2}$ can be conducted through
straightforward modifications of standard algorithms. The computational
details are discussed in the next Section.

##  4 Posterior computation

###  4.1 Estimation procedures

Consider the model in equation (11). The parameters
$\theta=(\alpha,\sigma,\phi)$ can be estimated by maximizing the likelihood in
equation (7), with $S(t;\theta)=S(t;\alpha,\sigma,\phi)$, $\beta_{1}<0$ and
$\beta_{2}\leq 0$. In practice, it may suffice to ignore these constraints and
apply routine algorithms for fitting logistic regression, as the unconstrained
maximum likelihood estimates typically satisfy the constraints. In this case,
the resulting estimate $\hat{\theta}$ has the following appealing property.

######  Proposition 4.

Let $\hat{\theta}=(\hat{\alpha},\hat{\sigma},\hat{\phi})$ be the unconstrained
maximizer of equation (7) under the three-parameter specification in (10), if
it exists. If $K_{n}=k$ is the number of discoveries within the data
$D_{1},\dots,D_{n}$, then the expectation $E(K_{n})$, evaluated at
$\hat{\theta}$, equals $k$.

In other words, the $n$th term of the smoothed accumulation curve $E(K_{n})$
matches the total number of distinct labels observed in the sequence when the
parameters are estimated through unconstrained maximum likelihood.

Although we can obtain confidence intervals and standard errors for the
parameters using the aforementioned maximum likelihood strategy, conducting
inferences in this manner ignores the parameter constraints. In contrast, a
fully Bayesian approach can easily incorporate them through a prior, such as
$\beta\sim N(\mu,\Sigma)\mathds{1}(\beta_{1}<0;\beta_{2}\leq 0).$ Under this
prior, a straightforward modification of the Pólya-gamma data-augmentation
strategy of Polson et al. (2013) can be used for posterior sampling. The
covariate-dependent regression detailed in equation (B.2) can be naturally
carried out in a similar manner. For additional details refer to the
Supplementary Material.

###  4.2 Simulations

We test our log-logistic models on four synthetic sequences of length
$n=90,000$. These sequences are simulated according to four different models:
i) Dirichlet process with $\alpha=30$, ii) Pitman–Yor process with $\alpha=30$
and $\sigma=0.25$, iii) Dirichlet-multinomial process with $\sigma=-0.25$ and
$\alpha=H|\sigma|$ and $H=5,000$, and iv) drawing the species from a Zipf
distribution with support $\\{1,\dots,H\\}$ with $H=5,000$ and shape parameter
$0.3$. Our log-logistic formulations display excellent performance in each
case even though the generating mechanisms are species sampling models. Cases
i)-ii) have $K_{\infty}=\infty$ while for iii)-iv) we get $K_{\infty}=5,000$.
We estimate the parameters of each model on the first third of each sequence,
comprising $30,000$ data points, and then assess predictive performance on the
remaining $60,000$ observations.

We chose truncated and independent normal priors centered at $0$ and with
standard deviation $10$. We run a Markov Chain Monte Carlo algorithm for a
total of $15,000$ iterations, discarding the first $5,000$ samples. We
estimate out-of-sample accumulation curves by averaging over the posterior
samples of $E(K_{n+m}\mid D_{1},\dots,D_{n})$, obtained as in Proposition 1.
To compare different log-logistic specifications, we use the Deviance
Information Criterion (dic) of Spiegelhalter et al. (2002). In Table 1 the dic
values and the absolute deviations between the predicted and the true values
of $K_{n+m}$ are reported, with $n=30,000$ and $m\in\\{n/3,n,2n\\}$.

Table 1: Model performances over simulated sequences of length 90,000. Estimates are based on the first 30,000 observations. |  | Posterior means |  | Average prediction error
---|---|---|---|---
Data | Model | $\alpha$ | $\sigma$ | $\phi$ | dic | $m=n/3$ | $m=n$ | $n=2n$
Dirichlet |  ll-1 | $33.47$ | - | - | $1,983.60$ | $2.38$ | $0.82$ | $1.75$
|  ll-2 | $31.07$ | $0.02$ | - | $1,985.34$ | $1.84$ | $0.57$ | $4.11$
|  ll-3 | $25.22$ | $0.06$ | $0.99$ | $1,986.59$ | $4.15$ | $6.58$ | $9.93$
Pitman-Yor |  ll-1 | $111.49$ | - | - | $5,187.99$ | $24.03$ | $51.93$ | $96.79$
|  ll-2 | $24.08$ | $0.20$ | - | $5,130.76$ | $3.37$ | $3.2$ | $0.49$
|  ll-3 | $20.67$ | $0.22$ | $0.99$ | $5,132.38$ | $9.13$ | $16.74$ | $44.49$
Dir-multinomial |  ll-1 | $760.92$ | - | - | $14,392.07$ | $63.18$ | $193.96$ | $327.3$
|  ll-2 | $2,239.38$ | $-0.12$ | - | $14,351.46$ | $19.51$ | $78.05$ | $127.24$
|  ll-3 | $1,512.80$ | $-0.07$ | $0.99$ | $14,350.51$ | $4.26$ | $0.51$ | $34.81$
Zipf |  ll-1 | $1,699.21$ | - | - | $18,669.83$ | $446.91$ | $1,104.64$ | $1,777.94$
|  ll-2 | $9.2\times 10^{5}$ | $-0.71$ | - | $17,495.11$ | $135.63$ | $299.57$ | $429.15$
|  ll-3 | $5,087.35$ | $-0.01$ | $0.99$ | $17,274.05$ | $17.63$ | $22.80$ | $24.10$

Throughout the rest of the paper, we refer to the Dirichlet process and two-
and three-parameter log-logistic models as ll-1, ll-2 and ll-3, respectively.
When data are generated according to a Dirichlet process, ll-1 has the lowest
dic and prediction errors. Prediction for ll-2 closely resembles that for
ll-1, as the posterior mean for $\sigma$ is close to $0$. When data are
generated according to a Pitman–Yor, ll-2 achieves the highest accuracy, as
expected, since our model also accounts for polynomial growth rates for
$K_{n}$.

When the true number of species is finite, our ll-3 model has much better
performance compared to ll-1 and ll-2. For the sequence from the Dirichlet-
multinomial process, there is little difference between ll-2 and ll-3 in terms
of dic. However, the increased flexibility of ll-3 leads to a much higher out-
of-sample accuracy. This behavior is even more evident in the sequence
generated from a Zipf distribution. In this case, all $5,000$ species are
observed in the first $30,000$ samples so that the true accumulation curve
$K_{n+m}$ is a horizontal line.

##  5 Fungal biodiversity application

We analyze data from a fungi biodiversity study in Finland (Abrego et al.,
2020). Each sample contains a large number of fungal dna barcode sequences
obtained either from air samples or soil samples. The total number of
sequences to be obtained by each sample can be controlled for when performing
the high-throughput sequencing, with the cost of sequencing increasing with
the number of sequences. To optimize the number of sequences per sample, one
would like to know how many species are to be missed for a particular
sequencing depth. This is the goal of our analysis.

The data consist of 174 different samples from different sites across five
cities in Finland. For each site, fungi samples are collected on the same
dates at two urban areas, one at the core and one at the edge of the city, and
two nearby natural areas, again with one at the core and one at the edge. Two
different sampling methods were used: i) through air, via a cyclone trap and
continuously for 24 hours, and ii) through soil, gathering a small portion of
soil close to the air trap. We exclude samples with less than $1,000$
sequences, as in such cases the samples lacked sufficient numbers of spores
for more comprehensive barcoding.

This leaves us with a total of 166 samples; the average number of barcoded DNA
sequences per sample is $111,980$ and the average number of species discovered
is $1,184$. We fit the three different log-logistic models to training data
containing the first third of the fungal barcode sequences, allowing each
sample to have its own parameters. Model fitting and prediction proceeded
exactly as in Section 4.2.

In Table 2 we report the percentage absolute errors between the predicted and
true values $K_{n}$ and $K_{n+m}$, averaged over the 166 samples. As in Figure
1, model ll-1 obtained poor in-sample fit. Model ll-2 has better in-sample
accuracy but exhibits explosive out-of-sample behavior. Model ll-3 is highly
accurate both in- and out-of-sample, while having the lowest dic. This is
confirmed by Figure 2, which displays the performance of the ll-3 model on the
same data as used in Figure 1. In Table 3 we report in-sample performance of
the three models using all the available DNA sequences. Again, ll-3 has the
lowest dic. The ll-1 model strongly over-predicts the initial part of the
curve as before. Although ll-2 represents a slight improvement over ll-1, the
three-parameter log-logistic ll-3 has uniformly better performance.

Table 2: Summary of the model performances across the 166 samples. Estimates use the first $1/3$ of the DNA sequences in each sample. |  | Fraction of the curve
---|---|---
Model | dic | $0.1$ | $0.25$ | $0.33$ | $0.50$ | $0.66$ | $0.75$ | $1.00$
|  | Average percentage errors
ll-1 | $13.61\times 10^{5}$ | $17.06$ | $3.54$ | $0.10$ | $3.31$ | $3.86$ | $3.76$ | $2.83$
ll-2 | $13.19\times 10^{5}$ | $3.37$ | $1.91$ | $0.10$ | $4.03$ | $8.58$ | $11.25$ | $20.25$
ll-3 | $13.17\times 10^{5}$ | $1.55$ | $0.98$ | $0.10$ | $1.94$ | $2.97$ | $3.5$ | $5.27$
Table 3: Summary of the model performances across the 166 samples. Estimates are based on all the available DNA sequences. |  | Fraction of the curve
---|---|---
Model | dic | $0.1$ | $0.25$ | $0.33$ | $0.50$ | $0.66$ | $0.75$ | $1.00$
|  | Average percentage errors
ll-1 | $18.29\times 10^{5}$ | $18.07$ | $4.94$ | $2.68$ | $3.01$ | $3.31$ | $2.91$ | $0.12$
ll-2 | $18.03\times 10^{5}$ | $3.34$ | $7.61$ | $8.6$ | $8.35$ | $6.5$ | $5.20$ | $0.16$
ll-3 | $17.78\times 10^{5}$ | $2.47$ | $2.13$ | $1.38$ | $0.67$ | $1.00$ | $1.06$ | $0.12$
Figure 2: Performance of ll-3 on the same data as Figure 1. The dots represent
the observed values. Left panel: the solid line is the predicted in-sample and
out-of-sample accumulation curve computed by averaging over posterior samples
of $E(K_{n})$ and $E(K_{n+m}\mid D_{1},\dots,D_{n})$, respectively. Black
dashed lines indicate the 95% posterior predictive credible intervals. Right
panel: the solid line is the predicted accumulation curve obtained averaging
posterior samples of $E(K_{n})$ Figure 3: Left panel: distribution of the
posterior mean sample species richness for the 166 samples. Right panel:
distribution of the posterior mean sample saturation for the 166 samples.

We rely on model ll-3 in performing inferences on i) the sample species
richness, which is the total number of species that can be detected through
barcoding within a sample, and ii) whether dna barcoding has reached
saturation at different sites, meaning that only very few species are missed.
To address i), we estimate the posterior mean $E(K_{\infty}\mid
D_{1},\ldots,D_{n})$ for each individual sample, which is guaranteed to be
finite. The results are reported in the left panel of Figure 3, which displays
the expected sample species richness for each of the 166 samples across site
characteristics. Air samples tend to contain more species, and there is some
evidence of greater species richness in natural environments, as reported by
Abrego et al. (2020).

For task ii), let $C_{n}=K_{n}/K_{\infty}\leq 1$ represent the saturation
level of a given sample after $n$ barcoded sequences. Differences across sites
can be evaluated via $E(C_{n}\mid D_{1},\ldots,D_{n})$, which represents the
posterior expected saturation level of a sample. Figure 3, right panel,
summarizes posterior mean saturation stratified by sampling site
characteristics. While there is some variability across sites, most of them
have a ratio greater than $0.95$. The results suggest that if additional dna
sequences are barcoded there is the opportunity to detect approximately
$2.5$-$5$% more species in each sample. We can estimate the number $m$ of
additional sequences that would need to be barcoded to reach a desired
saturation level $C_{n+m}$. For example, for the sample highlighted in Figures
1 and 2, the saturation level is $98.0$%. To achieve a saturation level of
$99.5$% we would need to barcode an additional $41,547$ sequences, an increase
of $43\%$.

## Acknowledgements

This project has received funding from the European Research Council under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 856506).

##  Appendix A Proofs

###### Proof of Theorem 1.

We first discuss the likelihood $\mathscr{L}(\alpha\mid X_{1},\ldots,X_{n})$.
Given the sequence of tags $X_{1},\ldots,X_{n}$, any exchangeable prediction
scheme defines a random partition $\Psi_{n}$ of the integers $\\{1,\dots,n\\}$
such that $i$ and $j$ belong to the same set in $\Psi_{n}$ if and only if
$X_{i}=X_{j}$. Let $\\{C_{1},\dots,C_{k}\\}$ be a partition of $n$ into $k$
groups, with $n_{j}=\text{card}(C_{j})$, $(j=1,\dots,k)$ being the cardinality
of $C_{j}$. The resulting law of the random partition $\Psi_{n}$ in the
Dirichlet process case equals
$\textrm{pr}(\Psi_{n}=\\{C_{1},\dots,C_{k}\\})=\alpha^{k}/(\alpha)_{n}\prod_{j=1}^{k}(n_{j}-1)!$.
This is the likelihood function for $\alpha$, so that $\mathscr{L}(\alpha\mid
X_{1},\dots,X_{n})\propto\alpha^{k}/(\alpha)_{n}$, which only depends on
$K_{n}=k$ and not on $n_{1},\dots,n_{k}$. By letting
$(\alpha)_{n}=\alpha\prod_{i=2}^{n}(\alpha+i-1)$, we get

$\log\mathscr{L}(\alpha\mid
X_{1},\ldots,X_{n})=(k-1)\log\alpha-\sum_{i=2}^{n}\log(\alpha+i-1)+c_{X},$
(13)

where $c_{X}$ is a constant not depending on $\alpha$. On the other hand, the
logarithm of the likelihood induced by the discovery indicators is equal to

$\log\mathscr{L}(\alpha\mid
D_{1},\ldots,D_{n})=\log\alpha\sum_{i=2}^{n}D_{i}-\sum_{i=2}^{n}\log(\alpha+i-1)+c_{D},$
(14)

with $c_{D}$ being a constant not depending on $\alpha$. Since
$\sum_{i=2}^{n}D_{i}=k-1$, one has that equation (13) and (14) are equal up to
an additive constant. Thus, the result follows. ∎

###### of Proposition 1.

This proof follows from the definition of a Poisson-binomial distribution. In
particular, for every $n,m\geq 1$, we have that
$K_{m}^{(n)}=K_{n+m}-K_{n}=\sum_{i=n+1}^{n+m}D_{i}$ is a sum of independent
indicators with discovery probabilities $\pi_{n+1}\ldots\pi_{n+m}$ directed by
$S(t;\theta)$. Hence, from Definition 1, $(K_{m}^{(n)}\mid
D_{1},\ldots,D_{n})\sim\textsc{pb}\\{S(n;\theta),\ldots,S(n+m-1;\theta)\\}$,
whose expected value is $E(K_{m}^{(n)}\mid
D_{1},\ldots,D_{n})=\sum_{j=1}^{m}S(j+n-1;\theta)$. Finally, if $K_{n}=k$, it
naturally follows that $E(K_{n+m}\mid D_{1},\ldots,D_{n})=E(K_{n}\mid
D_{1},\ldots,D_{n})+E(K_{m}^{(n)}\mid D_{1},\ldots,D_{n})=k+E(K_{m}^{(n)}\mid
D_{1},\ldots,D_{n}).$ ∎

###### of Proposition 2.

Recall that $K_{n}=\sum_{i=1}^{n}D_{i}$ is non-decreasing in $n$. Taking the
limit as $n\to\infty$, we have that $K_{n}\to
K_{\infty}=\sum_{i=1}^{\infty}D_{i}$ almost surely. Then,
$E(K_{\infty})=\sum_{i=1}^{\infty}S(i-1;\theta)$, as a consequence of the
monotone convergence theorem. Moreover, equation (8) follows from a simple
calculus inequality: as $S(t;\theta)$ is positive and strictly decreasing in
$t$, we have that

$\int_{0}^{n-1}S(t;\theta)\mathrm{d}t+S(n-1;\theta)\leq\sum_{i=1}^{n}S(i-1;\theta)\leq\int_{0}^{n-1}S(t;\theta)\mathrm{d}t+S(0;\theta),$

for every $n\geq 1$. Taking the limit for $n\to\infty$, one has that
$S(n-1;\theta)\to 0$. But then, as $S(0;\theta)=1$, we have that $E(T)\leq
E(K_{\infty})\leq E(T)+1$, with $E(K_{\infty})$ defined above and
$E(T)=\int_{0}^{\infty}S(t;\theta)\mathrm{d}t$. ∎

###### Proof of Corollary 1.

We begin by proving that $K_{\infty}=\infty$ if and only if
$E(K_{\infty})=\infty$. One side follows from the monotone convergence
theorem: if $K_{\infty}=\infty$, then necessarily
$\lim_{n\to\infty}E(K_{n})=E(K_{\infty})=\infty$ by the same argument in the
proof of Proposition 2. The other direction can be proved by contrapposition:
suppose that $K_{\infty}<\infty$. Then, there exists a positive constant
$M<\infty$ such that $K_{\infty}<M$ almost surely. This means that
$E(K_{\infty})<E(M)=M<\infty$. The rest of the claim naturally follows from
the inequality in equation (8) of Proposition 2. ∎

###### of Corollary 2.

To prove this claim, we rely on the limit comparison test for the ratio of two
series. In particular,

$\lim_{n\to\infty}\frac{S(n-1;\theta)\\{1-S(n-1;\theta)\\}}{S(n-1;\theta)}=1-\lim_{n\to\infty}S(n-1;\theta)=1.$

This implies that
$\textup{var}(K_{\infty})=\sum_{i=1}^{\infty}S(i-1;\theta)\\{1-S(i-1;\theta)\\}$
diverges if and only if $E(K_{\infty})=\sum_{i=1}^{\infty}S(i-1;\theta)$
diverges. Following the same argument in the proof of Corollary 1, having
$K_{\infty}<\infty$ almost surely implies that $E(K_{\infty})<\infty$, and in
turn $\textup{var}(K_{\infty})<\infty$. ∎

###### Proof of Theorem 2.

The first part of the theorem is a consequence of the strong law of large
numbers for the sum of independent random variables. In particular, let
$b_{n}=\int_{1}^{n}S(t-1;\theta)\mathrm{d}t$. Then, $b_{n}<b_{n+1}$ for every
$n$, and $b_{n}\to E(T)=\infty$ as $n\to\infty$. Since by assumption
$S(n-1;\theta)>S(n;\theta)$ and $b_{n}^{2}<b_{n+1}^{2}$ for every $n$, we have
that $\sum_{n=1}^{\infty}\textup{var}(D_{n})/b_{n}^{2}<\infty$, which holds by
the series convergence test, because

$\displaystyle\lim_{n\to\infty}\frac{\textup{var}(D_{n+1})}{b_{n+1}^{2}}\frac{b_{n}^{2}}{\textup{var}(D_{n})}$
$\displaystyle=\lim_{n\to\infty}\frac{S(n;\theta)\\{1-S(n;\theta)\\}}{S(n-1;\theta)\\{1-S(n-1;\theta)\\}}\frac{b_{n}^{2}}{b_{n+1}^{2}}<\lim_{n\to\infty}\frac{1-S(n;\theta)}{1-S(n-1;\theta)}=1.$

Hence, the above condition ensures that $\\{K_{n}-E(K_{n})\\}/b_{n}\to 0$
almost surely as $n\to\infty$ by the strong law of large numbers. This means
that $\lim_{n\to\infty}K_{n}/b_{n}=\lim_{n\to\infty}E(K_{n})/b_{n}=1$, almost
surely, as a consequence of Proposition 2.

The second part of the claim follows from Lyapunov’s central limit theorem.
Define $\sigma^{2}_{n}=\textup{var}(K_{n})$ for every $n$. As the discovery
indicators $(D_{n})_{n\geq 1}$ are all independent, we can prove the central
limit theorem for $K_{n}$ by showing that there exists a $\delta>0$ such that
$\lim_{n\to\infty}1/\sigma_{n}^{2+\delta}\sum_{i=1}^{n}E(|D_{i}-\pi_{i}|^{2+\delta})=0$,
where $\pi_{n}=S(n-1;\theta)$ is the discovery probability at every $n$. Fix
$\delta=2$. From the proofs of Corollaries 1 and 2, we have that
$K_{\infty}=\infty$ implies that $\lim_{n\to\infty}\sigma^{2}_{n}=\infty$.
Moreover, by looking at the fourth centered moment of a Bernoulli distribution
we have that

$\sum_{i=1}^{n}E(|D_{i}-\pi_{i}|^{4})=\sum_{i=1}^{n}\pi_{i}(1-\pi_{i})\\{1-3\pi_{i}(1-\pi_{i})\\}\leq\sum_{i=1}^{n}\pi_{i}(1-\pi_{i})=\sigma^{2}_{n},$

which leads to
$0\leq\lim_{n\to\infty}1/\sigma_{n}^{4}\sum_{i=1}^{n}E(|D_{i}-\pi_{i}|^{4})\leq\lim_{n\to\infty}1/\sigma_{n}^{2}=0$,
concluding the proof. ∎

###### Proof of Proposition 3.

This can be proved by means of the series convergence test. By the fact that

$\lim_{n\to\infty}\frac{S(n;\alpha,\sigma,\phi)}{S(n-1;\alpha,\sigma,\phi)}=\lim_{n\to\infty}\frac{\alpha\phi^{n}}{\alpha\phi^{n}+n^{1-\sigma}}\frac{\alpha\phi^{n-1}+(n-1)^{1-\sigma}}{\alpha\phi^{n-1}}=\phi,$

having $\phi<1$ implies that
$E(K_{\infty})=\sum_{i=1}^{\infty}S(i-1;\alpha,\sigma,\phi)<\infty$ almost
surely. But then, $K_{\infty}<\infty$ as well by the proof of Corollary 1. ∎

###### Proof of Theorems 3 and 4.

The proofs of Theorems 3 and 4 are presented together. The arguments we use
follow a similar line of reasoning as in Charalambides (2005). As a first
step, we prove the triangular recurrence in Theorem 3. Following Definition 2,
we can write
$\prod_{k=0}^{n}(\alpha+k^{1-\sigma}\phi^{-k})=(\alpha+n^{1-\sigma}\phi^{-n})\prod_{k=0}^{n-1}(\alpha+k^{1-\sigma}\phi^{-k}),$
for any $n\geq 1$, from which it follows that

$\displaystyle\sum_{k=0}^{n+1}\alpha^{k}\mathscr{C}_{n+1,k}(\sigma,\phi)=\sum_{k=1}^{n+1}\alpha^{k}\mathscr{C}_{n,k-1}(\sigma,\phi)+\sum_{k=0}^{n}\alpha^{k}n^{1-\sigma}\phi^{-n}\mathscr{C}_{n,k}(\sigma,\phi).$

Hence, all the coefficients associated to each $\alpha^{k}$ must coincide
under both sides of the above equation. This means that
$\mathscr{C}_{n+1,k}(\sigma,\phi)=\mathscr{C}_{n,k-1}(\sigma,\phi)+n^{1-\sigma}\phi^{-n}\mathscr{C}_{n,k}(\sigma,\phi)$.
As for the initial conditions, it is easy to check that they naturally follow
from Definition 2.

To prove the second part of Theorem 3, we start by considering
$\mathrm{pr}(K_{n}=k)$. Call $j_{1},\ldots,j_{n}$ a sequence of indexes such
that $D_{j_{s}}=1$ for $s=1,\ldots,k$, and $D_{j_{s}}=0$ for $s=k+1,\ldots,n$.
By independence of the indicators, the probability of such a configuration is

$\displaystyle\mathrm{pr}(D_{j_{1}}=0,\ldots,D_{j_{k}}=1,D_{j_{k+1}}=0,\ldots,D_{j_{n}}=0)=$
$\displaystyle\quad=\prod_{s=1}^{k}S(j_{s}-1;\alpha,\sigma,\phi)\prod_{s=k+1}^{n}\\{1-S(j_{s}-1;\alpha,\sigma,\phi)\\}=\frac{\alpha^{k}}{\prod_{i=0}^{n-1}(\alpha+i^{1-\sigma}\phi^{-i})}\prod_{j=1}^{n-k}i_{j}^{1-\sigma}\phi^{-i_{j}},$

where the product in the last equality follows from relabeling the indexes as
$i_{1}=j_{k+1}-1,\ldots,i_{n-k}=j_{n}-1$. Moreover, note that
$\\{i_{1},\ldots,i_{n-k}\\}$ is one of the $n-k$ possible combinations of the
$n-1$ positive integers $\\{1,\ldots,n-1\\}$ for which we obtain precisely $k$
discoveries, with $1\leq k\leq n$ and $n\geq 2$. Hence, summing over all the
possible combinations of $\\{i_{1},\ldots,i_{n-k}\\}$ leads us to the
probability

$\mathrm{pr}(K_{n}=k)=\frac{\alpha^{k}}{\prod_{i=0}^{n-1}(\alpha+i^{1-\sigma}\phi^{-i})}\sum_{(i_{1},\ldots,i_{n-k})}\prod_{j=1}^{n-k}i_{j}^{1-\sigma}\phi^{-i_{j}},$
(15)

for $1\leq k\leq n$ and $n\geq 2$. The object in equation (15) is a
probability mass function. This means that

$\sum_{k=0}^{n}\mathrm{pr}(K_{n}=k)=\sum_{k=0}^{n}\frac{\alpha^{k}}{\prod_{i=0}^{n-1}(\alpha+i^{1-\sigma}\phi^{-i})}\sum_{(i_{1},\ldots,i_{n-k})}\prod_{j=1}^{n-k}i_{j}^{1-\sigma}\phi^{-i_{j}}=1,$

recalling that $\mathrm{pr}(K_{n}=0)=0$. Rearranging the equality, one has
that

$\prod_{i=0}^{n-1}(\alpha+i^{1-\sigma}\phi^{-i})=\sum_{k=0}^{n}\alpha^{k}\sum_{(i_{1},\ldots,i_{n-k})}\prod_{j=1}^{n-k}i_{j}^{1-\sigma}\phi^{-i_{j}},$

which is the same polynomial expansion proposed in Definition 2. Hence, it
must be that

$\mathscr{C}_{n,k}(\sigma,\phi)=\sum_{(i_{1},\ldots,i_{n-k})}\prod_{j=1}^{n-k}i_{j}^{1-\sigma}\phi^{-i_{j}},$
(16)

again for $1\leq k\leq n$ and $n\geq 2$. This last equality proves the second
part of Theorem 3. Finally, Theorem 4 naturally follows by plugging equation
(16) into (15). ∎

###### Proof of Proposition 4.

To find the maximizer $\hat{\theta}=(\hat{\alpha},\hat{\sigma},\hat{\phi})$ of
the likelihood in equation (7) with the the three-parameter log-logistic
specification we rely on the first order condition with respect to $\alpha$.
In particular, the logarithm of the likelihood becomes

$\log\mathscr{L}(\alpha,\sigma,\phi\mid D_{1},\ldots
D_{n})=\log\alpha\sum_{i=2}^{n}D_{i}-\sum_{i=2}^{n}\log\\{\alpha\phi^{i-1}+(i-1)^{1-\sigma}\\}+c_{\sigma,\phi},$

where $c_{\sigma,\phi}$ is a constant not dependent on $\alpha$. Hence, the
first order condition with respect to $\alpha$ leads to

$\sum_{i=1}^{n}\frac{\alpha\phi^{i-1}}{\alpha\phi^{i-1}+(i-1)^{1-\sigma}}=\sum_{i=1}^{n}D_{i}=k=E(K_{n}).$

This equality must be maintained at the solution $\hat{\theta}$.

∎

##  Appendix B Supplementary material

###  B.1 Posterior inference for a single accumulation curve

In the following we describe the estimation of the parameters
$\theta=(\alpha,\sigma,\phi)$ under the three-parameter log-logistic
specification and using Markov Chain Monte Carlo. Let $(D_{n})_{n\geq 1}$ be a
sequence of discovery indicators with $D_{1}=1$ and

$\pi_{n+1}=\mathrm{pr}(D_{n+1}=1\mid
D_{1},\ldots,D_{n})=\frac{\alpha\phi^{n}}{\alpha\phi^{n}+n^{1-\sigma}},\qquad
n\geq 1,$

for $\alpha>0$, $\sigma<1$, $0<\phi\leq 1$ and $\pi_{1}=1$. As discussed in
the manuscript, this implies that

$\log\frac{\pi_{n+1}}{1-\pi_{n+1}}=\log\alpha-(1-\sigma)\log
n+(\log\phi)n=\beta_{0}+\beta_{1}\log{n}+\beta_{2}n,\qquad n\geq 1,$

with $\beta_{0}=\log\alpha$, $\beta_{1}=\sigma-1<0$ and
$\beta_{3}=\log\phi\leq 0$. These constraints are imposed through a truncated
normal prior, namely $\beta\sim
N(\mu,\Sigma)\mathds{1}(\beta_{1}<0;\beta_{2}\leq 0)$.

Samples from the posterior can be easily obtained via the Pólya-gamma data-
augmentation strategy introduced in Polson et al. (2013). This procedure
introduces Pólya-gamma distributed positive latent variables
$\omega=(\omega_{2},\ldots,\omega_{n})^{\mathrm{T}}$. The resulting full
conditional distributions for $\beta$ and $\omega$ are available in closed
form. Let $d=(d_{2},\ldots,d_{n})^{\mathrm{T}}$ be the observed values for the
discovery indicators $D_{2},\ldots,D_{n}$ and let $V$ be the design matrix,
with $n-1$ rows and $3$ columns and entries
$v_{i}=(1,\log{i},i)^{\mathrm{T}}$, for $i=1,\ldots,n-1$. Then, the full
conditional for $(\beta\mid\omega,d)$ is a multivariate truncated normal
distribution with parameters $\mu_{\omega}$ and $\Sigma_{\omega}$ equal to

$\Sigma_{\omega}=(V^{\mathrm{T}}\Omega
V+\Sigma^{-1}\mu)^{-1},\quad\mu_{\omega}=\Sigma_{\omega}(V^{\mathrm{T}}\kappa+\Sigma^{-1}\mu),$
(17)

with $\kappa=(d_{2}-1/2,\ldots,d_{n}-1/2)^{\mathrm{T}}$ and
$\Omega=\mathrm{diag}(\omega_{2},\ldots,\omega_{n})$. The algorithm below
outlines the sampling procedure.

Algorithm 1 Pólya-Gamma Gibbs sampler for single site accumulation curve

1:Set an initial value $\beta$ and set number of samples $R$

2:for $r=1$ to $r=R$ do

3: for $i=1$ to $i=n-1$ do

4: Sample
$(\omega_{i}\mid\beta)\sim\textrm{PolyaGamma}(1,v_{i}^{\mathrm{T}}\beta)$

5: Sample $(\beta\mid\omega,d)\sim
N(\mu_{\omega},\Sigma_{\omega})\mathds{1}({\beta_{1}<0,\beta_{2}\leq 0})$,
with $\mu_{\omega}$, $\Sigma_{\omega}$ as in (17)

6:Output: collection of $R$ samples for $\beta$

In our work we obtain samples from the multivariate truncated normal through
the efficient algorithm proposed in Botev (2017).

###  B.2 Posterior sampling for multi-site data

We now describe a Markov Chain Monte Carlo algorithm for Bayesian inference
for the covariate-dependent model described in the manuscript. Recall that we
are given a collection of $L$ accumulation curves $(K_{1n})_{n\geq
1},\dots,(K_{Ln})_{n\geq 1}$ observed up to the terms $n_{1},\dots,n_{L}$.
Each curve is associated to a set of covariates $z_{\ell}=(z_{\ell
1},\ldots,z_{\ell p})^{\mathrm{T}}$ for $\ell=1,\dots,L$. Future observations
correspond to new discoveries within the set of the considered $L$ curves, so
that new covariates values are not expected.

Let $(D_{\ell n})_{n\geq 1}$ be the sequence of discovery indicators for the
$\ell$th location, with probabilities $(\pi_{\ell n})_{n\geq 1}$. Hence, we
get

$\log\frac{\pi_{\ell n+1}}{1-\pi_{\ell n+1}}=\beta_{\ell 0}+\beta_{\ell
1}\log{n}+\beta_{\ell
2}n=z_{\ell}^{\mathrm{T}}\gamma_{0}+(z_{\ell}^{\mathrm{T}}\gamma_{1})\log{n}+(z_{\ell}^{\mathrm{T}}\gamma_{2})n,$

with $\gamma_{0},\gamma_{1},\gamma_{2}\in\mathds{R}^{p}$ coefficient vectors
such that $z_{\ell}^{\mathrm{T}}\gamma_{2}<0$ and
$z_{\ell}^{\mathrm{T}}\gamma_{2}\leq 0$ for every $\ell=1,\dots,L$. The above
specification is a logistic regression and therefore inference on the
parameters $\gamma=(\gamma_{0},\gamma_{1},\gamma_{2})^{\mathrm{T}}$ may be
conducted through a simple modification of Algorithm 1.

Let $N=\sum_{\ell=1}^{L}(n_{\ell}-1)$ and let $V$ be a design matrix with $N$
rows and $3p$ columns, with rows $v_{(\ell
i)}=(z^{\mathrm{T}}_{\ell},z^{\mathrm{T}}_{\ell}\log{i},z^{\mathrm{T}}_{\ell}i)^{\mathrm{T}}$
for $i=1,\ldots,n_{\ell}-1$ and $\ell=1,\dots,L$. Moreover, call
$d=(d_{1}^{\mathrm{T}},\ldots,d_{L}^{\mathrm{T}})^{\mathrm{T}}$ the realized
discoveries, with $d_{\ell}=(d_{\ell 2},\ldots,d_{\ell
n_{\ell}})^{\mathrm{T}}$ be the observed values for $D_{\ell 1},\ldots,D_{\ell
n}$, for every $\ell=1,\ldots,L$. As before, we can incorporate the
constraints $z^{\mathrm{T}}_{\ell}\gamma_{1}<0$ and
$z^{\mathrm{T}}_{\ell}\gamma_{2}\leq 0$ by assigning $\gamma$ a multivariate
truncated normal prior,

$\gamma\sim
N(\mu,\Sigma)\mathds{1}(z^{\mathrm{T}}_{\ell}\gamma_{1}<0;z^{\mathrm{T}}_{\ell}\gamma_{2}\leq
0;\ell=1,\ldots,L).$

Let $\omega$ be a $N$-dimensional vector of Pólya-gamma latent variables.
Then, the full conditional for $(\gamma\mid\omega,d)$ is a multivariate
truncated normal distribution with mean $\mu_{\omega}$ and covariance matrix
$\Sigma_{\omega}$ equal to equation (17), whereas $\Omega$ is a diagonal
matrix whose diagonal elements are those of the vector $\omega$.

In Algorithm 2 we employ a vanilla acceptance rejection sampler for the full
conditional $(\gamma\mid\omega,d)$. This is indeed a reasonable approach in
most practical settings, as the data usually support the required constraints,
leading to very high acceptance rates. If needed, suitable adaptations of the
ideas of Botev (2017) may be alternatively considered.

Algorithm 2 Pólya-Gamma Gibbs sampler for covariate-dependent accumulation
curves

1:Set an initial value $\gamma$ and set the number of samples $R$

2:for $r=1$ to $r=R$ do

3: for $\ell=1$ to $\ell=L$ do

4: Sample $(\omega_{(\ell i)}\mid\gamma)\sim\textrm{PolyaGamma}(1,v_{(\ell
i)}^{\mathrm{T}}\gamma),\quad i=1,\ldots,n_{\ell}-1$

5: Sample $(\gamma\mid\omega,d)\sim N(\mu_{\omega},\Sigma_{\omega})$, with
$\mu_{\omega}$, $\Sigma_{\omega}$ as in (17) until $\gamma$

6: satisfies $z^{\mathrm{T}}_{\ell}\gamma_{1}<0$ and
$z^{\mathrm{T}}_{\ell}\gamma_{2}\leq 0$ for every $\ell=1,\ldots,L$.

7:Output: collection of $R$ samples for $\gamma$

## References

  * Abrego et al. (2020) Abrego, N., B. Crosier, P. Somervuo, N. Ivanova, A. Abrahamyan, A. Abdi, K. Hämäläinen, K. Junninen, M. Maunula, J. Purhonen, and O. Ovaskainen (2020).  Fungal communities decline with urbanization – more in air than soil.  ISME J. 14, 2806–15.
  * Antoniak (1974) Antoniak, C. E. (1974, 11).  Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems.  Ann. Statist. 2(6), 1152–74.
  * Blackwell and MacQueen (1973) Blackwell, D. and J. B. MacQueen (1973).  Ferguson distributions via Pólya urn schemes.  Ann. Statist. 1(2), 353–55.
  * Botev (2017) Botev, Z. I. (2017).  The normal law under linear restrictions: simulation and estimation via minimax tilting.  J. R. Statist. Soc. B 79(1), 125–148.
  * Bunge and Fitzpatrick (1993) Bunge, J. and M. Fitzpatrick (1993).  Estimating the number of species: a review.  J. Am. Statist. Assoc. 88(421), 364–73.
  * Camerlenghi et al. (2020) Camerlenghi, F., B. Dumitrascu, F. Ferrari, B. E. Engelhardt, and S. Favaro (2020).  Nonparametric Bayesian multi-armed bandits for single cell experiment design.  Ann. Appl. Stat. In press.
  * Camerlenghi et al. (2018) Camerlenghi, F., A. Lijoi, and I. Prünster (2018).  Bayesian nonparametric inference beyond the Gibbs-type framework.  Scand. J. Statist. 45, 1062–91.
  * Chao and Shen (2004) Chao, A. and T.-J. Shen (2004, 09).  Nonparametric prediction in species sampling.  J. Agric. Biol. and Envir. Statist. 9, 253–69.
  * Charalambides (2005) Charalambides, C. A. (2005, 06).  Combinatorial Methods in Discrete Distributions.  Hoboken, NJ: Wiley.
  * De Blasi et al. (2015) De Blasi, P., S. Favaro, A. Lijoi, R. H. Mena, I. Prünster, and M. Ruggiero (2015).  Are Gibbs-type priors the most natural generalization of the Dirichlet process?  IEEE Trans. Pat. Anal. Mach. Intel. 37(2), 212–29.
  * Efron and Thisted (1976) Efron, B. and R. Thisted (1976).  Estimating the number of unseen species: how many words did Shakespeare know?  Biometrika 63(3), 435–47.
  * Favaro et al. (2009) Favaro, S., A. Lijoi, R. H. Mena, and I. Prünster (2009).  Bayesian non-parametric inference for species variety with a two-parameter Poisson–Dirichlet process prior.  J. R. Statist. Soc. B 71(5), 993–1008.
  * Favaro et al. (2012) Favaro, S., A. Lijoi, and I. Prünster (2012).  A new estimator of the discovery probability.  Biometrics 68(4), 1188–96.
  * Ferguson (1973) Ferguson, T. S. (1973).  A Bayesian analysis of some nonparametric problems.  Ann. Statist. 1(2), 209–30.
  * Fisher et al. (1943) Fisher, R. A., A. S. Corbet, and C. B. Williams (1943).  The relation between the number of species and the number of individuals in a random sample of an animal population.  J. of Anim. Ecol. 12(1), 42–58.
  * Gao et al. (2007) Gao, Z., C.-H. Tseng, Z. Pei, and M. Blaser (2007, 03).  Molecular analysis of human forearm superfical skin bacterial biota.  Proc. Natl. Acad. Sci. U.S.A. 104, 2927–32.
  * Gnedin and Pitman (2005) Gnedin, A. and J. Pitman (2005).  Exchangeable Gibbs partitions and Stirling triangles.  Zapiski Nauchnykh Seminarov, POMI 325, 83–102.
  * Good (1953) Good, I. J. (1953, 12).  The population frequencies of species and the estimation of population parameters.  Biometrika 40(3-4), 237–64.
  * Good and Toulmin (1956) Good, I. J. and G. H. Toulmin (1956, 06).  The number of new species, and the increase in population coverage, when a sample is increased.  Biometrika 43(1-2), 45–63.
  * Hong (2013) Hong, Y. (2013).  On computing the distribution function for the Poisson binomial distribution.  Comput. Stat. Data Anal. 59, 41–51.
  * Hughes et al. (2001) Hughes, J. B., J. J. Hellmann, T. H. Ricketts, and B. J. M. Bohannan (2001).  Counting the uncountable: statistical approaches to estimating microbial diversity.  Appl. Environ. Microbiol. 67(10), 4399–406.
  * Ionita-Laza et al. (2009) Ionita-Laza, I., C. Lange, and N. M. Laird (2009).  Estimating the number of unseen variants in the human genome.  Proc. Natl. Acad. Sci. U.S.A. 106(13), 5008–13.
  * Korwar and Hollander (1973) Korwar, R. M. and M. Hollander (1973).  Contributions to the theory of Dirichlet processes.  Ann. Probab. 1(4), 705–11.
  * Lee et al. (2013) Lee, J., F. A. Quintana, P. Müller, and L. Trippa (2013).  Defining predictive probability functions for species sampling models.  Stat. Sci. 28(2), 209–22.
  * Lijoi et al. (2007) Lijoi, A., R. H. Mena, and I. Prünster (2007).  Bayesian nonparametric estimation of the probability of discovering new species.  Biometrika 94(4), 769–86.
  * Lijoi et al. (2020) Lijoi, A., I. Prünster, and T. Rigon (2020).  The Pitman–Yor multinomial process for mixture modeling.  Biometrika In press.
  * Mao (2004) Mao, C. X. (2004).  Predicting the conditional probability of discovering a new class.  J. Am. Statist. Assoc. 99(468), 1108–18.
  * Ovaskainen et al. (2020) Ovaskainen, O., N. Abrego, P. Somervuo, I. Palorinne, B. Hardwick, J.-M. Pitkänen, N. R. Andrew, P. A. Niklaus, N. M. Schmidt, S. Seibold, J. Vogt, E. V. Zakharov, P. D. N. Hebert, T. Roslin, and N. V. Ivanova (2020).  Monitoring fungal communities with the global spore sampling project.  Front. Ecol. Evol. 7, 511.
  * Perman et al. (1992) Perman, M., J. Pitman, and M. Yor (1992).  Size-biased sampling of Poisson point processes and excursions.  Prob. Theory Rel. Fields 92(1), 21–39.
  * Pitman (1996) Pitman, J. (1996).  Some developments of the Blackwell-Macqueen urn scheme.  In T. S. Ferguson, L. S. Shapley, and J. B. MacQueen (Eds.), Statistics, Probability and Game Theory. Papers in honor of David Blackwell, Volume 30 of IMS Lecture notes, Monograph Series, pp. 245–67. Hayward: Institute of Mathematical Statistics.
  * Pitman and Yor (1997) Pitman, J. and M. Yor (1997).  The two-parameter Poisson–Dirichlet distribution derived from a stable subordinator.  Ann. Prob. 25(2), 855–900.
  * Polson et al. (2013) Polson, N., J. Scott, and J. Windle (2013, 05).  Bayesian inference for logistic models using Pólya-gamma latent variables.  J. Am. Statist. Assoc. 108, 1339–49.
  * Quintana et al. (2020) Quintana, F. A., P. Müller, A. Jara, and S. N. MacEachern (2020).  The dependent Dirichlet process and related models.  arXiv:2007.06129.
  * Shen et al. (2003) Shen, T.-J., A. Chao, and C.-F. Lin (2003).  Predicting the number of new species in further taxonomic sampling.  Ecology 84(3), 798–804.
  * Spiegelhalter et al. (2002) Spiegelhalter, D. J., N. G. Best, B. P. Carlin, and A. Van Der Linde (2002).  Bayesian measures of model complexity and fit.  J. R. Statist. Soc. B 64(4), 583–639.
  * Thisted and Efron (1987) Thisted, R. and B. Efron (1987).  Did Shakespeare write a newly-discovered poem?  Biometrika 74(3), 445–55.
  * Zabell (1982) Zabell, S. L. (1982, 12).  W. E. Johnson’s "Sufficientness" postulate.  Ann. Statist. 10(4), 1090–99.


[affiliation=1,*]Chao-WeiHuang [affiliation=2,*]HuiLu
[affiliation=3]HongyuGong [affiliation=3]HirofumiInaguma
[affiliation=3]IliaKulikov [affiliation=3]RuslanMavlyutov
[affiliation=3]SravyaPopuri

# Investigating Decoder-only Large Language Models for Speech-to-text
Translation

###### Abstract

Large language models (LLMs), known for their exceptional reasoning
capabilities, generalizability, and fluency across diverse domains, present a
promising avenue for enhancing speech-related tasks. In this paper, we focus
on integrating decoder-only LLMs to the task of speech-to-text translation
(S2TT). We propose a decoder-only architecture that enables the LLM to
directly consume the encoded speech representation and generate the text
translation. Additionally, we investigate the effects of different parameter-
efficient fine-tuning techniques and task formulation. Our model achieves
state-of-the-art performance on CoVoST 2 and FLEURS among models trained
without proprietary data. We also conduct analyses to validate the design
choices of our proposed model and bring insights to the integration of LLMs to
S2TT.00footnotetext: ∗Work done during internship at Meta AI

###### keywords:

speech-to-text translation, large language models

##  1 Introduction

The task of speech-to-text translation (S2TT) involves converting audio
signals in one language into text in another, which is crucial for enabling
cross-lingual communication. Traditionally, S2TT has employed a cascaded
architecture with separate automatic speech recognition (ASR) and machine
translation (MT) components [1]. Recently, the emerging end-to-end (E2E)
approach, which integrates audio encoding and text decoding into a single
process, has gained popularity for the benefits of error propagation
mitigation and latency reduction [2, 3]. While it has achieved significant
performance improvement, S2TT still suffers from poor out-of-domain
generalization and failure to capture nuanced details, e.g., slangs and
cultural differences [4].

Large language models (LLMs) have emerged as powerful techniques for natural
language processing (NLP) due to their excellent reasoning capabilities and
generalizability. They excel at generating text for a wide range of tasks
based on large-scale pre-training [5, 6], instruction fine-tuning [7], and
reinforcement learning from human feedback [8, 9]. LLMs are also known for
their fluency and diverse domain coverage, which could potentially mitigate
the generalization gap for S2TT models. However, it is still under-explored as
to how LLMs should be integrated to improve S2TT performance.

In this paper, we aim to examine various aspects of adapting decoder-only LLMs
to S2TT, including architectural design, parameter-efficient fine-tuning, and
taks formulations. We propose a decoder-only architecture that directly
consumes continuous speech representation instead of discretized tokens. Our
proposed model achieves state-of-the-art S2TT performance without relying on
large amount of proprietary data. Furthermore, we analyze design choices of
each aspect of our experimental pipeline. Our contribution can be summarized
as the following:

  * •

We propose a decoder-only architecture for integrating LLMs to S2TT.

  * •

Our proposed model outperforms state-of-the-art S2TT models on CoVoST 2 and
FLEURS without training on proprietary data.

  * •

We conduct analyses to validate our design choices, which we hope could
facilitate future research on S2TT with LLMs.

Figure 1: Illustration of our proposed decoder-only architecture.

##  2 Related Work

###  2.1 Speech-to-text Translation

Speech-to-text translation has seen significant progress, especially for end-
to-end models. To solve the data scarcity issue of training end-to-end models,
multiple large-scale datasets have been collected, e.g., MuST-C [10], CoVoST
[11], Common Voice [12], and VoxPopuli [13]. Recent studies have started to
focus on multilingual S2TT, where a single end-to-end model supports multiple
translation directions [2]. The advent of pretrained models in language [5, 6]
and speech [14, 15] have facilitated new state-of-the-art models that
leveraged the pretrain-then-finetune paradigm [3, 16].

Our paper studies the integration of decoder-only LLMs to S2TT, which is still
under-explored due to their new architecture and emerging capabilities.

###  2.2 Speech and Audio LLMs

With the emergence of large language models, studies have explored applying
them to different modalities. LTU [17] fine-tuned LLMs on diverse audio
datasets, thus enabling LLMs to reason given audio inputs. Furthermore,
various works have explored extending the instruction-following capability of
LLMs to speech and audio inputs [18, 19, 20]. While these methods make it
possible for LLMs to handle a variety of speech and audio tasks, their
performance on individual tasks often falls short of that achieved by
specialized models.

Another line of research focuses on adapting LLMs to a specific speech or
audio task. Recent works have examined the integration of LLMs to automatic
speech recognition, demonstrating their potential in understanding the content
of speech [21, 22]. Similar to our work, AudioPaLM [23], Speech-LLaMA [24],
and SALM [25] aimed at leveraging LLMs to improve the state-of-the-art S2TT
performance. AudioPaLM proposed to adapt LLMs to speech by discretizing speech
representations and treat the discrete tokens as additional text tokens. Such
method has two drawbacks, as shown in the original paper: 1) its performance
is highly dependent on the quality of the speech encoder, and 2) the
discretization makes fine-tuning the speech encoder hard, which requires fine-
tuning the speech encoder with ASR first [23]. Our paper demonstrates that
using continuous speech representations mitigates these issues, achieving
better performance while being simpler. Speech-LLaMA and SALM both proposed
briding LLMs and speech encoders with a modality adaptor and fine-tunes LLMs
via LoRA [26]. Additionally, Speech-LLaMA introduced CTC compressor to shorten
the speech input. Our paper adopts a simpler length adaptor in our
architecture, and applies LNA fine-tuning [3] and demonstrates that it
outperforms LoRA significantly.

##  3 Our Method

In this section, we introduce the task formulations (§3.1), the architectural
designs of our model (§3.2), how the model is trained (§3.3), and parameter-
efficient fine-tuning techniques (§3.4).

###  3.1 Task Formulations

The task of speech-to-text translation is to translate the source speech input
$S$ into the corresponding target translation $Y=\\{y_{1},\cdots,y_{M}\\}$
which is in the target language. Following prior work [23], we define two
formulations of our S2TT model: 1) the standard formulation where the model
generates the target sequence directly $f\colon S\rightarrow Y$, and 2) the
chained formulation where the model first generates the transcription in the
source language then the translation in the target language
$f_{\text{chain}}\colon S\rightarrow\\{Y_{\text{ASR}},Y\\}$, where
$Y_{\text{ASR}}$ denotes the transcription of the source speech. It is also
common to include ASR during training as an auxiliary task, which is
formulated as $f_{\text{ASR}}\colon S\rightarrow Y_{\text{ASR}}$. Therefore,
we include $f$, $f_{\text{chain}}$, and $f_{\text{ASR}}$ during training for
multi-task training, and perform either $f$ or $f_{\text{chain}}$ during
inference.

###  3.2 Architecture

Our model consists of a speech encoder and a text decoder, both using the
Transformer architecture [27]. An illustration of the overall architecture is
shown in Figure 1.

Our speech encoder is based on W2v-BERT [15], a self-supervised pre-trained
speech encoder. For a given speech input $S$, we first convert the speech
signal to fbank features with 80 mel banks, a context window of 25 ms, and a
stride of 10 ms. The speech encoder $E_{s}$ encodes the fbank features
$F=\\{F_{1},\cdots,F_{n}\\}$ to their corresponding hidden representations
$E_{s}(F)$, where $n$ denotes the sequence length of the fbank features.
Speech frames are typically much more granular than text tokens. Therefore, we
employ a length adapter on top of the speech encoder to reduce the length of
the speech representations. The length adapter consists of a single
1-dimensional convolutional layer with a filter size and stride of $k$, which
reduces the length of the speech representations by $k$-fold.

The text decoder is based on LLaMA-2 [9], a decoder-only large language model
pre-trained on 2 trillion text tokens with a language modeling objective. The
speech inputs and text inputs are encoded with their corresponding encoders,
i.e., speech encoder for speech inputs and text embedding layer for text
inputs. Subsequently, the encoded representations are concatenated and fed to
the transformer decoder. In other words, we treat the encoded speech
representations $S$ the same as the text embeddings, without discretizing them
as done in prior work [23]. A triangular mask is appied to the self-attention
layers to restrict tokens from atteding to latter positions. More formally,
given an interleaving sequence of text and speech sequences
$X=\\{X^{1},F,X^{2}\\}$, where $X^{i}=\\{x^{i}_{i},\cdots,x^{i}_{|x^{i}|}\\}$
denotes a text sequence, $X^{1}$ denotes the prefix text, and $X^{2}$ denotes
the suffix text. After encoding, the input sequence to the transformer decoder
will be $\mathbf{X}=\\{\text{Emb}(X^{1}),E_{s}(F),\text{Emb}(X^{2})\\}$, where
Emb denotes the text embedding layer. Note that we flatten the sequences in
$\mathbf{X}$ before processing them with the decoder. Finally, we apply a
linear transformation to the decoder outputs to obtain the logits for
predicting the next token $\mathbf{O}=W^{\top}D(\mathbf{X})$, where $D$
denotes the transformer decoder and $W\in\mathbb{R}^{h\times|V|}$ is a
trainable matrix where $|V|$ denotes the vocabulary size.

Ar Ca Cy De Es Et Fa Fr Id It Ja Lv Mn Nl Pt Ru Sl Sv Ta Tr Zh Avg Trained
with Proprietary Data Whisper-large [28] 39.7 31.8 18.0 21.5 36.3 15.0 36.4
48.1 30.9 26.1 0.1 13.9 41.2 19.3 51.6 43.3 21.6 40.1 42.9 4.2 28.3 29.1 USM-M
[29] - - - - - - - - - - - - - - - - - - - - - 30.7 Speech-LLaMA [24] 28.2 - -
27.1 27.9 18.7 - 25.2 - 25.9 19.9 - - 36.5 32.0 36.8 22.7 29.0 - - 12.3 -
AudioPaLM [23] 48.7 38.4 25.5 13.7 43.4 30.0 44.8 56.2 44.3 25.9 7.6 35.0 48.3
29.4 57.3 55.6 42.6 44.2 53.3 9.0 41.0 37.8 Trained with Public Data Only
XLS-R [16] 17.1 33.8 9.4 14.0 33.6 11.1 37.6 16.5 34.9 3.5 1.6 19.5 31.7 12.9
41.8 39.5 19.6 39.2 29.6 0.5 16.7 22.1 ComSL-large [30] - - - - - - - - - - -
- - - - - - - - - - 31.5 AudioPaLM† [23] - - - - - - - - - - - - - - - - - - -
- - 33.1 W2vBERT+NLLB 42.0 38.4 18.3 52.0 39.3 23.6 41.3 47.3 39.4 18.1 3.5
18.4 43.0 27.2 50.8 51.7 36.9 42.2 40.6 6.2 33.2 34.0 Ours 45.8 39.5 22.4 56.9
41.2 20.4 44.5 54.5 42.9 24.4 0.9 21.9 46.8 26.3 56.1 53.3 42.7 45.1 53.7 5.3
34.4 37.1

Table 1: Main results on the X-En test sets of CoVoST 2 (%). We report corpus
BLEU scores computed with SacreBLEU. The best results among models trained
with public data are bolded. †The result reported in the AudioPaLM paper [23]
when trained on only public datasets.

###  3.3 Training

As described above, we include three formulations, i.e., $f$,
$f_{\text{chain}}$, and $f_{\text{ASR}}$, for multi-task training. To let our
model distinguish among tasks, we provide different instructions in natural
language for each task $t$. The instructions include a description of the
task, the source language, and the target language. We format the instruction
$I$ and the source speech $S$ into the input sequence $X$ with a template. The
target sequence for training is formatted as:

$Y^{\prime}=\begin{cases}\text{Translation: }Y&\text{if }t=f\\\
\text{Transcription: }Y_{\text{ASR}}&\text{if }t=f_{\text{ASR}}\\\
\text{Transcription: }Y_{\text{ASR}}\text{ Translation: }Y&\text{if
}t=f_{\text{chain}}.\end{cases}$

Given a source speech $S$, an instruction $I$, and the formatted target
sequence $Y^{\prime}$, the training objective is to minimize the S2TT loss:

$\mathcal{L}(S,Y^{\prime})=-\frac{1}{M^{\prime}}\sum_{i=1}^{M^{\prime}}\text{log}P(y^{\prime}_{i}\mid
S,I,Y^{\prime}_{<i})$

where $M^{\prime}$ denotes the length of $Y^{\prime}$ and
$P(y^{\prime}_{i}\mid S,I,Y^{\prime}_{<i})$ denotes the probability of
$y^{\prime}_{i}$ predicted by the model given the source speech and the prior
tokens $Y^{\prime}_{<i}$ in the target sequence. The predicted probability is
obtained by applying the softmax function to the logits $\mathbf{O}$.

###  3.4 Parameter-efficient Fine-tuning

Large language models have billions of parameters, making it computationally
expensive and inefficient to fine-tune all of the parameters during training.
It is common to apply parameter-efficient fine-tuning techniques when fine-
tuning LLMs on downstream tasks to improve efficiency and mitigate
catastrophic forgetting. To this end, we employ and compare two parameter-
efficient fine-tuning techniques in this paper: LNA fine-tuning [3] and Low
Rank Adaptation (LoRA) [26].

####  3.4.1 LNA Fine-tuning

LayerNorm and Attention (LNA) fine-tuning adapts pretrained language and
speech models to S2TT by fine-tuning only the layer normalization and the
multi-head attention layers [3]. This method greatly reduces the number of
trainable parameters during fine-tuning and avoids catastrophic forgetting,
thus improving the downstream performance for multilingual speech-to-text
translation. Since the pretrained language model we use is a decoder-only
transformer model, we apply LNA fine-tuning and fine-tune only the layer
normalization and the self-attention layers in the transformer decoder.

####  3.4.2 Low Rank Adaptation (LoRA)

LoRA injects trainable rank decomposition matrices into the projections layers
of a transformer model, which serves as a residual path in addition to a
projection layer. During fine-tuning, only the decomposition matrices are
updated, while all of the pretrained parameters are frozen. Thus, the number
of trainable parameters is significantly reduced. The decomposition matrices
can be merged into the original projection matrix after fine-tuning.
Therefore, there is no additional computation nor additional parameters
compared to the pretrained transformer model during inference, making LoRA a
common technique for adapting large language models efficiently.

##  4 Experiments

###  4.1 Experimental Setup

We train and evaluate our models on publicly available datasets. For training,
we use CoVoST2 [11], Common Voice 11 [12], and VoxPopuli [13] datasets.
CoVoST-2 is a speech-to-text translation dataset consisting of 21 languages.
The dataset includes human-labeled translation pairs from 21 languages to
English (X-En), and from English to 15 languages (En-X). Common Voice is a
collection of speech-text pairs where the speech was recorded by annotators
given the text transcription. VoxPopuli consists of speech from the European
Parliament with the corresponding transcriptions and interpretations in 15
languages.

We conduct in-domain evaluation on the test sets of CoVoST 2. Additionally, we
perform zero-shot evaluation on FLEURS [4], a dataset that aims to evaluate
the out-of-domain generalizability of speech translation models. Note that for
all datasets, we only use the directions that are present in CoVoST2. We
report BLEU scores from SacreBLEU and additionally the model-based COMET score
with the model wmt22-comet-da [31].

###  4.2 Implementation Details

We employ a pretrained W2v-BERT [15] model that was released in [32] with 600M
parameters that is pretrained on 4 million hours of speech data with a self-
supervised objective as the speech encoder. The text decoder is initialized
with LLaMA2-7B-chat [9].

We implement our models, training, and evaluation procedures with the Fairseq2
library111https://github.com/facebookresearch/fairseq2. During training, the
effective batch size is set to 800K speech frames, or 8000 seconds of speech
inputs. We optimize the model with the AdamW optimizer and set the learning
rate to 1e-4. The learning rate is warmed up for 5000 steps and linearly
decayed until the maximum number of steps is reached, which is set to 60000.
We fine-tune all parameters of the speech encoder and apply parameter-
efficient fine-tuning methods to the text decoder. All experiments are
conducted on 32 NVIDIA A100 GPUs.

CoVoST 2 FLEURS BLEU COMET BLEU COMET Encoder-Decoder W2vBERT+NLLB 33.8 80.1
22.7 76.4 W2vBERT+LLaMA2 33.3 79.7 18.2 74.2 Decoder-only Ours 36.3 81.4 23.4
77.6

Table 2: Results of different architectures (%). We report the average BLEU
score and COMET score on the 21 X-En directions on CoVoST 2 and FLEURS.

Rank Layers CoVoST 2 FLEURS Freeze W2vBERT - - 15.6 3.8 Freeze decoder - -
29.9 19.0 LoRA 8 q, v 32.1 20.8 32 q, k, v, o 32.7 21.1 32 All linear 33.0
21.4 LNA - - 36.3 23.4

Table 3: Results of different parameter-efficient fine-tuning methods (%).
Rank and Layers refer to the configuration of LoRA. The notations q, k, v, o
denote the query, key, value, output layers of the self-attention layers
respectively.

###  4.3 Baseline Methods

We compare our model with various state-of-the-art baselines that were trained
on the same set of public datasets as our method, i.e., CoVoST 2, Common
Voice, and VoxPopuli. XLS-R [16] is a self-supervised cross-lingual speech
representation model. ComSL [30] conducts self-training on the Common Voice
dataset. Additionally, we implement an encoder-decoder baseline with W2vBERT
as the speech encoder and NLLB [33] 1.3B as the text decoder.

We also compare our model with models trained with proprietary data. Whisper
[28] trains a robust speech recognition and translation model with large
amounts of weak supervisions. USM [29] is an universal speech model pretrained
with 12 million hours of speech data. Speech-LLaMA [24] shares a similar
architecture with our model and was trained with in-house data and LoRA [26].
AudioPaLM [23] is the state-of-the-art method on CoVoST2 which is trained on
proprietary data. We also include a variant of AudioPaLM that is trained on
public datasets only for a fair comparison, which is reported in the paper
[23].

###  4.4 Results

The main results on CoVoST 2 are reported in Table 1. Our model achieves an
average BLEU score of 37.1, which is the new state-of-the-are performance
among models trained with public data only. Notably, our model outperforms the
AudioPaLM variant which was trained on only public datasets, demonstrating the
superiority of our proposed method. When compared to models trained with
proprietary data, our model outperforms all of them and achieves comparable
performance to AudioPaLM. These results demonstrate that our method integrates
LLMs to S2TT efficiently and effectively.

| CoVoST 2 | FLEURS
---|---|---
Ours | 37.1 | 23.4
\- $f_{\text{chain}}$ | 35.4 | 22.4
\- $f_{\text{ASR}}$ | 36.4 | 22.7
\- $f_{\text{ASR}}\mathbin{\&}f_{\text{chain}}$ | 35.8 | 22.5
Table 4: Results of formulation ablation (%).

##  5 Discussion

In this section, we conduct various experiments to analyze and discuss the
details of our proposed method.

###  5.1 Architectural Design

With decoder-only LLMs, it is unclear as to which architecture performs the
best for S2TT. We compare our decoder-only architecture with encoder-decoder
models, with NLLB [33] and LLaMA-2 [9] as the text decoder. As shown in Table
2, our model significantly outperforms the encoder-decoder counterpart on both
CoVoST 2 and FLEURS. Furthermore, encoder-decoder with LLaMA 2 even
underperforms NLLB, demonstrating that encoder-decoder architecture are
unsuitable for decoder-only LLMs. We hypothesize that it is the newly
introduced encoder-decoder attention layers which are not pretrained that
degrade the performance of encoder-decoder models.

###  5.2 Parameter-efficient Fine-tuning

We compare LNA fine-tuning, LoRA, and the effect of freezing pretrained
models. As shown in Table 3, LNA fine-tuning significantly outperforms LoRA
with various configurations. This result suggests that adopting LoRA, as done
in prior work such as Speech-LLaMA [24], is suboptimal for S2TT. Freezing the
text decoder during fine-tuning yields even worse performance than LoRA,
demonstrating the importance of fine-tuning the text decoder. Finally,
freezing the speech encoder results in detrimental performance degradation.
This result shows that fine-tuning the speech encoder is crucial for aligning
the speech representation with the text inputs. We hypothesize that this leads
to the underperformance of AudioPaLM with encoders that are not fine-tuned
with ASR [23], since the discretization of speech representations makes fine-
tuning the speech encoder non-trivial.

###  5.3 Ablation of Formulations

Table 4 shows the results of various combination of the formulations. Removing
either $f_{\text{ASR}}$ or $f_{\text{chain}}$ degrades the S2TT performance.
Notably, training with $f$ and $f_{\text{ASR}}$ slightly underperforms $f$,
showing that multi-task training with ASR does not always improve performance.

##  6 Conclusion

In this paper, we propose a decoder-only architecture that adapts a decoder-
only LLM to the speech-to-text translation task. Our proposed method is simple
and effective, achieving state-of-the-art performance and is comparable to the
best-performing proprietary model. We conduct additional analyses to examine
the effect of different design choices regarding architectural design,
parameter-efficient fine-tuning, and task formulations. We hope that our
findings could facilitate future work on leveraging LLMs in the S2TT task.

## References

  * [1] H. Ney, ``Speech translation: Coupling of recognition and translation,'' in _1999 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings. ICASSP99 (Cat. No. 99CH36258)_ , vol. 1. IEEE, 1999, pp. 517–520.
  * [2] H. Inaguma, K. Duh, T. Kawahara, and S. Watanabe, ``Multilingual end-to-end speech translation,'' in _2019 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2019, pp. 570–577.
  * [3] X. Li, C. Wang, Y. Tang, C. Tran, Y. Tang, J. Pino, A. Baevski, A. Conneau, and M. Auli, ``Multilingual speech translation from efficient finetuning of pretrained models,'' in _Proceedings of ACL-IJCNLP 2021_ , 2021, pp. 827–838.
  * [4] A. Conneau, M. Ma, S. Khanuja, Y. Zhang, V. Axelrod, S. Dalmia, J. Riesa, C. Rivera, and A. Bapna, ``Fleurs: Few-shot learning evaluation of universal representations of speech,'' in _2022 IEEE Spoken Language Technology Workshop (SLT)_. IEEE, 2023, pp. 798–805.
  * [5] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, ``BERT: Pre-training of deep bidirectional transformers for language understanding,'' in _Proceedings NAACL-HLT 2019_ , 2019, pp. 4171–4186.
  * [6] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell _et al._ , ``Language models are few-shot learners,'' _Advances in neural information processing systems_ , vol. 33, pp. 1877–1901, 2020.
  * [7] J. Wei, M. Bosma, V. Y. Zhao, K. Guu, A. W. Yu, B. Lester, N. Du, A. M. Dai, and Q. V. Le, ``Finetuned language models are zero-shot learners,'' _arXiv:2109.01652_ , 2021.
  * [8] L. Ouyang, J. Wu, X. Jiang, D. Almeida, C. Wainwright, P. Mishkin, C. Zhang, S. Agarwal, K. Slama, A. Ray _et al._ , ``Training language models to follow instructions with human feedback,'' _Advances in Neural Information Processing Systems_ , vol. 35, pp. 27 730–27 744, 2022.
  * [9] H. Touvron, L. Martin, K. Stone, P. Albert, A. Almahairi, Y. Babaei, N. Bashlykov, S. Batra, P. Bhargava, S. Bhosale _et al._ , ``Llama 2: Open foundation and fine-tuned chat models,'' _arXiv:2307.09288_ , 2023.
  * [10] M. A. Di Gangi, R. Cattoni, L. Bentivogli, M. Negri, and M. Turchi, ``MuST-C: a Multilingual Speech Translation Corpus,'' in _Proceedings of NAACL-HLT 2019_ , 2019, pp. 2012–2017.
  * [11] C. Wang, A. Wu, and J. Pino, ``Covost 2 and massively multilingual speech-to-text translation,'' _arXiv:2007.10310_ , 2020.
  * [12] R. Ardila, M. Branson, K. Davis, M. Kohler, J. Meyer, M. Henretty, R. Morais, L. Saunders, F. Tyers, and G. Weber, ``Common voice: A massively-multilingual speech corpus,'' in _Proceedings of LREC 2020_ , 2020, pp. 4218–4222.
  * [13] C. Wang, M. Riviere, A. Lee, A. Wu, C. Talnikar, D. Haziza, M. Williamson, J. Pino, and E. Dupoux, ``VoxPopuli: A large-scale multilingual speech corpus for representation learning, semi-supervised learning and interpretation,'' in _Proceedings of ACL-IJCNLP 2021_ , 2021, pp. 993–1003.
  * [14] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, ``wav2vec 2.0: A framework for self-supervised learning of speech representations,'' _Advances in neural information processing systems_ , vol. 33, pp. 12 449–12 460, 2020.
  * [15] Y.-A. Chung, Y. Zhang, W. Han, C.-C. Chiu, J. Qin, R. Pang, and Y. Wu, ``W2v-bert: Combining contrastive learning and masked language modeling for self-supervised speech pre-training,'' in _2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2021, pp. 244–250.
  * [16] A. Babu, C. Wang, A. Tjandra, K. Lakhotia, Q. Xu, N. Goyal, K. Singh, P. von Platen, Y. Saraf, J. Pino, A. Baevski, A. Conneau, and M. Auli, ``XLS-R: Self-supervised Cross-lingual Speech Representation Learning at Scale,'' in _Proc. Interspeech 2022_ , 2022, pp. 2278–2282.
  * [17] Y. Gong, H. Luo, A. H. Liu, L. Karlinsky, and J. R. Glass, ``Listen, think, and understand,'' in _The Twelfth International Conference on Learning Representations_ , 2024.
  * [18] M. Wang, W. Han, I. Shafran, Z. Wu, C.-C. Chiu, Y. Cao, N. Chen, Y. Zhang, H. Soltau, P. K. Rubenstein _et al._ , ``Slm: Bridge the thin gap between speech and text foundation models,'' in _ASRU 2023_. IEEE, 2023, pp. 1–8.
  * [19] Y. Chu, J. Xu, X. Zhou, Q. Yang, S. Zhang, Z. Yan, C. Zhou, and J. Zhou, ``Qwen-audio: Advancing universal audio understanding via unified large-scale audio-language models,'' _arXiv:2311.07919_ , 2023.
  * [20] C. Tang, W. Yu, G. Sun, X. Chen, T. Tan, W. Li, L. Lu, Z. Ma, and C. Zhang, ``Salmonn: Towards generic hearing abilities for large language models,'' _arXiv:2310.13289_ , 2023.
  * [21] Y. Fathullah, C. Wu, E. Lakomkin, J. Jia, Y. Shangguan, K. Li, J. Guo, W. Xiong, J. Mahadeokar, O. Kalinli _et al._ , ``Prompting large language models with speech recognition abilities,'' _arXiv:2307.11795_ , 2023.
  * [22] W. Yu, C. Tang, G. Sun, X. Chen, T. Tan, W. Li, L. Lu, Z. Ma, and C. Zhang, ``Connecting speech encoder and large language model for asr,'' _arXiv:2309.13963_ , 2023.
  * [23] P. K. Rubenstein, C. Asawaroengchai, D. D. Nguyen, A. Bapna, Z. Borsos, F. d. C. Quitry, P. Chen, D. E. Badawy, W. Han, E. Kharitonov _et al._ , ``Audiopalm: A large language model that can speak and listen,'' _arXiv:2306.12925_ , 2023.
  * [24] J. Wu, Y. Gaur, Z. Chen, L. Zhou, Y. Zhu, T. Wang, J. Li, S. Liu, B. Ren, L. Liu _et al._ , ``On decoder-only architecture for speech-to-text and large language model integration,'' in _2023 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2023, pp. 1–8.
  * [25] Z. Chen, H. Huang, A. Andrusenko, O. Hrinchuk, K. C. Puvvada, J. Li, S. Ghosh, J. Balam, and B. Ginsburg, ``Salm: Speech-augmented language model with in-context learning for speech recognition and translation,'' _arXiv:2310.09424_ , 2023.
  * [26] E. J. Hu, Y. Shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, and W. Chen, ``LoRA: Low-rank adaptation of large language models,'' in _International Conference on Learning Representations_ , 2022.
  * [27] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, ``Attention is all you need,'' _Advances in neural information processing systems_ , vol. 30, 2017.
  * [28] A. Radford, J. W. Kim, T. Xu, G. Brockman, C. McLeavey, and I. Sutskever, ``Robust speech recognition via large-scale weak supervision,'' in _International Conference on Machine Learning_. PMLR, 2023, pp. 28 492–28 518.
  * [29] Y. Zhang, W. Han, J. Qin, Y. Wang, A. Bapna, Z. Chen, N. Chen, B. Li, V. Axelrod, G. Wang _et al._ , ``Google usm: Scaling automatic speech recognition beyond 100 languages,'' _arXiv:2303.01037_ , 2023.
  * [30] C. Le, Y. Qian, L. Zhou, S. Liu, Y. Qian, M. Zeng, and X. Huang, ``Comsl: A composite speech-language model for end-to-end speech-to-text translation,'' _Advances in Neural Information Processing Systems_ , vol. 36, 2024.
  * [31] R. Rei, J. G. C. de Souza, D. Alves, C. Zerva, A. C. Farinha, T. Glushkova, A. Lavie, L. Coheur, and A. F. T. Martins, ``COMET-22: Unbabel-IST 2022 submission for the metrics shared task,'' in _Proceedings of the Seventh Conference on Machine Translation (WMT)_ , 2022, pp. 578–585.
  * [32] L. Barrault, Y.-A. Chung, M. C. Meglioli, D. Dale, N. Dong, P.-A. Duquenne, H. Elsahar, H. Gong, K. Heffernan, J. Hoffman _et al._ , ``Seamlessm4t-massively multilingual & multimodal machine translation,'' _arXiv:2308.11596_ , 2023.
  * [33] M. R. Costa-jussà, J. Cross, O. Çelebi, M. Elbayad, K. Heafield, K. Heffernan, E. Kalbassi, J. Lam, D. Licht, J. Maillard _et al._ , ``No language left behind: Scaling human-centered machine translation,'' _arXiv:2207.04672_ , 2022.