# A Bayesian Nonparametric Approach to Species Sampling Problems with Ordering

Cecilia Balocchi<EMAIL_ADDRESS>Department of Economics and
Statistics, University of Torino, Torino, Italy.  Federico Camerlenghi Also
affiliated to Collegio Carlo Alberto, Piazza V. Arbarello 8, Torino, and
BIDSA, Bocconi University, Milano, Italy<EMAIL_ADDRESS>Department of Economics, Management and Statistics, University of Milano-
Bicocca, Milano, Italy.  Stefano Favaro Also affiliated to Collegio Carlo
Alberto, Piazza V. Arbarello 8, Torino, and IMATI-CNR “Enrico Magenes”, Milan,
Italy<EMAIL_ADDRESS>Department of Economics and Statistics,
University of Torino, Torino, Italy.

###### Abstract

Species-sampling problems (SSPs) refer to a vast class of statistical problems
that, given an observable sample from an unknown population of individuals
belonging to some species, call for estimating (discrete) functionals of the
unknown species composition of additional unobservable samples. A common
feature of SSPs is the invariance with respect to species labelling, i.e.
species’ labels are immaterial in defining the functional of interest, which
is at the core of the development of the Bayesian nonparametric (BNP) approach
to SSPs under the popular Pitman-Yor process (PYP) prior. In this paper, we
consider SSPs that are not “invariant” to species labelling, in the sense that
an ordering or ranking is assigned to species’ labels, and we develop a BNP
approach to such problems. In particular, inspired by the population genetics
literature on age-ordered alleles’ compositions, with a renowned interest in
the frequency of the oldest allele, we study the following SSP with ordering:
given an observable sample from unknown population of individuals belonging to
some species (alleles), with species’ labels being ordered according to
weights (ages), estimate the frequencies of the first $r$ order species’
labels in an enlarged sample obtained by including additional unobservable
samples. Our BNP approach relies on an ordered version of the PYP prior, which
leads to an explicit posterior distribution of the first $r$ order
frequencies, with corresponding estimates being simple and computationally
efficient. We apply our approach to the analysis of genetic variation, showing
its effectiveness in the estimation of the frequency of the oldest allele, and
then discuss other applications in the contexts of citations to academic
articles and online purchases of items.

Keywords: Bayesian nonparametrics; exchangeable partition probability
function; first $r$ order frequency; ordered Pitman-Yor process prior; species
sampling problems; population genetics.

##  1 Introduction

Species sampling problems (SSPs) refer to a vast class of statistical
problems, of which the estimation of the number of unseen species is arguably
the most popular example (Good and Toulmin, 1956; Efron and Thisted, 1976;
Bunge and Fitzpatrick, 1993; Lijoi et al., 2007; Orlitsky et al., 2016).
Consider $n\geq 1$ observable samples from a generic population of
individuals, with each individual taking a value in a (possibly infinite)
discrete space of symbols or species. The unseen-species problem assumes that
observable samples are modeled as a random sample $(X_{1},\ldots,X_{n})$ from
an unknown discrete distribution $p$, and calls for estimating

$K_{m}^{(n)}=|\\{X_{n+1},\ldots,X_{n+m}\\}\setminus\\{X_{1},\ldots,X_{n}\\}|,$
(1)

namely the number of hitherto unseen symbols that would be observed if $m\geq
1$ additional samples $(X_{n+1},\ldots,X_{n+m})$ were collected from the same
distribution. In general, SSPs comprise of generalizations or refinements of
the unseen-species problem, calling for the estimation of (discrete)
functionals of the species’ composition of unobservable samples, e.g. missing
mass, discovery probabilities, unseen species with prevalences and coverages
of prevalence. We refer to Deng et al. (2019) and Balocchi et al. (2022) for
up-to-date reviews of SSPs, both in methods and applications, mostly in the
field of biological sciences but also in machine learning, electrical
engineering, theoretical computer science and information theory.

A common feature of SSPs is that species’ labels identifying the $X_{i}$’s are
immaterial in the definition of the functional of interest, as for instance in
(1), thus making SSPs “invariant” to species labelling. Such a feature is at
the core of the development of the Bayesian nonparametric (BNP) approach to
SSPs (Lijoi et al., 2007, 2008; Favaro et al., 2009, 2013), which relies on
the specification of a (nonparametric) prior $\mathscr{P}$ for the unknown
distribution $p$, i.e.

$\displaystyle X_{i}\,|\,P$
$\displaystyle\quad\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\quad
P\qquad i=1,\ldots,n,$ (2) $\displaystyle P$
$\displaystyle\quad\sim\quad\mathscr{P}.$

Species sampling models (SSMs) (Pitman, 1996) provide a natural choice for
$\mathscr{P}$, including the Dirichlet process (DP) prior (Ferguson, 1973) and
the Pitman-Yor process (PYP) prior (Pitman and Yor, 1997). Under the BNP model
(2) with a SSM for $P$, the random sample $(X_{1},\ldots,X_{n})$ induces a
random partition $\tilde{\pi}_{n}$ of $[n]=\\{1,\ldots,n\\}$ whose blocks are
the (equivalence) classes of the equivalence relation $i\sim j\iff
X_{i}=X_{j}$ almost surely. In particular, $\tilde{\pi}_{n}$ is exchangeable
(Pitman, 2006, Chapter 2), namely its distribution is such that the
probability of any partition of $[n]$ with $k$ blocks of frequencies
$(n_{1},\ldots,n_{k})$ is a symmetric function of compositions
$(n_{1},\ldots,n_{k})$ of $[n]$. The exchangeability of $\tilde{\pi}_{n}$
implies that blocks’ labels are immaterial, and therefore it legitimates the
BNP approach to SSPs under SSMs.

###  1.1 Our contributions

In this paper, we consider SSPs that are not “invariant” with respect to
species labelling, in the sense that an ordering or ranking is assigned to
species’ labels, and we develop a BNP approach to such problems. Under the
infinitely-many neutral alleles model for the evolution of genetic populations
(Ewens, 1972; Kingman, 1975; Watterson, 1977; Griffiths, 1979), the work of
Donelly and Tavaré (1986) first investigated the alleles’ composition of a
random sample from the population by also taking into account the ages of
alleles, namely the times elapsed since the first time each allele first
appeared in the sample. This study led to an age-ordered version of the random
partition induced by the DP prior, where species are alleles and species’
labels are ordered according to the age of alleles in such a way that the
smaller the order the older the allele. Besides providing distributional
properties of the age-ordered random partition, Donelly and Tavaré (1986)
applied such a model to answer a critical question raised in Crow (1972): “Is
the most frequent allele the oldest?”. Under the infinitely-many neutral
alleles model, Donelly and Tavaré (1986) came up with a positive answer to
such a question, showing that the probability that an allele represented $i$
times in sample of size $n$ is the oldest is $i/n$ (Kelly, 1977; Watterson,
1977). The question of Crow (1972) is to some extent a SSP with ordering, as
the object of interest involves a species’ label with a precise order, namely
the species’ label of order $1$ that corresponds to the oldest allele.

Inspired by the seminal work of Donelly and Tavaré (1986), we study the
following SSP with ordering: assuming $n$ observable samples to be modeled as
a random sample $((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ from an unknown
discrete distribution $q$, with the $T_{i}$’s being weights inducing an
ordering among species’ labels identifying the $X_{i}$’s, estimate the
frequencies of the first $r$ order species’ labels in an enlarged sample
obtained by including $m$ additional unobservable samples
$((T_{n+1},X_{n+1}),\ldots,(T_{n+m},X_{n+m}))$ from the same distribution.
Within the population genetic setting of Donelly and Tavaré (1986), this SSP
corresponds to estimate the frequency of the $r$ oldest alleles in a sample of
size $(n+m)$ based on $n$ observable samples. We introduce a BNP approach to
estimate the first $r$ order frequencies, which relies on the use of spatial
neutral to the right SSMs (James, 2006), or ordered SSMs, as prior
distributions for $q$. The most popular ordered SSM is the ordered DP prior,
which is known to induce the age-ordered random partition of Donelly and
Tavaré (1986) (Gnedin and Pitman, 2005). Here, we consider the more general
ordered PYP prior (Gnedin and Pitman, 2005), and we determine the posterior
distribution of the first $r$ order frequencies; then, a BNP estimator is
proposed in terms of the posterior mean, which results to be of easy
implementation and computationally efficient. Of special interest is the case
$r=1$ that, in the setting of Donelly and Tavaré (1986), leads to an estimate
of the frequency of the oldest allele.

We present an empirical validation of the effectiveness of our BNP approach,
both on synthetic and real data. It is natural for us to focus on applications
to genetic data, with the weights $T_{i}$’s being interpreted as the ages of
the alleles $X_{i}$’s. Modeling the interplay between the alleles’ composition
of a genetic population and the age of alleles dates back to the 1970s and the
1980s, and nowadays the genealogical structure of alleles is well recognized
as a fundamental aspect in many inferential (decision) processes of population
genetics. In particular, studying genetic variation while incorporating the
information on the variants’ age allows to enhance the investigation of
several problems, such as analyzing and comparing populations structure,
detecting which samples are related, studying demographic history, and
learning about genetic susceptibility to disease (Mathieson and McVean, 2014).
For example, it enables researchers to use variants’ age distribution to
compare populations, to differentiate age distributions in pathogenic and
benign variants, and to learn about genealogical history (Albers and McVean,
2020). Here, we apply our BNP methodology to the problem of estimating the
frequency of the oldest allele, using genetic variation data from the 1000
Genomes Project (1000 Genomes Project Consortium, 2015) and variants’ age
estimates from the Human Genome Dating Project (Albers and McVean, 2020).

Besides population genetics, we discuss other two contexts where SSPs with
ordering arise naturally: citations to academic articles and online purchases
of items. In the context of citations to academic articles, with articles
being ordered according to their publication’s dates, one may be interested in
the frequency of citations to the oldest paper. Citation data are often
studied in the framework of citation networks, thus allowing to study the
movement of ideas in academic fields or examine scholars’ influence (Portenoy
et al., 2017). Incorporating knowledge of articles’ age permits the
investigation of the effects of time on the number of citations (Hajra and
Sen, 2005) and to answer questions such as “are older papers more frequently
cited than newer papers?”. Assuming that each citation represents an
observation, that the article cited represents a species, and that the order
of the cited article is determined by its publication’s date, we may apply our
BNP approach to predict the frequency of citations to the oldest article in
future observations. In the context of online purchases of items, species are
represented by the items purchased, with items being ordered according to
their costs, as well as by a generic independent measure of popularity. In
particular, we may apply our BNP approach to study of the distribution of the
most popular item in future purchases, which is particularly relevant in order
to plan changes in marketing strategies.

###  1.2 Organization of the paper

The paper is structured as follows. In Section 2 we present the ordered PYP
prior, and review its sampling structure in terms of the induced sampling
formulae and predictive distributions. In Section 3 we provide the posterior
distribution of the first $r$ order frequencies, with emphasis on the special
case $r=1$, and obtain corresponding estimators. Section 4 contains numerical
illustrations of our BNP approach, both on synthetic and real data, whereas in
Section 5 we discuss our work and some directions for future research.
Additional numerical illustrations on genetic data, an illustration in the
context of citations to academic articles, and the proofs of our results are
deferred to the Supplementary Material.

##  2 The ordered PYP

To introduce the ordered PYP, it is useful to recall the PYP and its sampling
structure. Let $P$ be a PYP with parameter $\alpha\in[0,1)$ and
$\theta>-\alpha$ on a measurable space $\mathbb{X}$. That is $P=\sum_{i\geq
1}P_{i}\delta_{S_{i}}$, where: i) $P_{1}=V_{1}$ and $P_{i}=V_{i}\prod_{1\leq
j\leq i-1}(1-V_{j})$ with $(V_{i})_{i\geq 1}$ being independent Beta random
variables with parameter $(1-\alpha,\theta+i\alpha)$; ii) $(S_{i})_{i\geq 1}$
be random variables, independent of the $V_{i}$’s, and independent and
identically distributed according to a non-atomic distribution $\nu$ on
$\mathbb{X}$ (Perman et al., 1992; Pitman, 1995). Because of the (almost sure)
discreteness of $P$, a random sample $(X_{1},\ldots,X_{n})$ from $P$ induces a
random partition $\tilde{\pi}_{n}$ of $[n]$ into $K_{n}\leq n$ blocks,
labelled by $\\{X_{1}^{\ast},\ldots,X_{K_{n}}^{\ast}\\}$, with corresponding
frequencies $N_{j,n}=|i\in[n]\text{ : }X_{i}=X_{j}^{\ast}|$ for
$j=1,\ldots,K_{n}$ and such that $N_{j,n}\geq 1$ and $\sum_{1\leq j\leq
K_{n}}N_{j,n}=n$. In particular, if we set $(a)_{(r)}=\prod_{0\leq i\leq
r-1}(a+i)$ for any $a\geq 0$ and $r\in\mathbb{N}_{0}$, then the probability of
any partition of $[n]$ with $k$ blocks of frequencies $(n_{1},\ldots,n_{k})$
is

$\Pi_{k}^{(n)}(n_{1},\ldots,n_{k})=\frac{\prod_{i=1}^{k}(\theta+(i-1)\alpha)}{(\theta)_{(n)}}\prod_{i=1}^{k}(1-\alpha)_{(n_{i}-1)}.$
(3)

Equation (3) is referred to as the exchangeable partition probability function
(EPPF), a concept introduced in Pitman (1995) as a development of earlier
results in Kingman (1978). For $\alpha=0$ the PYP reduces to DP, and hence (3)
reduces to the Ewens sampling formula (Ewens, 1972). See Pitman (2006, Chapter
3 and Chapter 4) for a detailed account on the EPPF (3).

The predictive distribution of the PYP provides a generative scheme for the
random partition $\tilde{\pi}_{n}$. This is typically stated in terms of the
Chinese Restaurant Process (Pitman, 2006, Chapter 3), which is a sequential
construction of the random partition $\tilde{\pi}_{n}$ through the metaphor of
customers (observations) sitting at tables (species) of a restaurant. Under
Chinese Restaurant Process, the first customer $X_{1}$ arrives and is assigned
to a table. After $n$ customers $(X_{1},\ldots,X_{n})$ have arrived and have
been assigned to $k$ tables $\\{X_{1}^{\ast},\ldots,X_{k}^{\ast}\\}$, with
$n_{i}$ being the number of customers at table $i=1,\ldots,k$, the customer
$X_{n+1}$ arrives and

  * i)

she will sit at a (“new”) table $X^{\ast}$, that is a table not already
occupied, with a probability

$p^{\tiny{(new)}}=\frac{\theta+k\alpha}{\theta+n};$ (4)

  * ii)

she will sit at a table $X^{\ast}_{j}$ that has been already occupied, for
$j=1,\ldots,k$, with a probability

$p^{\tiny{(old)}}_{j}=\frac{n_{j}-\alpha}{\theta+n}.$ (5)

We refer to Pitman (2006, Chapter 3 and Chapter 4) for a detailed account on
Chinese Restaurant Process and its generalizations to SSMs. In particular, the
PYP is the sole SSM for which $p^{\tiny{(new)}}$ depends only on $(n,k)$ and
$p^{\tiny{(old)}}_{j}$ depends only on $(n,n_{j})$ (Bacallado et al., 2017).

Equation (3) is a symmetric function of compositions $(n_{1},\ldots,n_{k})$ of
$[n]$, that is the random partition $\tilde{\pi}_{n}$ induced by the PYP is an
exchangeable random partition (Pitman, 2006, Chapter 2). The ordered PYP is a
discrete random probability measure generalizing the PYP, in the sense that
random sampling from the ordered PYP allows to couple each species’ label with
a corresponding order (Gnedin and Pitman, 2005; James, 2006). If we denote by
$Q$ the ordered PYP with parameter $\alpha\in[0,1)$ and $\theta>0$, then a
random sample from $Q$ may be represented as
$((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$, with the $T_{i}$’s being the weights
inducing an ordering among species’ labels identifying the $X_{i}$’s. Because
of the (almost sure) discreteness of $Q$, the random sample
$((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ from $Q$ induces a random partition of
$[n]$ into $K_{n}\leq n$ blocks, labelled by a $K_{n}$-tuple
$((T^{\ast}_{1},X_{1}^{\ast}),\ldots,(T^{\ast}_{K_{n}},X_{K_{n}}^{\ast}))$
that is ordered according to the $T^{\ast}_{j}$’s in such a way that
$T^{\ast}_{1}>\cdots>T^{\ast}_{K_{n}}$, with corresponding ordered frequencies
$M_{j,n}=|{i\in[n]\text{ : }(T_{i},X_{i})=(T^{\ast}_{j},X^{\ast}_{j})}|$ for
$j=1,\ldots,K_{n}$ and such that $M_{j,n}\geq 1$ and $\sum_{1\leq j\leq
K_{n}}M_{j,n}=n$. Species’ labels $X_{j}^{\ast}$’s are thus ordered with
respect to the decreasing ordering of the their weights $T^{\ast}_{j}$’s,
namely the larger the weight the smaller the order, such that $M_{j,n}$ is the
frequency of the species’ label of order $j$ that corresponds to the $j$-th
largest weight $T^{\ast}_{j}$. An analogous construction follows for ordered
SSMs (James, 2006).

The predictive distribution of the ordered PYP $Q$ may be stated as an ordered
version of the Chinese Restaurant Process, with tables ordered according to
weights (James, 2006). Under the ordered Chinese Restaurant Process, the first
$n$ customers $((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ arrive and they are
assigned to the $k$ ordered tables
$((T^{\ast}_{1},X_{1}^{\ast}),\ldots,(T^{\ast}_{k},X_{k}^{\ast}))$, with the
table of order $j$ corresponding to the $j$-th largest weight $T^{\ast}_{j}$.
Hereinafter, we denote by $m_{j}$ the number of customers seated at the table
of order $j$ for $j=1,\ldots,k$, and we set $r_{j}=m_{j}+\cdots+m_{k}$, for
any $j=1,\ldots,k$, and $r_{k+1}=0$. Then, the customer $(T_{n+1},X_{n+1})$
arrives and

  * i)

she will sit at a “new” table $(T^{\ast},X^{\ast})$ of order $j=1,\ldots,k+1$,
that is a table not already occupied and whose order $j$ is determined through
the weight $T^{\ast}$, with a probability

$q^{\tiny{(new)}}_{j}=\frac{\theta+\alpha
r_{j}}{(1+r_{j})(\theta+n)}\prod_{i=1}^{j-1}\frac{r_{i}(\alpha
r_{i+1}+\alpha+\theta m_{i})}{(r_{i}+1)(\alpha r_{i+1}+\theta m_{i})};$ (6)

  * ii)

she will sit at a table $(T^{\ast}_{j},X^{\ast}_{j})$ that has been already
occupied, for $j=1,\ldots,k$, with a probability

$q^{\tiny{(old)}}_{j}=\frac{r_{j}(m_{j}-\alpha)(\alpha r_{j+1}+\theta
m_{j}+\theta)}{(1+r_{j})(\theta+n)(\alpha r_{j+1}+\theta
m_{j})}\prod_{i=1}^{j-1}\frac{r_{i}(\alpha r_{i+1}+\alpha+\theta
m_{i})}{(r_{i}+1)(\alpha r_{i+1}+\theta m_{i})}.$ (7)

We refer to Gnedin and Pitman (2005) and James (2006) for a detailed account
on (6) and (7). In particular, the predictive distribution of the ordered DP
arises from (6) and (7) by setting $\alpha=0$.

Equation (6) and Equation (7) provide a generative scheme for the random
partition induced by the ordered PYP. In particular, according to such a
generative scheme: i) if the $(n+1)$-th customer $(T_{n+1},X_{n+1})$ sits at
the new table $(T^{\ast},X^{\ast})$, which happens with probability (6), then
the order of such a table with respect to the ordering of the already occupied
tables $((T^{\ast}_{1},X_{1}^{\ast}),\ldots,(T^{\ast}_{k},X_{k}^{\ast}))$ is
determined by $T^{\ast}$, thus possibly changing the ordering of occupied
tables by shifting the order of tables $(T^{\ast}_{j},X^{\ast}_{j})$’s with
weights smaller than $T^{*}$; ii) if the $(n+1)$-th customer
$(T_{n+1},X_{n+1})$ sits at a table $(T^{\ast}_{j},X^{\ast}_{j})$ that is
already occupied, which happens with probability (7), then the order of such a
table with respect to the ordering of the already occupied tables is
determined by $T^{\ast}_{j}$, thus not changing the ordering of occupied
tables. In other terms, a new customer siting at a new table may determine a
change in the ordering of the already occupied tables, whereas a new customer
siting at a table already occupied does not determine a change in the ordering
of the already occupied tables.

Gnedin and Pitman (2005) and James (2006) first investigated distributional
properties of the random partition induced by the ordered PYP, and introduced
the notion of ordered EPPF. Generalizing the definition of EPPF, the ordered
EPPF is defined as the probability of any ordered partition of $[n]$ with $k$
blocks of frequencies $(m_{1},\ldots,m_{k})$. Gnedin and Pitman (2005) showed
that the ordered PYP induces a random partition whose ordered EPPF is

$\Phi^{(n)}_{k}(m_{1},\ldots,m_{k})=\frac{\prod_{i=1}^{k}\frac{\theta
m_{i}+\alpha
r_{i+1}}{r_{i}}}{(\theta)_{(n)}}\prod_{i=1}^{k}(1-\alpha)_{(m_{i}-1)}.$ (8)

See also Gnedin (2010), and references therein, for a comprehensive account on
ordered EPPFs and generalizations thereof. Note that the EPPF (3) can be
recovered from the ordered EPPF (8) by summing over the set $S_{k}$ of all
possible permutations of the $k$ blocks, that is

$\Pi_{n}^{(k)}(m_{1},\ldots,m_{k})=\sum_{\pi\in
S_{k}}\Phi^{(n)}_{k}(m_{\pi(1)},\ldots,m_{\pi(k)}).$ (9)

See Section A.1 of the Supplementary Material for details on Equation (9). The
distribution of the age-ordered partition of Donelly and Tavaré (1986) arises
from (8) by setting $\alpha=0$, where species’ labels are ordered according to
weights $T_{i}$’s that are interpreted as the ages of alleles. Another special
case of (8) is obtained by setting $\alpha\in(0,1)$ and $\theta=0$ (Favaro and
James, 2016).

By a direct application of the ordered EPPF (8), one may compute the
probability $P_{n}(i;\alpha,\theta)$ that a species with frequency $i$ has
species’ label of order $1$, i.e. the species’ label corresponding to the
largest weight $T^{\ast}_{1}$. For $\alpha=0$, Donelly and Tavaré (1986)
computed such a probability, showing that it is independent of $\theta$ and
also an increasing (linear) function of $i$, i.e.

$P_{n}(i;0,\theta)=\frac{i}{n}.$ (10)

Within the population genetic setting of Donelly and Tavaré (1986), Equation
(10) shows that the most frequent allele is the oldest allele. In general, for
any $\alpha\in[0,1)$ and $\theta>0$ it holds

$P_{n}(i;\alpha,\theta)=\frac{\alpha
n+i(\theta-\alpha)}{n}\mathds{E}\left[\frac{1}{\theta+\alpha K_{n-i}}\right],$
(11)

where $K_{n-i}$ is the number of distinct species in $(n-i)$ random samples
for the ordered PYP $Q$, with the proviso $K_{0}=0$ (Pitman, 2006, Chapter 3).
See Section A.2 of the Supplementary Material for the proof of Equation (11).
It is easy to show that (11) reduces to $\eqref{oldest_dp}$ for $\alpha=0$.
Differently from (10), the probability (11) depends on $(\alpha,\theta)$ and,
most importantly, it is no more an increasing (linear) function of $i$. Figure
1 shows that the probability (11) may increase or decrease in $i$ according to
the value of $(\alpha,\theta)$; for instance, for $\alpha\in(0,1)$ and
$\theta=0$ the probability (11) is the product of a term decreasing in $i$ and
a term increasing in $i$.

Figure 1: The probability $P_{n}(i;\alpha,\theta)$, where $n=1000$, as a
function of the frequency $i$, for different values of $\alpha$ and $\theta$.
Each panel corresponds to a different $\theta$, from top left to bottom right,
we have: $\theta=1,10,100,500$.

##  3 BNP inference for the first $r$ order frequencies

As SSPs, SSPs with ordering assume $n\geq 1$ observable samples from a generic
population of individuals, with each individual taking a value in a (possibly
infinite) discrete space of symbols, and then consider $m\geq 1$ additional
unobservable sample from the same population. The critical difference between
SSPs and SSPs with ordering lies in the definition of the (discrete)
functional of interest: while in SSPs such a functional is “invariant” with
respect to species ordering and deals with the species’ composition of the
additional samples, in SSPs with ordering the functional is not “invariant”
with respect to species ordering and deals with the species’ composition of
both the additional samples and the enlarged sample. The estimation of the
first $r$ order frequencies is arguably the most natural example of SSPs, as
it is the unseen-species problem in the context of SSPs. Assuming $n$
observable samples to be modeled as a random sample
$((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ from the ordered PYP $Q$, i.e.

$\displaystyle(T_{i},X_{i})\,|\,Q$
$\displaystyle\quad\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\quad
Q\qquad i=1,\ldots,n,$ (12) $\displaystyle Q$
$\displaystyle\quad\sim\quad\text{PYP}(\alpha,\theta),$

with $\text{PYP}(\alpha,\theta)$ being the law of the ordered PYP with
parameter $(\alpha,\theta)$, we introduce a BNP approach to estimate the
frequencies of the first $r$ order species in an enlarged sample obtained by
collecting $m$ additional samples
$((T_{n+1},X_{n+1}),\ldots,(T_{n+m},X_{n+m}))$ from the same $Q$.

###  3.1 Posterior distributions for the first $r$ order frequencies

We start by introducing a marginal distribution of the random partition
induced by $Q\sim\quad\text{PYP}(\alpha,\theta)$. In particular, for any
$n\geq 1$ let $((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ be a random sample under
the BNP model (12), such that the sample features $K_{n}=k$ distinct species
with corresponding ordered frequencies $\mathbf{M}_{n}=\mathbf{m}$.
Hereinafter, for the sake of simplicity in notation, we denote by
$|\mathbf{m}|_{1:r}$ the sum of the first $r$ elements of $\mathbf{m}$, i.e.
$|\mathbf{m}|_{1:r}=\sum_{1\leq j\leq r}m_{j}$, with the proviso
$|\mathbf{m}|_{1:0}=0$. For any index $r\in\\{1,\ldots,n\\}$ such that
$r\leq|\mathbf{m}|_{1:r}\leq n-k+r$, if we set

$C_{r,n}(\alpha,\theta,\mathbf{m})=\prod_{j=1}^{r}\frac{[\alpha(n-|\mathbf{m}|_{1:j})+\theta
m_{j}]}{n-|\mathbf{m}|_{1:j-1}}(1-\alpha)_{(m_{j}-1)}$

then

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}\geq r]$ (13)
$\displaystyle\quad={n\choose
m_{1},\ldots,m_{r},n-|\mathbf{m}_{1:r}|}\frac{C_{r,n}(\alpha,\theta,\mathbf{m})}{(\theta+n-|\mathbf{m}|_{1:r})_{|\mathbf{m}|_{1:r}}}.$

See Section A.4 of the Supplementary Material for the proof of Equation (13).
Equation (13) generalizes Donelly and Tavaré (1986, Proposition 6.1), which is
recovered from (13) by letting $\alpha\to 0$. For $r=1$, Equation (13)
provides the distribution of first order frequency, i.e.

$\text{Pr}[M_{1,n}=m_{1}]={n\choose m_{1}}\frac{\alpha(n-m_{1})+\theta
m_{1}}{n(\theta+n-m_{1})_{(m_{1})}}(1-\alpha)_{(m_{1}-1)}.$ (14)

Within the population genetic setting of Donelly and Tavaré (1986), Equation
(14) with $\alpha=0$ provides the distribution of the frequency of the oldest
allele (Kelly, 1977; Watterson, 1977).

Now, we can state our main results on the posterior distribution of the first
$r$ order frequencies, and in particular for $r=1$. Let
$((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ be a random sample from the ordered PYP
$Q$, and let $((T_{n+1},X_{n+1}),\ldots,(T_{n+m},X_{n+m}))$ be an additional
random sample from the same ordered PYP $Q$. We denote by $K_{m}^{(n)}$ the
number of distinct species in the sample
$((T_{n+1},X_{n+1}),\ldots,(T_{n+m},X_{n+m}))$ that are not in
$((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$, i.e. $K_{m}^{(n)}=K_{n+m}-K_{n}$, and
we denote by $W_{i,n+m}$ the frequency of the specie’s label of order $i$ in
the enlarged sample $((T_{1},X_{1}),\ldots,(T_{n+m},X_{n+m}))$, for
$i=1,\ldots,K_{n+m}$. To determine the distribution of the ordered frequencies
$\mathbf{W}_{n+m}=(W_{1,n+m},\ldots,W_{K_{n+m},n+m})$, it is useful to set

$A_{r}=\\{\text{species' labels with order }1,\ldots,r\text{ are new}\\},$
(15)

i.e., the event that the species’ labels with higher weights have not been
recorded in the first sample, and

$B_{r}=\\{\text{species' labels with order }1,\ldots,r\text{ are old}\\},$
(16)

i.e. the event that the observations with higher weights have been recorded in
the initial sample.

######  Theorem 1.

Let $((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ be a random sample under the BNP
model (12), such that the sample features $K_{n}=k$ distinct species with
corresponding ordered frequencies $\mathbf{M}_{n}=\mathbf{m}$. Moreover, let
$((T_{n+1},X_{n+1}),\ldots,(T_{n+m},X_{n+m}))$ be an additional random sample
under the same BNP model (12) such that the enlarged sample
$((T_{1},X_{1}),\ldots,(T_{n+m},X_{n+m}))$ features $K_{m+n}$ distinct species
with corresponding ordered frequencies $\mathbf{W}_{n+m}$, and set
$K_{m}^{(n)}=K_{n+m}-K_{n}$. If $A_{r}$ and $B_{r}$ are the events defined in
(15) and (16), respectively, then it holds:

  * i)

for $r\in\\{1,\ldots,n+m\\}$ such that $r\leq|\mathbf{w}|_{1:r}\leq m$,

$\displaystyle\text{Pr}[A_{r},W_{1,n+m}=w_{1},\ldots,W_{r,n+m}=w_{r},K_{m}^{(n)}\geq
r\,|\,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$ (17) $\displaystyle\quad={m\choose
w_{1},\ldots,w_{r},m-|\mathbf{w}|_{1:r}}\frac{C_{r,n+m}(\alpha,\theta,\mathbf{w})}{(\theta+n+m-|\mathbf{w}|_{1:r})_{(|\mathbf{w}|_{1:r})}};$

  * ii)

for $r\in\\{1,\ldots,k\\}$ such that $0\leq|\mathbf{w}|_{1:r}\leq m$,

$\displaystyle\text{Pr}[B_{r},W_{1,n+m}=w_{1}+m_{1},\ldots,W_{r,n+m}=w_{r}+m_{r}\,|\,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
(18) $\displaystyle\quad\times{m\choose
w_{1},\ldots,w_{r},m-|\mathbf{w}|_{1:r}}\frac{\frac{C_{r,n+m}(\alpha,\theta,\mathbf{w}+\mathbf{m})}{(\theta+n+m-|\mathbf{w}+\mathbf{m}|_{1:r})_{(|\mathbf{w}+\mathbf{m}|_{1:r})}}}{\frac{C_{r,n}(\alpha,\theta,\mathbf{m})}{(\theta+n-|\mathbf{m}|_{1:r})_{(|\mathbf{m}|_{1:r})}}}$

See Section A.5 of the Supplementary Material for the proof of Theorem 1.
Theorem 1 may be viewed as the posterior counterpart of Equation (13), with
respect to an initial observable sample
$((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$. In particular, Equation (17) and
Equation (18) provide two posterior distributions of the first $r$ order
frequencies under the events $A_{r}$ and $B_{r}$, respectively, for $r\geq 1$.
Equation (17) provides the posterior distribution of the first $r$ order
frequencies having species’ labels not belonging to the additional observable
samples; that is the ordering of species’ labels in the initial sample
$((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ is changed according to the additional
sample $((T_{n+1},X_{n+1}),\ldots,(T_{n+m},X_{n+m}))$. Equation (18) provides
the posterior distribution of the first $r$ order frequencies having species’
labels belonging to the additional observable samples; that is the ordering of
species’ labels in the initial sample $((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$
is not changed according to the additional sample
$((T_{n+1},X_{n+1}),\ldots,(T_{n+m},X_{n+m}))$. BNP estimators of the first
$r$ order frequencies, with respect to a squared loss function, are obtained
in terms of the corresponding posterior expectations, i.e. the vectors of
expected values with respect to the posterior distributions (17) and (18). As
a corollary of Theorem 1, we obtain the posterior distributions of the
frequency of order $1$.

######  Corollary 1.

Let $((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ be a random sample under the BNP
model (12), such that the sample features $K_{n}=k$ distinct species with
corresponding ordered frequencies $\mathbf{M}_{n}=\mathbf{m}$. Moreover, let
$((T_{n+1},X_{n+1}),\ldots,(T_{n+m},X_{n+m}))$ be an additional random sample
under the same BNP model (12) such that the enlarged sample
$((T_{1},X_{1}),\ldots,(T_{n+m},X_{n+m}))$ features $K_{m+n}$ distinct species
with corresponding ordered frequencies $\mathbf{W}_{n+m}$, and set
$K_{m}^{(n)}=K_{n+m}-K_{n}$. If $A_{1}$ and $B_{1}$ are the events defined in
(15) and (16), respectively, then it holds:

  * i)

$\displaystyle\text{Pr}[A_{1},W_{1,n+m}=w_{1},K_{m}^{(n)}\geq
1\,|\,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$ (19)
$\displaystyle\quad=\binom{m}{w_{1}}\frac{\alpha(n+m-w_{1})+\theta
w_{1}}{(n+m)(\theta+n+m-w_{1})_{(w_{1})}}(1-\alpha)_{(w_{1}-1)};$

  * ii)

$\displaystyle\text{Pr}[B_{1},W_{1,n+m}=w_{1}+m_{1}\,|\,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
(20)
$\displaystyle\quad=\binom{m}{w_{1}}\frac{\frac{\alpha(n+m-w_{1}-m_{1})+\theta(w_{1}+m_{1})}{(n+m)(\theta+n+m-w_{1}-m_{1})_{w_{1}+m_{1}}}}{\frac{\alpha(n-m_{1})+\theta
m_{1}}{n(\theta+n-m_{1})_{(m_{1})}}}(m_{1}-\alpha)_{(w_{1})};$

  * iii)

$\displaystyle\text{Pr}[W_{1,n+m}=w\,|\,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
(21) $\displaystyle\quad=\frac{[\alpha(n+m-w)+\theta
w](1-\alpha)_{(w-1)}}{(n+m)(\theta+n+m-w)_{(w)}}$
$\displaystyle\quad\quad\times\Big{[}\mathds{1}_{\\{1,\ldots,m\\}}(w)\binom{m}{w}+\mathds{1}_{\\{m_{1},\ldots,m_{1}+m\\}}(w)\binom{m}{w-m_{1}}\frac{n(\theta+n-m_{1})_{(m_{1})}}{[\alpha(n-m_{1})+\theta
m_{1}](1-\alpha)_{(m_{1}-1)}}\Big{]}.$

The proofs of Equation (19) and Equation (20) follow by a direct application
of Theorem 1, that is by setting $r=1$ in Theorem 1, whereas Equation (21)
follows by combining (19) and (20). Within the population genetic setting of
Donelly and Tavaré (1986), Equation (21) with $\alpha=0$ provides the
posterior distribution of the frequency of the oldest alleles. By exploiting
Corollary 1, BNP estimators of the frequency of order $1$, with respect to a
squared loss function, are obtained in terms of the expected values of the
posterior distributions (19), (20) and (21). Here, we report the BNP estimator
with respect to the posterior distribution (21) and we refer to Section A.6 of
the Supplementary Material for the BNP estimators with respect to the
posterior distributions (19) and (20), In particular, if we set

$\displaystyle C(\alpha,\theta,n,m,m_{1})$
$\displaystyle\quad=[\alpha(n+m-m_{1})+\theta
m_{1}]\left[m_{1}+m\frac{m_{1}-\alpha}{\theta+n-\alpha}\right]$
$\displaystyle\quad\quad+\left[m(\theta-\alpha)\frac{m_{1}-\alpha}{\theta+n-\alpha}\right]\left[(m_{1}+1)+(m-1)\frac{m_{1}+1-\alpha}{\theta+n+1-\alpha}\right]$

then

$\displaystyle\mathds{E}[W_{1,n+m}\,|\,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
(22)
$\displaystyle\quad=\frac{m(\theta+n+1-\alpha)_{(m)}}{(n+m)(\theta+n)_{(m)}}\frac{\theta+n\alpha}{\theta+n+1-\alpha}+\frac{n(\theta+n-\alpha)_{(m)}}{(n+m)(\theta+n)_{(m)}}\frac{C(\alpha,\theta,n,m,m_{1})}{\alpha(n-m_{1})+\theta
m_{1}}.$

We refer to Section A.7 of the Supplementary Material for explicit expressions
of the posterior probabilities of the events $A_{1}$ and $B_{1}$, which
complete the main result of Corollary 1.

###  3.2 Estimation of prior’s parameter $(\alpha,\theta)$

To perform the posterior inference described above using the ordered PYP
prior, the prior’s parameter $(\alpha,\theta)$ need to be specified in a
suitable way. Here we consider an empirical Bayes approach, in which the
parameters are estimated given an observed sample of $n$ observations. In
particular, we consider two classes of general methods: maximum likelihood
estimation methods, and methods of moments combined with least square. For the
first class of methods, we find the optimal values of $\alpha$ and $\theta$
that maximize the ordered EPPF given the species’ frequencies. EPPFs from
different models can be used; we consider the EPPF corresponding to the PYP,
the ordered PYP and the ordered DP The sets of prior’s parameters obtained by
optimizing these EPPF are respectively denoted as stdPYP, ordPYP and ordDP.
For the second class of methods, we consider a statistic of interest and match
their population first moment with the observed sample statistic, for a grid
of sample size values. Specifically, we consider $1\leq n_{1}\leq\ldots\leq
n_{d}=n$ equally spaced; for each $n_{i}$, we compute the discrepancy between
the first moment and the observed statistics for the first $n_{i}$
observations; finally, we minimize the sum of the squared discrepancies, where
the sum is for $i=1,\ldots,d$. We consider the prior’s parameters that
minimize the least square for the frequency of the first ordered species
$M_{1,n}$ (lsX1) and the ones that minimize the least squares for the number
of distinct species $K_{n}$ (denoted with lsK). For example, for the number of
unique species $K_{n}$, the objective function to be minimized is
$\sum_{i=1}^{d}\left(\mathds{E}[K_{n_{i}}]-K_{n_{i}}\right)^{2}$.

##  4 Numerical illustrations

We empirically study the ordered PYP and assess the performance of our BNP
approach for estimating the frequency of the first ordered species, using
synthetic and real data. Moreover, we empirically study distributional
properties of the ordered random partitions induced by the model, in
particular focusing on the species’ ordering distribution. In the synthetic
data, we compare the performance when the data is generated from the model,
i.e. the PYP prior is correctly specified, and when the data is generated from
different distributions, thus under model misspecification. For the
application of our model to real data, we analyze genetic variation using
samples from the 1000 Genome Project (1000 Genomes Project Consortium, 2015),
which we combine with variants’ age estimates obtained from the Human Genome
Dating Project (Albers and McVean, 2020). Studying genetic variation is of
great importance to investigate population structure, to detect related
samples, to investigate demographic history, and to learn about the risk of
diseases and different quantitative traits. By also incorporating information
on the variants’ age, similar and further issues can be assessed, such as
using the variants’ age distribution to compare populations, differentiating
age distributions in pathogenic and benign variants, and learning about
genealogical history (Albers and McVean, 2020). Here we assess the performance
of our method on predicting frequencies of variants ordered by their age,
focusing in particular on the oldest variant. Code is available at
https://github.com/cecilia-balocchi/OrderedSSP.

###  4.1 Preliminaries

While the distributional properties of the random partition induced by the PYP
are well-known in the BNP literature, the properties induced by the ordered
PYP are less understood. Because of the marginality property (9), some
distributional properties of the ordered PYP are equivalent to the ones of the
PYP. For example, the distribution of the number of distinct species (or
clusters) induced by the ordered PYP is equal to the one induced by the PYP,
reported in equation (25). However, other properties that relate to the
ordering on the species are not as well understood. In particular, we are
interested in learning the behavior of the distribution that assigns the order
to a new cluster. We aim to characterize it using simple descriptive features.
We achieve this goal by studying the predictive distribution (6) that assigns
the $n+1$ observation to a new species of order $j$, and to examine its
behavior marginalizing on all the configurations of partitions of $n$
observations into $K_{n}$ species. In other words, we empirically study the
marginal distribution that a new species after $n$ observations is assigned
order $j$, given that the previous $n$ observations are partitioned into
$K_{n}$ species, for $j=1,\ldots,K_{n}+1$. This distribution is “marginal”
compared to (6), because it does not condition on the frequencies of the
ordered species.

In Figure 2 we represent the marginal ordering distribution for a new species
given an observed sample of size $n=10$, for different configurations
(parameters) of the PYP prior. In particular, from left to right, we represent
the distribution under the ordered DP ($\theta>0,\alpha=0$, first panel), the
ordered PYP with $\alpha<\theta$, $\alpha=\theta$, and $\alpha>\theta$ and the
ordered $\alpha$-stable process ($\theta=0,\alpha>0$, last panel). Figure 2
shows that the ordering distribution changes depending on the parameters
$\theta$ and $\alpha$: for $\theta>\alpha$ the probability that a new cluster
is assigned to order $j$ decreases with $j$ for each previous number of
species $K_{n}$ (first and second panels from the left) and the decreasing
trend in stronger when the difference $\theta-\alpha$ is large. For
$\theta=\alpha$ (third panel from the left) the trend is constant over $j$,
for all $K_{n}$. For $\theta<\alpha$ we see instead that the trend is
increasing with $j$, for each $K_{n}$ (fourth panel). The last panel shows
that in the case of the $\alpha$-stable process ($\theta=0<\alpha$) the trend
is increasing from $j=2$ and constant afterwards. These intuitions are useful
for constructing a distribution for ordered partitions that has properties
similar to those induced by the ordered PYP prior.

Figure 2: Order distribution for a new species, given a sample of $n=10$
observations, divided into $K_{n}=1,\ldots,10$ species (each color represents
a different value of $K_{n}$). Colored curves represent the marginal
distribution and the points show the variation across different partition
configurations. The dashed black line represents the order distribution for a
new species, marginally on the number of previous species $K_{n}$.

###  4.2 Analysis of synthetic data

Assuming synthetic data, we consider the problem of making inference on some
quantities of interest with our BNP approach. In particular, we study the
performance on an additional sample of size $m$ of the posterior mean
predictors for the total number of species $K_{n+m}$, the frequency of the
oldest cluster $W_{1,n+m}$, and the frequency $W_{1,n+m}$ conditionally on the
knowledge that either the event $A_{1}$ or $B_{1}$ happened ($W_{1,n+m}|A_{1}$
and $W_{1,n+m}|B_{1}$). For notational simplicity we will remove the
dependence on $n+m$. We compare the performance of the ordered PYP model, with
prior’s parameters estimated with the methods described in Section 3.2. We
generate the synthetic data under different scenarios. We first consider a
framework where the model is correctly specified, i.e. the data is generated
from the ordered PYP prior. We then focus on a framework where the model is
misspecified, as the data is generated from different distributions. We
evaluate the performance of our estimators derived from the ordered PYP model,
comparing the behavior when the prior’s parameters are estimated with the
methods described in Section 3.2. In Section B.1 of the Supplementary Material
we report numeric results related to the evaluation described in the following
sections.

####  4.2.1 Inference of the first ordered frequency under correct
specification

We assume the data from the ordered PYP prior, i.e. under correct model
specification. We compute the absolute percentage error for our four
quantities of interest, across 100 simulations. Specifically, we consider 100
datasets of size $n=500$, generated from the ordered PYP prior with randomly
sampled parameters $(\theta,\alpha)$; for each dataset, we consider 25
additional datasets of size $m=5000$, and compute the median prediction error
across the 25 additional datasets.

The first panel from the left of Figure 3 shows that estimating the parameters
using a least squares method for the number of distinct species $K$ (lsK)
achieves the best predictive performance, followed by estimating the
parameters maximizing the EPPF of the standard PYP (stdPYP) and of the ordered
PYP (ordPYP). For these methods the median percentage error is less than 0.05
($5\%$), suggesting quite an accurate prediction. The second panel shows the
absolute percentage error for the frequency of the oldest species $W_{1}$,
with the best predictive performance achieved by lsX1, with $65\%$ median
error. Note that the percentage errors are larger, and we have worse
predictions compared to the number of distinct species $K$. This is most
likely due to the added complexity of predicting if event $A_{1}$ (the first
order species is a new species) or $B_{1}$ (the first order species was
observed in the initial sample) will occur. In the former case, the frequency
$W_{1}$ will most likely be small, while in the latter case $W_{1}$ will take
on very large values. This discrepancy causes the error to be much larger. For
example, ordPYP achieves a median error of $123\%$. When predicting the
conditional frequencies $W_{1}|A_{1}$ and $W_{1}|B_{1}$, the errors are
smaller, and the parameter-estimating method ordPYP achieves the best
performance with median percentage errors respectively of $50\%$ and $48\%$.

Figure 3: Predictive performance for our BNP approach for the total number of
distinct species $K$, the frequency of the oldest cluster $W_{1}$, and the
conditional frequencies $W_{1}|A_{1}$ and $W_{1}|B_{1}$ (shown in different
columns). The boxplots display the median absolute percentage error between
the predicted posterior mean and the true value, across several datasets
simulated from the ordered PYP prior. Various parameter-estimating methods are
compared.

####  4.2.2 Inference of the first ordered frequency under misspecification

We assume the data not from the ordered PYP prior, i.e. under model
misspecification. We aim at evaluating how our method performs in such adverse
conditions, which are in general also to be expected from real data. We
consider several different data-generating processes. In all of them, we
generate the species order independently from the observations’ species (or
cluster) assignment. For the clustering distribution we consider: (a) the
standard DP, (b) the standard PYP, and (c) the (infinite support) Zipf
distribution (also known as the zeta distribution), where each observation is
associated to a integer (sampled from the Zipf distribution) and clusters are
formed by aggregating observations mapped to the same integer. Given that
these distributions do not induce an order on the species, we need to sample
an order for each new species. For the ordering distribution we consider: (1)
the alpha-stable distribution (induced by the ordered $\alpha$-stable process,
a special case of the ordered PYP with $\theta=0$ where the predictive
distribution (6) simplifies significantly), and (2) what we call the arrival-
weighted distribution. The latter considers the ordering induced by the
cluster arrival and introduces a random component by sampling for each new
cluster an exponentially distributed weight with mean given by the arrival
order, and ordering the species according to the weight. We generate 100
datasets of size $n=500$ using randomly sampled parameters, and for each of
them consider 25 additional datasets of size $m=5000$; we consider the median
absolute percentage errors across the 25 additional datasets.

Figure 4 displays the predictive performance results for the number of species
$K$ in the first row, comparing the performance of the different parameter-
estimation methods across the different data-generating distributions.
Similarly to the results presented in Section 4.2.1, the predictions using the
stdPYP and lsK parameter-estimation methods achieve the best performance, with
the former performing better for data generated from a DP and PYP clustering
distribution, and the latter being best for data generated from a Zipf
distribution. This is not surprising, given that the Zipf distribution has
behavior similar to the standard PYP, but might require a more “robust”
methods of estimation. The second row of figure 4 shows the error committed
when predicting the frequency $W_{1}$. The range of percentage error committed
is again larger than the one for $K$.Under the alpha-stable ordering, the best
performance are achieved by ordPYP and lsX1 (with median error ranging from
$61\%$ when the DP clustering-distribution, and $126\%$ under the PYP
clustering-distribution). The predictions worsen when data is generated from
an arrival-weighted distribution, with the best performances achieved by the
ordPYP and lsX1 when data are generated from the DP (median errors of
$400\%$), and by lsK and stdPYP when data are generated from the PYP or Zipf
(best median errors of $880\%$ and $260\%$). Finally, the third and the fourth
rows show the predictive performance for the conditional frequencies
$W_{1}|A_{1}$ and $W_{1}|B_{1}$, and again the predictions for this quantities
are much more precise, with smallest median errors ranging between $30\%$ and
$160\%$. The best performances for $W_{1}|A_{1}$ are achieved mostly by lsX1
and ordPYP, while for $W_{1}|B_{1}$ by stdPYP, ordDP and lsK.

Figure 4: Predictive performance for our BNP approach, for the total number of
clusters $K$, the frequency of the oldest cluster $W_{1}$, and the conditional
frequencies $W_{1}|A_{1}$ and $W_{1}|B_{1}$ (shown in different rows). The
boxplots display the median absolute percentage error between the predicted
posterior mean and the true value, across several datasets simulated under
different distributions (shown in different columns). We compare various
parameter-estimating methods.

Overall, when the model is misspecified, lsX1 seems to be more robust for the
estimation of both the marginal distribution and the conditional distribution
of $W_{1}$, but it performs quite poorly for the prediction of $K$. While the
parameter-estimation method ordPYP does not always produce the best
predictions, it seems to provide a good balance between learning the number of
species $K$ and the frequency of the first order species $W_{1}$.
Alternatively, separate parameter-estimation methods could be used for
predictions of the two quantities.

###  4.3 Analysis of genetic variation data

We employ our method to analyze genetic variation data from the 1000 Genomes
Project. We examine single nucleotide polymorphisms (SNPs) corresponding to
certain genes for a sample of 2548 individuals. In particular, we consider a
variant to be present for a given individual and SNP locus if the present DNA
base is different from the reference genome in either allele. We combine the
variants from all individuals and all SNPs corresponding to a certain gene, to
create our basic sample, where the variant location (i.e. the SNP) represents
the species each sample unit belongs to. To associate a ordering to each
species, we used SNPs age estimates from Human Genome Dating project (Albers
and McVean, 2020), and discard any variant for which no age information is
available. We study two datasets, formed by the genetic variants corresponding
respectively to the BRCA2 and the EDAR genes, for which we analyze
respectively 1482 and 1073 unique SNPs. We focus on the predictive performance
for the number of distinct variants $K$ and the frequency of the oldest
variant $W_{1}$, using the different parameter-estimation methods for our BNP
approach. We repeat our analysis for 100 different training and testing sets,
randomly sampled so that the testing set size is approximately 20 times larger
than the training set size.

Figure 5: Prediction performance across several training-testing sets for
BRCA2 (panel one and three) and EDAR (panels two and four).

In Figure 5 we report the percentage errors for predicting $K$ and $W_{1}$ for
the datasets corresponding to the BRCA2 and the EDAR genes. For the prediction
of the number of species, we find the methods ordPYP, stdPYP and lsK provide
the best results. In terms of predicting the frequency of the oldest variant,
ordDP and lsX1 perform better. We note that the performance of ordPYP is
comparable to that of ordDP in the BRCA2 dataset, as the former often
estimated values of $\alpha$ close to zero; however, the two methods achieved
different performance in the EDAR dataset, due to the fact that ordPYP
estimated positive values of $\alpha$. Under further inspection, we note that
lsX1 has performance comparable to ordDP in the prediction of $W_{1}$ for the
EDAR dataset because it often estimates values of $\alpha$ very close to zero.
It is also surprising to see that these datasets seem to have power-law
behavior in terms of the number of species, for which values of $\alpha$
greater than 0 provide more accurate predictions, but the best predictions for
the frequency of the oldest variant $W_{1}$ are produced by methods that
estimate $\alpha=0$. This is probably due to the model being not correctly
specified for these data. Thus we recommend analyzing the two predictions
problems separately.

Figure 6: Prediction for the curve for BRCA2, for $K$ (top panels) and $W_{1}$
(bottom panels).

Our method can also be used to analyze the predictions as functions of the
number of samples, rather than simply make a prediction for the entire test
set. Using the multiple training and testing splits method, we also looked at
the performance in the prediction of the curve of $K$ and $W_{1}$. Figure 6
reports the results for each parameter-estimating method, applied to the BRCA2
dataset. Similar plots are reported in Section B.2 of the Supplementary
Material for the EDAR gene dataset. For simplicity of visualization, we depict
in black the actual curve of $K$ (or respectively $W_{1}$) observed in the
training-testing split that most closely represents the median error achieved
by each method. The prediction curve (depicted in red) is also the one
corresponding to the “median” training-testing split, while the red bands
represent the $95\%$ confidence bands. We report the empirical confidence
bands, computed using the empirical quantiles of the curve estimates across
the various training-testing splits. For the estimation of the number of
distinct species $K$, the top panel of Figure 6 confirm that stdPYP and lsK
have the best performance, with quite accurate prediction and with the
represented curve falling within the 95% bands. In the other parameter-
estimating methods the prediction is much worse, and often the curve does not
fall within the confidence bands. In terms of predictions for $W_{1}$, the
results are consistent with the expectations from the predictions on the whole
test set, with ordDP, ordPYP and lsX1 performing better, with very accurate
prediction of $W_{1}$ as a function of the number of additional samples.

##  5 Discussion

A common feature of SSPs is the invariance with respect to species labelling,
i.e. species’ labels are immaterial in defining the functional of interest,
which is at the core of the development of the BNP approach to SSPs under the
popular PYP prior (Lijoi et al., 2007, 2008; Favaro et al., 2009, 2013). In
this paper, we considered SSPs that are not “invariant” to species labelling,
in the sense that an ordering or ranking is assigned to species’ labels, and
we developed a BNP approach to such problems. Inspired by the seminal work of
Donelly and Tavaré (1986) on age-ordered alleles’ compositions, with a
renowned interest in the frequency of the oldest allele (Crow, 1972; Kelly,
1977; Watterson, 1977), we studied the following SSP with ordering: given an
observable sample from unknown population of individuals belonging to some
species (alleles), with species’ labels being ordered according to some
weights (ages), estimate the frequencies of the first $r$ order species’
labels in an enlarged sample obtained by including additional unobservable
samples from the same population. Our BNP approach relies on the ordered PYP
prior, which leads to an explicit posterior distribution of the first $r$
order frequencies, with corresponding estimates being simple and
computationally efficient. We presented an empirical validation of our
approach on both synthetic and real data. In particular, we applied our
approach to the analysis of the genetic variation, showing its effectiveness
in the estimation of the frequency of the oldest allele.

The sampling structure of the ordered PYP has been first presented in Gnedin
and Pitman (2005) and James (2006), and since then no other works have further
investigated such a sampling structure in BNPs. That is, differently from the
sampling formulae arising from the PYP prior, to date the sampling formulae
arising from the ordered PYP prior appear to be unknown in the BNP literature.
Our work highlights the great potential of the ordered PYP in BNP inference
for SSPs with ordering, paving the way to some directions for future research.
First, in our work we considered the problem of estimating the frequency of
the first $r$ order species, which is arguably the most natural SSP with
ordering; other (discrete) functionals of the ordered species’ composition of
unobservable samples may be of interest, e.g. the number of unseen species
with order less than $r$ that would be observed in additional samples, and the
number of unseen species observed with a certain prevalence and order less
than $r$ in the enlarged sample. Second, while SSPs with ordering have a
natural motivation in population genetics, they may be of interest in
different fields; in Section B.3 of the Supplementary Material we present an
application in the context of citations to academic articles, where each cited
article is a species whose order in determined by the publication date;
another application is in the context of online purchasing of items, with
order being determined by a suitable measure of popularity assigned to items.

Interest in ordered SSMs is not limited to BNP inference for SSPs with
ordering. Ordered SSMs, and in particular the ordered PYP prior, may be also
applied to a setting in which species are not explicitly observed, but need to
be inferred. This is the case of Bayesian mixture models (Frühwirth-Schnatter
et al., 2019), where each observation is assumed to be assigned to a latent
component (species’ label), which characterizes some features of the
distribution of the observations. In such a setting, ordered SSMs may be used
to specify the distribution of the mixture’s probabilities of the latent
components, with such a components being ordered according to some weights
(James, 2006a). Ordered SSMs admit a natural extension to the features
sampling framework, which generalizes the species sampling framework by
allowing each observation to belong to multiple species, now called features.
Feature sampling problems first appeared in ecology for modeling the presence
or absence of an animal in a trap, and their importance has grown dramatically
in recent years driven by applications in biological and physical sciences
(Masoero et al., 2021, 2022). In such a context, the Beta process prior
(Broderick et al., 2013) is the most popular nonparametric prior for modeling
the unknown feature’s composition of the population. The definition of an
ordered version of the Beta process prior, and generalizations thereof, would
be the starting point to introduce and investigate feature sampling problems
with ordering.

## Acknowledgments

The authors thank Paul Jenkins for useful discussions on the use of age-
ordered random partitions in population genetics. Cecilia Balocchi, Federico
Camerlenghi and Stefano Favaro received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme under grant agreement No 817257. Federico Camerlenghi and Stefano
Favaro gratefully acknowledge the support from the Italian Ministry of
Education, University and Research (MIUR), “Dipartimenti di Eccellenza” grant
2018-2022. Federico Camerlenghi is a member of the Gruppo Nazionale per
l’Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the
Istituto Nazionale di Alta Matematica (INdAM).

## References

  * 1000 Genomes Project Consortium (2015) 1000 Genomes Project Consortium. (2015). A global reference for human genetic variation. Nature 526, 68.
  * Albers and McVean (2020) Albers, P. K. and McVean, G. (2020). Dating genomic variants and shared ancestry in population-scale sequencing data. PLoS biology 18, e3000586.
  * Bacallado et al. (2017) Bacallado, S., Battiston, M., Favaro, S. and Trippa, L. (2015). Sufficientness postulates for Gibbs-type priors and hierarchical generalizations. Statistical Science 32, 487–500.
  * Balocchi et al. (2022) Balocchi, C., Favaro, S. and Naulet, Z. (2022). Bayesian nonparametric inference for “species-sampling” problems. Preprint available at: http://arxiv.org/abs/2203.06076.
  * Bunge and Fitzpatrick (1993) Bunge, J. and Fitzpatrick, M. (1993) Estimating the number of species: a review. Journal of the American Statistical Association 88, 364–373.
  * Broderick et al. (2013) Broderick, T., Jordan, M.I. and Pitman, J. (2013). Cluster and feature modeling from combinatorial stochastic processes. Statistical Science 28, 289-–312.
  * Charalambides (2005) Charalambides, C. (2005) Combinatorial methods in discrete distributions. Wiley.
  * Crow (1972) Crow, J.F. (1972) The dilemma of nearly neutral mutations; how important are they for the evolution and human welfare. Journal of Heredity 63, 306–316.
  * Deng et al. (2019) Deng, C. Daley, T., De Sena Brandine, G. and Smith, A.D. (2019). Molecular heterogeneity in large-scale biological data: techniques and applications. Annual Review of Biomedical Data Science 2, 39–67
  * Donelly and Tavaré (1986) Donelly, P. and Tavaré, S. (1986) The age of alleles and a coalescent. Advances in Applied Probability 18, 1–19
  * Efron and Thisted (1976) Efron, B. and Thisted, R. (1976). Estimating the number of unseen species: How many words did Shakespeare know? Biometrika 63, 435–447.
  * Ewens (1972) Ewens, W.J. (1972) The sampling theory of selectively neutral alleles. Theoretical population biology 3, 87–112.
  * Favaro and James (2016) Favaro, S. and James, L.F. (2016). Relatives of the Ewens sampling formula in Bayesian nonparametrics. Statistical Science 31, 30–33.
  * Favaro et al. (2013) Favaro, S., Lijoi, A. and Prünster, I. (2013). Conditional formulae for Gibbs-type exchangeable random partitions. The Annals of Applied Probability 23, 1721–1754.
  * Favaro et al. (2009) Favaro, S., Lijoi, A., Mena, R.H. and Prünster, I. (2009). Bayesian nonparametric inference for species variety with a two parameter Poisson-Dirichlet process prior. Journal of the Royal Statistical Society Series B 71, 992–1008.
  * Ferguson (1973) Ferguson, T.S. (1973). A Bayesian analysis of some nonparametric problems. The Annals of Statistics 1, 209–230.
  * Frühwirth-Schnatter et al. (2019) Frühwirth-Schnatter, S., Celeux, G. and Robert, C.P. (2019) Handbook of mixture analysis. Chapman and Hall/CRC.
  * Gehrke et al. (2003) Gehrke, J., Ginsparg, P. and Kleinberg, J. (2003). Overview of the 2003 KDD Cup. textitAcm Sigkdd Explorations Newsletter, 5, 149–151.
  * Gnedin (2010) Gnedin, A (2010). Regeneration in random combinatorial structures. Probability Surveys 7, 105–156.
  * Gnedin and Pitman (2005) Gnedin, A. and Pitman, J. (2005). Regenerative composition structures. The Annals of Probability 33, 445-479.
  * Good and Toulmin (1956) Good, I.J. and Toulmin, G.H. (1956). The number of new species, and the increase in population coverage, when a sample is increased. Biometrika 43, 45–63.
  * Griffiths (1979) Griffiths, R.C. (1979). Exact sampling distributions from the infinite neutral alleles model. Advances in Applied Probability 11, 326–354.
  * Hajra and Sen (2005) Hajra, K. B. and Sen, P. (2005). Aging in citation networks. Physica A: Statistical Mechanics and its Applications 346, 44-48.
  * James (2006) James, L.F. (2006). Poisson calculus for spatial neutral to the right processes. The Annals of Statistics 34, 416–440.
  * James (2006a) James, L.F. (2006). Spatial neutral to the right species sampling mixture models. Preprint arXiv:math/0604266v2
  * Ji and Jin (2016) Ji, P. and Jin, J. (2016). Coauthorship and citation networks for statisticians. The Annals of Applied Statistics, 10, 1779–1812.
  * Kelly (1977) Kelly, F.P. (1977) Exact results for the Moran neutral allele model. Advances in Applied Probability 9, 197–201.
  * Kingman (1975) Kingman, J.F.C. (1975) Random discrete distributions. Journal of the Royal Statistical Society Series B 37, 1–22.
  * Kingman (1978) Kingman, J.F.C. (1978) The representation of partition structures. Journal of the London Mathematical Society 18, 374–380.
  * Lijoi et al. (2007) Lijoi, A., Mena, R.H. and Prünster, I. (2007). Bayesian nonparametric estimation of the probability of discovering new species. Biometrika 94, 769–786.
  * Lijoi et al. (2008) Lijoi, A., Prünster, I. and Walker, S.G. (2008). Bayesian nonparametric estimators derived from conditional Gibbs structures. The Annals of Applied Probability 18, 1519–1547.
  * Masoero et al. (2022) Masoero, L., Camerlenghi, F., Favaro, S. and Broderick, T. (2020). More for less: predicting and maximizing genetic variant discovery via Bayesian nonparametrics. Biometrika 109, 17–32.
  * Masoero et al. (2021) Masoero, L., Camerlenghi, F., Favaro, S. and Broderick, T. (2020). Scaled process priors for Bayesian nonparametric estimation of the unseen genetic variation. Preprint arXiv:2106.15480.
  * Mathieson and McVean (2014) Mathieson, I. and McVean, G. (2014). Demography and the age of rare variants. PLoS genetics, 10, e1004528.
  * Orlitsky et al. (2016) Orlitsky, A., Suresh, A.T. and Wu, Y. (2017). Optimal prediction of the number of unseen species. Proceeding of the National Academy of Sciences of USA 113, 13283–13288.
  * Perman et al. (1992) Perman, M., Pitman, J. and Yor, M. (1992). Size-biased sampling of Poisson point processes and excursions. Probability Theory and Related Fields 92, 21–39.
  * Pitman (1995) Pitman, J. (1995). Exchangeable and partially exchangeable random partitions. Probability Theory and Related Fields 102, 145–158.
  * Pitman (1996) Pitman, J. (1996). Some developments of the Blackwell-MacQueen urn scheme. In Statistics, Probability and Game Theory, Ferguson, T.S., Shapley, L.S. and MacQueen, J.B. Eds., Institute of Mathematical Statistics.
  * Pitman (2006) Pitman, J. (2006). Combinatorial stochastic processes. Lecture Notes in Mathematics, Springer Verlag.
  * Pitman and Yor (1997) Pitman, J. and Yor, M. (1997). The two parameter Poisson-Dirichlet distribution derived from a stable subordinator. The Annals of Probability 25, 855–900.
  * Portenoy et al. (2017) Portenoy, J., Hullman, J. and West, J. D. (2017). Leveraging citation networks to visualize scholarly influence over time. Frontiers in Research Metrics and Analytics, 2, 8.
  * Watterson (1977) Watterson, G.A (1977). Is the most frequent allele the oldest? Theoretical Population Biology 11, 141–160.

##  Appendix A Proofs and details

###  A.1 Proof of (9)

The proof of (9) is based on the following identity

$\sum_{\pi\in
S_{k}}\theta^{-1}\prod_{j=1}^{k}\frac{\alpha(m_{\pi(j+1)}+\ldots+m_{\pi(k)})+\theta
m_{\pi(j)}}{(m_{\pi(j)}+\ldots+m_{\pi(k)})}=\prod_{i=1}^{k-1}(\theta+i\alpha).$
(23)

that we are going to show by induction. When $k=1$ this is trivially true. Now
assume that (inductive hypothesis) Equation (23) is true with $k-1$ in place
of $k$, thus one has

$\displaystyle\frac{1}{\theta}\sum_{\pi\in
S_{k}}\prod_{j=1}^{k}\frac{\alpha(m_{\pi(j+1)}+\ldots+m_{\pi(k)})+\theta
m_{\pi(j)}}{(m_{\pi(j)}+\ldots+m_{\pi(k)})}$
$\displaystyle\qquad=\frac{1}{\theta}\sum_{\ell=1}^{k}\sum_{\pi\in
S_{k}^{(\ell)}}\prod_{j=1}^{k}\frac{\alpha(m_{\pi(j+1)}+\ldots+m_{\pi(k)})+\theta
m_{\pi(j)}}{(m_{\pi(j)}+\ldots+m_{\pi(k)})}$

where $S_{k}^{(\ell)}$ denotes all the permutations $\pi$ satisfying
$\pi(1)=\ell$. As a consequence one has

$\displaystyle\frac{1}{\theta}\sum_{\pi\in
S_{k}}\prod_{j=1}^{k}\frac{\alpha(m_{\pi(j+1)}+\ldots+m_{\pi(k)})+\theta
m_{\pi(j)}}{(m_{\pi(j)}+\ldots+m_{\pi(k)})}$
$\displaystyle\qquad=\frac{1}{\theta}\sum_{\ell=1}^{k}\sum_{\pi\in
S_{k}^{(\ell)}}\prod_{j=2}^{k}\frac{\alpha(m_{\pi(j+1)}+\ldots+m_{\pi(k)})+\theta
m_{\pi(j)}}{(m_{\pi(j)}+\ldots+m_{\pi(k)})}\cdot\frac{\alpha(n-m_{\ell})+\theta
m_{\ell}}{n}$
$\displaystyle\qquad=\frac{1}{\theta}\sum_{\ell=1}^{k}\frac{\alpha(n-m_{\ell})+\theta
m_{\ell}}{n}\sum_{\pi\in
S_{k}^{(\ell)}}\prod_{j=2}^{k}\frac{\alpha(m_{\pi(j+1)}+\ldots+m_{\pi(k)})+\theta
m_{\pi(j)}}{(m_{\pi(j)}+\ldots+m_{\pi(k)})}.$

We now use the inductive hypothesis for $k-1$ to evaluate the last sum over
partitions, and we obtain

$\displaystyle\frac{1}{\theta}\sum_{\pi\in
S_{k}}\prod_{j=1}^{k}\frac{\alpha(m_{\pi(j+1)}+\ldots+m_{\pi(k)})+\theta
m_{\pi(j)}}{(m_{\pi(j)}+\ldots+m_{\pi(k)})}=\sum_{\ell=1}^{k}\frac{\alpha(n-m_{\ell})+\theta
m_{\ell}}{n}\prod_{i=1}^{k-2}(\theta+i\alpha)$
$\displaystyle\qquad\qquad=\prod_{i=1}^{k-2}\cdot\left[k\alpha+\frac{1}{n}(\theta-\alpha)\sum_{j=1}^{k}m_{j}\right]=\prod_{i=1}^{k-2}(\theta+i\alpha)\cdot[\theta+\alpha(k-1)]=\prod_{i=1}^{k-1}(\theta+i\alpha)$

thus, (23) follows by induction. Now, Equation (9) is an immediate consequence
of this equality and the expression of the ordered EPPF (8).
∎

###  A.2 Proof of Equation (11)

First of all we observe that if $i=n$, then the probability
$P_{n}(n;\alpha,\theta)=1$, indeed it is the only species in the sample and it
is both the youngest and the oldest. We now assume that $i<n$ and we would
like to evaluate the probability of the following event

$\mathscr{E}:=\\{\text{a species represented $i$ times in a sample of size $n$
has highest weight}\\}.$

In order to evaluate $\text{Pr}[\mathscr{E}]$, we denote by $\tilde{\pi}_{n}$
the random partition generated by the sample of size $n$, which will be
characterized by a species having frequency $i$, while the remaining $n-i$
observations are partitioned into $K_{n-i}\geq 1$ distinct values with
frequencies $w_{1},\ldots,w_{k}$. Note that $\tilde{\pi}_{n}$ does not take
into account the ordering of the species. Thus we get

$\text{Pr}[\mathscr{E}]=\mathds{E}[\text{Pr}[\mathscr{E}|\tilde{\pi}_{n}]]=\mathds{E}\left[\frac{\sum_{\pi\in
S_{K_{n-i}}}\Phi_{K_{n-i}+1}^{(n)}(i,w_{\pi(1)},\ldots,w_{\pi(K_{n-i})})}{\frac{(1-\alpha)_{(i-1)}\prod_{j=1}^{K_{n-i}}(1-\alpha)_{(w_{j}-1)}}{(\theta+1)_{(n-1)}}\prod_{s=1}^{K_{n-i}+1-1}(\theta+\alpha
s)}\right]$ (24)

where the denominator is nothing but the EPPF of the PYP, whereas the
numerator is the sum of all the order EPPF of the PYP in which the oldest
species equals a specific species having frequency $i$. Note that the sum in
(24) is over all possible permutations of the remaining $K_{n-i}$ distinct
species, thus the probability $\text{Pr}[\mathscr{E}]$ equals

$\displaystyle\mathds{E}\left[\frac{\frac{(1-\alpha)_{(i-1)}\prod_{j=1}^{K_{n-i}}(1-\alpha)_{(w_{j}-1)}}{(\theta+1)_{(n-1)}}\sum_{\pi\in
S_{K_{n-i}}}\theta^{-1}\prod_{j=1}^{K_{n-i}}\frac{\alpha(w_{\pi(j+1)}+\ldots+w_{\pi(K_{n-i})})+\theta
w_{\pi(j)}}{(w_{\pi(j)}+\ldots+w_{\pi(K_{n-i})})}\cdot\frac{\alpha(n-i)+\theta
i}{n}}{\frac{(1-\alpha)_{(i-1)}\prod_{j=1}^{K_{n-i}}(1-\alpha)_{(w_{j}-1)}}{(\theta+1)_{(n-1)}}\prod_{s=1}^{K_{n-i}}(\theta+\alpha
s)}\right]$ $\displaystyle=\mathds{E}\left[\frac{\sum_{\pi\in
S_{K_{n-i}}}\theta^{-1}\prod_{j=1}^{K_{n-i}}\frac{\alpha(w_{\pi(j+1)}+\ldots+w_{\pi(K_{n-i})})+\theta
w_{\pi(j)}}{(w_{\pi(j)}+\ldots+w_{\pi(K_{n-i})})}\cdot\frac{\alpha(n-i)+\theta
i}{n}}{\prod_{s=1}^{K_{n-i}}(\theta+\alpha s)}\right]$
$\displaystyle=\mathds{E}\left[\frac{\prod_{s=1}^{K_{n-i}-1}(\theta+\alpha
s)}{\prod_{s=1}^{K_{n-i}}(\theta+\alpha s)}\cdot\frac{\alpha(n-i)+\theta
i}{n}\right]=\frac{\alpha
n+i(\theta-\alpha)}{n}\cdot\mathds{E}\left[\frac{1}{\theta+\alpha
K_{n-i}}\right]$

where we used the explicit expression of the ordered EPPF and the
marginalization result (23). Note that the expected value in the last equation
is made w.r.t. the distribution of $K_{n-i}$ (number of distinct types in a
sample of size $n-i$).
∎

###  A.3 Distributions of the number of distinct species

In this section we remind the distributions of $K_{n}$ and of $K_{m}^{(n)}$,
which has been used in proofs and in the main paper. The distribution of the
number of distinct species $K_{n}$ is contained in Pitman (2006, Chapter 3).
Pitman (2006) has shown that that for any $k=1,\ldots,n$

$\text{Pr}[K_{n}=k]=\frac{\left(\frac{\theta}{\alpha}\right)_{(k)}}{(\theta)_{(n)}}\mathscr{C}(n,k;\alpha),$
(25)

where $\mathscr{C}(n,k;\alpha)=(k!)^{-1}\sum_{0\leq i\leq k}{k\choose
i}(-1)^{i}(-i\alpha)_{(n)}$ denotes the generalized factorial coefficient,
with the proviso that $\mathscr{C}(0,0;\alpha)=1$ and
$\mathscr{C}(m,0;\alpha)=0$ (Charalambides, 2005, Chapter 3). See also Pitman
(2006, Chapter 2) for details.

Now we remind the posterior distribution of $K_{m}^{(n)}=K_{n+m}-K_{n}$, that
has been determined in Lijoi et al. (2007, Proposition 1). For any
$s=0,\ldots,m$ one has

$\text{Pr}[K_{m}^{(n)}=s\,|\,K_{n}=k]=\frac{\left(k+\frac{\theta}{\alpha}\right)_{(s)}}{(\theta+n)_{(m)}}\frac{\mathscr{C}(m,s;\alpha;-n+k\alpha)}{\alpha^{s}}$
(26)

where $\mathscr{C}(n,k;\alpha,\gamma)=(k!)^{-1}\sum_{0\leq i\leq k}{k\choose
i}(-1)^{i}(-i\alpha-\gamma)_{(n)}$ denotes the non-centered (shifted)
generalized factorial coefficient, with the proviso that
$\mathscr{C}(0,0;\alpha,\gamma)=1$ and
$\mathscr{C}(n,0;\alpha,\gamma)=(-\gamma)_{(n)}$ (Charalambides, 2005, Chapter
3).

###  A.4 Proof of Equation (13)

Recall the ordered EPPF in (8), thus we can evaluate the following probability

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}=k]=\sum_{(\star)}\binom{n}{m_{1},\ldots,m_{k}}\Phi_{k}^{(n)}(m_{1},\ldots,m_{k})$
$\displaystyle\qquad\sum_{(\star)}\binom{n}{m_{1},\ldots,m_{r},m_{r+1},\ldots,m_{k}}\Phi_{k}^{(n)}(m_{1},\ldots,m_{r},m_{r+1},\ldots,m_{k})$

where the sum runs over all the vector of positive integers
$(m_{r+1},\ldots,m_{k})$ with $m_{i}\geq 1$ and
$m_{r+1}+\cdots+m_{k}=n-m_{1}-\cdots-m_{r}$. Using the explicit expression of
the previous ordered EPPF we get

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}=k]$
$\displaystyle\qquad=\sum_{(\star)}\binom{n}{m_{1},\ldots,m_{r},m_{r+1},\ldots,m_{k}}\frac{\prod_{j=r+1}^{k}(1-\alpha)_{(m_{j}-1)}\prod_{j=1}^{r}(1-\alpha)_{(m_{j}-1)}}{(\theta+1)_{(n-1)}}$
$\displaystyle\qquad\qquad\times\theta^{-1}\prod_{j=r+1}^{k}\frac{\alpha
r_{j+1}+\theta m_{j}}{r_{j}}\frac{\prod_{j=1}^{r}[\alpha(n-m_{1}-\cdots-
m_{j})+\theta m_{j}]}{n(n-m_{1})\cdots(n-m_{1}-\cdots-m_{r-1})}$
$\displaystyle\quad=\frac{n!}{(\theta+1)_{(n-1)}}\frac{1}{(n-m_{1}-\cdots-
m_{r})!}\frac{\prod_{j=1}^{r}[\alpha(n-m_{1}-\cdots-m_{j})+\theta
m_{j}]}{n(n-m_{1})\cdots(n-m_{1}-\cdots-
m_{r-1})}\cdot\prod_{j=1}^{r}\frac{(1-\alpha)_{(m_{j}-1)}}{m_{j}!}$
$\displaystyle\qquad\qquad\times\sum_{(\star)}\binom{n-m_{1}-\cdots-
m_{r}}{m_{r+1},\ldots,m_{k}}\theta^{-1}\prod_{j=r+1}^{k}(1-\alpha)_{(m_{j}-1)}\prod_{j=r+1}^{k}\frac{\alpha(m_{j+1}+\cdots+m_{k})+\theta
m_{j}}{m_{j}+\cdots+m_{k}}.$

If we divide the last sum by the coefficient $(\theta+1)_{(n-m_{1}-\ldots-
m_{r}-1)}$, we can recognize the probability $\text{Pr}[K_{n-m_{1}-\ldots-
m_{r}}=k-r]$, thus one has

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}=k]=\frac{n!}{(\theta+1)_{(n-1)}}\frac{(\theta+1)_{(n-m_{1}-\ldots-
m_{r}-1)}}{(n-m_{1}-\cdots-m_{r})!}$
$\displaystyle\qquad\times\frac{\prod_{j=1}^{r}[\alpha(n-m_{1}-\cdots-
m_{j})+\theta m_{j}]}{n(n-m_{1})\cdots(n-m_{1}-\cdots-
m_{r-1})}\cdot\prod_{j=1}^{r}\frac{(1-\alpha)_{(m_{j}-1)}}{m_{j}!}\text{Pr}[K_{n-m_{1}-\cdots-
m_{r}}=k-r].$

Now we can exploit (25) to get

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}=k]=\frac{n!}{(\theta+1)_{(n-1)}}\frac{\mathscr{C}(n-m_{1}-\ldots-
m_{r},k-r;\alpha)}{(n-m_{1}-\cdots-m_{r})!}$
$\displaystyle\qquad\times\frac{\prod_{j=1}^{r}[\alpha(n-m_{1}-\cdots-
m_{j})+\theta m_{j}]}{n(n-m_{1})\cdots(n-m_{1}-\cdots-
m_{r-1})}\cdot\prod_{j=1}^{r}\frac{(1-\alpha)_{(m_{j}-1)}}{m_{j}!}\frac{\prod_{i=1}^{k-r-1}(\theta+i\alpha)}{\alpha^{k-r}}.$

By dividing the previous expression by the probability $\text{Pr}[K_{n}=k]$ in
(25), one may now obtain the following conditional distribution

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r}\,|\,K_{n}=k]$ (27)
$\displaystyle\quad={n\choose
m_{1},\ldots,m_{r},n-|\mathbf{m}|_{1:r}}\frac{\displaystyle\prod_{j=1}^{r}\frac{[\alpha(n-|\mathbf{m}|_{1:j})+\theta
m_{j}](1-\alpha)_{(m_{j}-1)}}{n-|\mathbf{m}|_{1:j-1}}}{\displaystyle\frac{\mathscr{C}(n,k;\alpha)}{\left(\frac{\theta}{\alpha}+k-r\right)_{(r)}\mathscr{C}(n-|\mathbf{m}|_{1:r},k-r;\alpha)}}.$

In order to obtain (13), note that:

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}\geq
r]=\sum_{k=r}^{s}\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}=k]$

where $s=r+n-m_{1}-\cdots-m_{r}$, thus, we get

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}\geq
r]=\sum_{k=r}^{s}\frac{n!}{(\theta+1)_{(n-1)}}\frac{(\theta+1)_{(n-m_{1}-\ldots-
m_{r}-1)}}{(n-m_{1}-\cdots-m_{r})!}$
$\displaystyle\qquad\times\frac{\prod_{j=1}^{r}[\alpha(n-m_{1}-\cdots-
m_{j})+\theta m_{j}]}{n(n-m_{1})\cdots(n-m_{1}-\cdots-
m_{r-1})}\cdot\prod_{j=1}^{r}\frac{(1-\alpha)_{(m_{j}-1)}}{m_{j}!}\text{Pr}[K_{n-m_{1}-\cdots-
m_{r}}=k-r]$

where we observe that the sum over $k$ of the probability
$\text{Pr}[K_{n-m_{1}-\cdots-m_{r}}=k-r]$ is equal to $1$, as a consequence we
obtain

$\displaystyle\text{Pr}[M_{1,n}=m_{1},\ldots,M_{r,n}=m_{r},K_{n}\geq
r]=\frac{n!}{(\theta+1)_{(n-1)}}\frac{(\theta+1)_{(n-m_{1}-\ldots-
m_{r}-1)}}{(n-m_{1}-\cdots-m_{r})!}$
$\displaystyle\qquad\times\frac{\prod_{j=1}^{r}[\alpha(n-m_{1}-\cdots-
m_{j})+\theta m_{j}]}{n(n-m_{1})\cdots(n-m_{1}-\cdots-
m_{r-1})}\cdot\prod_{j=1}^{r}\frac{(1-\alpha)_{(m_{j}-1)}}{m_{j}!}$

and the thesis follows.
∎

###  A.5 Proof of Theorem 1

We first focus on the proof of (17). To this end, evaluate

$\displaystyle\text{Pr}[A_{r},W_{1,n+m}=w_{1},\ldots,W_{r,n+m}=w_{r},K_{m}^{(n)}=s|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\frac{\text{Pr}[A_{r},W_{1,n+m}=w_{1},\ldots,W_{r,n+m}=w_{r},K_{m}^{(n)}=s,K_{n}=k,M_{i,n}=m_{i},\;i=1,\ldots,k]}{\text{Pr}[K_{n}=k,M_{i,n}=m_{i},\;i=1,\ldots,k]}$

where the numerator may be evaluated along similar lines as in Section A.4 of
the Supplementary Material. In fact one has

$\displaystyle\text{Pr}[A_{r},W_{1,n+m}=w_{1},\ldots,W_{r,n+m}=w_{r},K_{m}^{(n)}=s,K_{n}=k,M_{i,n}=m_{i},\;i=1,\ldots,k]$
$\displaystyle\qquad=\sum_{(\star)}\binom{n}{n_{r+1},\ldots,n_{k+s}}\binom{m}{w_{1},\ldots,w_{r},w_{r+1},\ldots,w_{k+s}}\frac{\prod_{j=r+1}^{k+s}(1-\alpha)_{(w_{j}+n_{j}-1)}}{(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\qquad\times\prod_{j=1}^{r}(1-\alpha)_{(w_{j}-1)}\theta^{-1}\prod_{j=r+1}^{k+s}\frac{\alpha(w_{j+1}+n_{j+1}+\cdots+w_{k+s}+n_{k+s})+\theta(w_{j}+n_{j})}{w_{j}+n_{j}+\cdots+w_{k+s}+n_{k+s}}$
$\displaystyle\qquad\qquad\qquad\times\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-\cdots-
w_{j})+\theta w_{j})}{(n+m)(n+m-w_{1})\cdots(n+m-w_{1}-\cdots-w_{r-1})}$

where the sum runs over all the vectors $(w_{r+1},\ldots,w_{k+s})$ and
$(n_{r+1},\ldots,n_{k+s})$ such that $w_{i}\geq 0$ and
$w_{r+1}+\cdots+w_{k+s}=m-w_{1}-\cdots-w_{r}$, whereas
$n_{r+1}+\cdots+n_{k+s}=n$ and $n_{i_{j}}=m_{j}$ for $j=1,\ldots,k$ and
$i_{1}<\cdots<i_{k}$ the other values of $n_{i}$ are equal to $0$; finally
$w_{i}+n_{i}\geq 1$ for any $i=r+1,\ldots,k+s$. Thus

$\displaystyle\text{Pr}[A_{r},W_{1,n+m}=w_{1},\ldots,W_{r,n+m}=w_{r},K_{m}^{(n)}=s,K_{n}=k,M_{i,n}=m_{i},\;i=1,\ldots,k]$
$\displaystyle\qquad=\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-\cdots-
w_{j})+\theta w_{j})}{(n+m)(n+m-w_{1})\cdots(n+m-w_{1}-\cdots-
w_{r-1})}\prod_{j=1}^{r}\frac{(1-\alpha)_{(w_{j}-1)}}{w_{j}!}$
$\displaystyle\qquad\times\frac{m!}{(m-w_{1}-\cdots-
w_{r})!(\theta+1)_{(n+m-1)}}\sum_{(\star)}\binom{n}{n_{r+1},\ldots,n_{k+s}}\binom{m-w_{1}-\cdots-
w_{r}}{w_{r+1},\ldots,w_{k+s}}$
$\displaystyle\qquad\qquad\times\theta^{-1}\prod_{j=r+1}^{k+s}(1-\alpha)_{(w_{j}+n_{j}-1)}\prod_{j=r+1}^{k+s}\frac{\alpha(w_{j+1}+n_{j+1}+\cdots+w_{k+s}+n_{k+s})+\theta(w_{j}+n_{j})}{w_{j}+n_{j}+\cdots+w_{k+s}+n_{k+s}}$

if we divide the last sum by the coefficient $(\theta+1)_{(n+m-w_{1}-\cdots-
w_{r}-1)}$ we obtain the following

$\begin{split}&\frac{1}{(\theta+1)_{(n+m-w_{1}-\cdots-
w_{r}-1)}}\sum_{(\star)}\binom{n}{n_{r+1},\ldots,n_{k+s}}\binom{m-w_{1}-\cdots-
w_{r}}{w_{r+1},\ldots,w_{k+s}}\\\
&\quad\times\theta^{-1}\prod_{j=r+1}^{k+s}(1-\alpha)_{(w_{j}+n_{j}-1)}\\\
&\qquad\qquad\times\prod_{j=r+1}^{k+s}\frac{\alpha(w_{j+1}+n_{j+1}+\cdots+w_{k+s}+n_{k+s})+\theta(w_{j}+n_{j})}{w_{j}+n_{j}+\cdots+w_{k+s}+n_{k+s}}\\\
&\qquad=\text{Pr}[K_{m-w_{1}-\cdots-
w_{r}}^{(n)}=s-r|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]\\\
&\qquad\qquad\times\text{Pr}[K_{n}=k,\mathbf{M}_{n}=\mathbf{m}].\end{split}$
(28)

Substituting (28) in the expression of the joint probability under study, we
obtain

$\displaystyle\text{Pr}[A_{r},W_{1,n+m}=w_{1},\ldots,W_{r,n+m}=w_{r},K_{m}^{(n)}=s,K_{n}=k,M_{i,n}=m_{i},\;i=1,\ldots,k]$
$\displaystyle\qquad=\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-\cdots-
w_{j})+\theta w_{j})}{(n+m)(n+m-w_{1})\cdots(n+m-w_{1}-\cdots-
w_{r-1})}\prod_{j=1}^{r}\frac{(1-\alpha)_{(w_{j}-1)}}{w_{j}!}$
$\displaystyle\qquad\qquad\times\frac{m!(\theta+1)_{(n+m-w_{1}-\cdots-
w_{r}-1)}}{(m-w_{1}-\cdots-w_{r})!(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\qquad\qquad\times\text{Pr}[K_{m-w_{1}-\cdots-
w_{r}}^{(n)}=s-r|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad\qquad\qquad\qquad\times\text{Pr}[K_{n}=k,\mathbf{M}_{n}=\mathbf{m}].$

As a consequence:

$\displaystyle\text{Pr}[A_{r},W_{1,n+m}=w_{1},\ldots,W_{r,n+m}=w_{r},K_{m}^{(n)}=s|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-\cdots-
w_{j})+\theta w_{j})}{(n+m)(n+m-w_{1})\cdots(n+m-w_{1}-\cdots-
w_{r-1})}\prod_{j=1}^{r}\frac{(1-\alpha)_{(w_{j}-1)}}{w_{j}!}$
$\displaystyle\qquad\qquad\times\frac{m!(\theta+1)_{(n+m-w_{1}-\cdots-
w_{r}-1)}}{(m-w_{1}-\cdots-w_{r})!(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\qquad\qquad\times\text{Pr}[K_{m-w_{1}-\cdots-
w_{r}}^{(n)}=s-r|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$

where the last probability can be evaluated resorting to the posterior
distribution of $K_{m-w_{1}-\cdots-w_{r}}^{(n)}$ in (26). Summing over
$s=r,\ldots,r+m-w_{1}-\cdots-w_{r}$ the previous expression, we obtain (17).

We now move to the proof of (18). Proceeding along the same lines as in the
first part of the proof, we have

$\displaystyle\text{Pr}[B_{r},W_{1,n+m}=w_{1}+m_{1},\ldots,W_{r,n+m}=w_{r}+m_{r},K_{m}^{(n)}=s|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\quad=\frac{\text{Pr}[B_{r},W_{1,n+m}=w_{1}+m_{1},\ldots,W_{r,n+m}=w_{r}+m_{r},K_{m}^{(n)}=s,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]}{\text{Pr}[K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]}$

where the numerator may be evaluated along similar lines as in the proof of
Proposition 1. In fact, as for the numerator, one has

$\displaystyle\text{Pr}[B_{r},W_{1,n+m}=w_{1}+m_{1},\ldots,W_{r,n+m}=w_{r}+m_{r},K_{m}^{(n)}=s,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\sum_{(\star)}\binom{n}{m_{1},\ldots,m_{r},n_{r+1},\ldots,n_{k+s}}\binom{m}{w_{1},\ldots,w_{r},w_{r+1},\ldots,w_{k+s}}\frac{\prod_{j=r+1}^{k+s}(1-\alpha)_{(w_{j}+n_{j}-1)}}{(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\qquad\times\prod_{j=1}^{r}(1-\alpha)_{(w_{j}+m_{j}-1)}\theta^{-1}\prod_{j=r+1}^{k+s}\frac{\alpha(w_{j+1}+n_{j+1}+\cdots+w_{k+s}+n_{k+s})+\theta(w_{j}+n_{j})}{w_{j}+n_{j}+\cdots+w_{k+s}+n_{k+s}}$
$\displaystyle\qquad\qquad\qquad\times\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-m_{1}-\cdots-
w_{j}-m_{j})+\theta(w_{j}+m_{j}))}{(n+m)(n+m-w_{1}-m_{1})\cdots(n+m-w_{1}-m_{1}-\cdots-
w_{r-1}-m_{r-1})}$

where the sum runs over all the vectors $(w_{r+1},\ldots,w_{k+s})$ and
$(n_{r+1},\ldots,n_{k+s})$ such that $w_{i}\geq 0$ and
$w_{r+1}+\cdots+w_{k+s}=m-w_{1}-\cdots-w_{r}$, whereas
$n_{r+1}+\cdots+n_{k+s}=n-m_{1}-\cdots-m_{r}$ and $n_{i_{j}}=m_{j}$ for
$j=1,\ldots,k-r$ and $i_{1}<\cdots<i_{k-r}$ the other values of $n_{i}$ are
equal to $0$; finally $w_{i}+n_{i}\geq 1$ for any $i=r+1,\ldots,k+s$. Thus

$\displaystyle\text{Pr}[B_{r},W_{1,n+m}=w_{1}+m_{1},\ldots,W_{r,n+m}=w_{r}+m_{r},K_{m}^{(n)}=s,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-m_{1}-\cdots-
w_{j}-m_{j})+\theta(w_{j}+m_{j}))}{(n+m)(n+m-w_{1}-m_{1})\cdots(n+m-w_{1}-m_{1}-\cdots-
w_{r-1}-m_{r-1})}$
$\displaystyle\qquad\times\prod_{j=1}^{r}\frac{(1-\alpha)_{(w_{j}+m_{j}-1)}}{w_{j}!m_{j}!}\frac{m!n!}{(m-w_{1}-\cdots-
w_{r})!(n-m_{1}-\cdots-m_{r})!(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\qquad\times\sum_{(\star)}\binom{n-m_{1}-\cdots-
m_{r}}{n_{r+1},\ldots,n_{k+s}}\binom{m-w_{1}-\cdots-
w_{r}}{w_{r+1},\ldots,w_{k+s}}\theta^{-1}\prod_{j=r+1}^{k+s}(1-\alpha)_{(w_{j}+n_{j}-1)}$
$\displaystyle\qquad\qquad\qquad\times\prod_{j=r+1}^{k+s}\frac{\alpha(w_{j+1}+n_{j+1}+\cdots+w_{k+s}+n_{k+s})+\theta(w_{j}+n_{j})}{w_{j}+n_{j}+\cdots+w_{k+s}+n_{k+s}}$

if we divide the last sum by the coefficient
$(\theta+1)_{(n+m-w_{1}-m_{1}-\cdots-w_{r}-m_{r}-1)}$ we obtain the following

$\begin{split}&\frac{1}{(\theta+1)_{(n+m-w_{1}-m_{1}-\cdots-
w_{r}-m_{r}-1)}}\sum_{(\star)}\binom{n-m_{1}-\cdots-
m_{r}}{n_{r+1},\ldots,n_{k+s}}\binom{m-w_{1}-\cdots-
w_{r}}{w_{r+1},\ldots,w_{k+s}}\\\
&\quad\times\theta^{-1}\prod_{j=r+1}^{k+s}(1-\alpha)_{(w_{j}+n_{j}-1)}\\\
&\qquad\qquad\times\prod_{j=r+1}^{k+s}\frac{\alpha(w_{j+1}+n_{j+1}+\cdots+w_{k+s}+n_{k+s})+\theta(w_{j}+n_{j})}{w_{j}+n_{j}+\cdots+w_{k+s}+n_{k+s}}\\\
&\qquad=\text{Pr}[K_{m-w_{1}-\cdots-w_{r}}^{(n-m_{1}-\cdots-
m_{r})}=s|K_{n-m_{1}-\cdots-m_{r}}=k-r,M_{i,n}=m_{i},\;i=r+1,\ldots,k]\\\
&\qquad\qquad\times\text{Pr}[K_{n-m_{1}-\cdots-
m_{r}}=k-r,M_{i,n}=m_{i},\;i=r+1,\ldots,k].\end{split}$ (29)

Substituting (29) in the expression of the joint probability under study, we
obtain

$\displaystyle\text{Pr}[B_{r},W_{1,n+m}=w_{1}+m_{1},\ldots,W_{r,n+m}=w_{r}+m_{r},K_{m}^{(n)}=s,K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-m_{1}-\cdots-
w_{j}-m_{j})+\theta(w_{j}+m_{j}))}{(n+m)(n+m-w_{1}-m_{1})\cdots(n+m-w_{1}-m_{j}-\cdots-
w_{r-1}-m_{r-1})}$
$\displaystyle\qquad\quad\times\prod_{j=1}^{r}\frac{(1-\alpha)_{(w_{j}+m_{j}-1)}}{w_{j}!m_{j}!}\cdot\frac{m!n!(\theta+1)_{(n+m-w_{1}-m_{1}-\cdots-
w_{r}-m_{r}-1)}}{(m-w_{1}-\cdots-w_{r})!(n-m_{1}-\cdots-
m_{r})!(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\quad\quad\qquad\times\text{Pr}[K_{m-w_{1}-\cdots-
w_{r}}^{(n-m_{1}-\cdots-m_{r})}=s|K_{n-m_{1}-\cdots-
m_{r}}=k-r,M_{i,n}=m_{i},\;i=r+1,\ldots,k]$
$\displaystyle\qquad\quad\quad\quad\qquad\times\text{Pr}[K_{n-m_{1}-\cdots-
m_{r}}=k-r,M_{i,n}=m_{i},\;i=r+1,\ldots,k].$

As a consequence:

$\displaystyle\text{Pr}[B_{r},W_{1,n+m}=w_{1}+m_{1},\ldots,W_{r,n+m}=w_{r}+m_{r},K_{m}^{(n)}=s|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-m_{1}-\cdots-
w_{j}-m_{j})+\theta(w_{j}+m_{j}))}{(n+m)(n+m-w_{1}-m_{1})\cdots(n+m-w_{1}-m_{j}-\cdots-
w_{r-1}-m_{r-1})}$
$\displaystyle\qquad\quad\times\prod_{j=1}^{r}\frac{(1-\alpha)_{(w_{j}+m_{j}-1)}}{w_{j}!m_{j}!}\cdot\frac{m!n!(\theta+1)_{(n+m-w_{1}-m_{1}-\cdots-
w_{r}-m_{r}-1)}}{(m-w_{1}-\cdots-w_{r})!(n-m_{1}-\cdots-
m_{r})!(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\quad\quad\qquad\times\text{Pr}[K_{m-w_{1}-\cdots-
w_{r}}^{(n-m_{1}-\cdots-m_{r})}=s|K_{n-m_{1}-\cdots-
m_{r}}=k-r,M_{i,n}=m_{i},\;i=r+1,\ldots,k]$
$\displaystyle\qquad\quad\quad\quad\qquad\times\frac{\text{Pr}[K_{n-m_{1}-\cdots-
m_{r}}=k-r,M_{i,n}=m_{i},\;i=r+1,\ldots,k]}{\text{Pr}[K_{n}=k,M_{i,n}=m_{i},\;i=1,\ldots,k]}.$

Thus we have to evaluate the last ratio in the previous expression, this
calculation can be easily addressed resorting to the available expression of
the ordered EPPF for the PYP (8). Indeed we get:

$\displaystyle\frac{\text{Pr}[K_{n-m_{1}-\cdots-
m_{r}}=k-r,M_{i,n}=m_{i},\;i=r+1,\ldots,k]}{\text{Pr}[K_{n}=k,M_{i,n}=m_{i},\;i=1,\ldots,k]}=\frac{\binom{n-m_{1}-\cdots-
m_{r}}{m_{r+1},\ldots,m_{k}}}{\binom{n}{m_{1},\ldots,m_{k}}}$
$\displaystyle\qquad\times\frac{(\theta+1)_{(n-1)}}{(\theta+1)_{(n-m_{1}-\cdots-
m_{r}-1)}}\cdot\frac{1}{\prod_{j=1}^{r}(1-\alpha)_{(m_{j}-1)}}\cdot\frac{n(n-m_{1})\cdots(n-m_{1}-\cdots-
m_{r-1})}{\prod_{j=1}^{r}(\alpha(n-m_{1}-\cdots-m_{j})+\theta m_{j})}.$

Substituting the previous equality in the expression of the conditional
probability under study we get:

$\displaystyle\text{Pr}[B_{r},W_{1,n+m}=w_{1}+m_{1},\ldots,W_{r,n+m}=w_{r}+m_{r},K_{m}^{(n)}=s|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\frac{\prod_{j=1}^{r}(\alpha(n+m-w_{1}-m_{1}-\cdots-
w_{j}-m_{j})+\theta(w_{j}+m_{j}))}{(n+m)(n+m-w_{1}-m_{1})\cdots(n+m-w_{1}-m_{j}-\cdots-
w_{r-1}-m_{r-1})}$
$\displaystyle\qquad\quad\times\prod_{j=1}^{r}\frac{(1-\alpha)_{(w_{j}+m_{j}-1)}}{w_{j}!(1-\alpha)_{(m_{j}-1)}}\cdot\frac{m!(\theta+1)_{(n+m-w_{1}-m_{1}-\cdots-
w_{r}-m_{r}-1)}(\theta+1)_{(n-1)}}{(m-w_{1}-\cdots-
w_{r})!(\theta+1)_{(n+m-1)}(\theta+1)_{(n-m_{1}-\cdots-m_{r}-1)}}$
$\displaystyle\qquad\quad\quad\qquad\times\text{Pr}[K_{m-w_{1}-\cdots-
w_{r}}^{(n-m_{1}-\cdots-m_{r})}=s|K_{n-m_{1}-\cdots-
m_{r}}=k-r,M_{i,n}=m_{i},\;i=r+1,\ldots,k]$
$\displaystyle\qquad\quad\quad\quad\qquad\times\frac{n(n-m_{1})\cdots(n-m_{1}-\cdots-
m_{r-1})}{\prod_{j=1}^{r}(\alpha(n-m_{1}-\cdots-m_{j})+\theta m_{j})}.$

Summing overall the possible values of $s\in\\{0,\ldots,m-w_{1}-\cdots-
w_{r}\\}$ we finally get the expression (18).
∎

###  A.6 Proof of Equation (22)

We focus on the proof of (22) by evaluating the two expected values

$\mathds{E}[W_{1,n+m}\mathds{1}_{A_{1}}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
(30)

and

$\mathds{E}[W_{1,n+m}\mathds{1}_{B_{1}}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}].$
(31)

To calculate (30), we need to evaluate the posterior expected value of
$W_{1,n+m}$ on the event $A_{1}$:

$\displaystyle\mathds{E}[W_{1,n+m}\mathds{1}_{A_{1}}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\sum_{w=1}^{m}w\binom{m}{w}\frac{\alpha(n+m-w)+\theta
w}{n+m}\cdot\frac{(1-\alpha)_{(w-1)}(\theta+1)_{(n+m-w-1)}}{(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad=\frac{m}{(n+m)(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\qquad\times\sum_{w=0}^{m-1}\binom{m-1}{w}[\alpha(n+m-w-1)+\theta(w+1)](1-\alpha)_{(w)}(\theta+1)_{(n+m-w-2)}$
$\displaystyle\qquad=\frac{m}{(n+m)\Gamma(\theta+n+m)\Gamma(1-\alpha)}(S_{1}+S_{2})$
(32)

where $S_{1}$ and $S_{2}$ are the two sums defined as follows

$\displaystyle
S_{1}:=\sum_{w=0}^{m-1}\binom{m-1}{w}[\alpha(n+m-1)+\theta]\Gamma(w-\alpha+1)\Gamma(\theta+n+m-q-1)$
(33) $\displaystyle
S_{2}:=\sum_{w=0}^{m-1}\binom{m-1}{w}w(\theta-\alpha)\Gamma(w-\alpha+1)\Gamma(\theta+n+m-w-1).$
(34)

We now concentrate on the two sums separately. In order to evaluate the first
one (33) we exploit the integral representation of the beta function:

$\displaystyle S_{1}$
$\displaystyle=[\alpha(n+m-1)+\theta]\Gamma(\theta+n+m-\alpha)\sum_{w=0}^{m-1}\binom{m-1}{w}B(w-\alpha+1,\theta+n+m-w-1)$
$\displaystyle=[\alpha(n+m-1)+\theta]\Gamma(\theta+n+m-\alpha)$
$\displaystyle\qquad\qquad\times\int_{0}^{1}x^{-\alpha}(1-x)^{\theta+n+m-2}\sum_{w=0}^{m-1}\binom{m-1}{w}\frac{x^{w}}{(1-x)^{w}}\mathrm{d}x$
$\displaystyle=[\alpha(n+m-1)+\theta]\Gamma(\theta+n+m-\alpha)\int_{0}^{1}x^{-\alpha}(1-x)^{\theta+n+m-2}\Big{(}1+\frac{x}{1-x}\Big{)}^{m-1}\mathrm{d}x$
$\displaystyle=[\alpha(n+m-1)+\theta]\Gamma(\theta+n+m-\alpha)\Gamma(\theta+n+m-\alpha)B(1-\alpha,\theta+n).$

The evaluation of the sum (34) proceeds in a similar fashion. We first observe
that for $w=0$ the summand in $S_{2}$ is zero, thus we can write:

$\displaystyle S_{2}$
$\displaystyle=(\theta-\alpha)(m-1)\sum_{w=0}^{m-2}\binom{m-2}{w}\Gamma(w-\alpha+2)\Gamma(\theta+n+m-w-2)$
$\displaystyle=(\theta-\alpha)(m-1)\Gamma(\theta+n+m-\alpha)\sum_{w=0}^{m-2}\binom{m-2}{w}B(w-\alpha+2,\theta+n+m-w-2)$
$\displaystyle=(\theta-\alpha)(m-1)\Gamma(\theta+n+m-\alpha)\int_{0}^{1}\sum_{w=0}^{m-2}\binom{m-2}{w}x^{w-\alpha+1}(1-x)^{\theta+n+m-w-3}\mathrm{d}x$
$\displaystyle=(\theta-\alpha)(m-1)\Gamma(\theta+n+m-\alpha)\int_{0}^{1}x^{-\alpha+1}(1-x)^{\theta+n+m-3}\Big{(}1+\frac{x}{1-x}\Big{)}^{m-2}\mathrm{d}x$
$\displaystyle=(\theta-\alpha)(m-1)\Gamma(\theta+n+m-\alpha)B(2-\alpha,\theta+n).$

We can now put together the expressions we determined for $S_{1}$ and $S_{2}$
in (32) to get

$\displaystyle\mathds{E}[W_{1,n+m}\mathds{1}_{A_{1}}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\frac{m\Gamma(\theta+n+m-\alpha)\Gamma(\theta+n)}{(n+m)\Gamma(\theta+n+1-\alpha)\Gamma(\theta+n+m)}$
$\displaystyle\qquad\qquad\qquad\times\left\\{\alpha(n+m-1)+\theta+\frac{(\theta-\alpha)(m-1)(1-\alpha)}{\theta+n+1-\alpha}\right\\}$
$\displaystyle\qquad=\frac{m}{n+m}\frac{(\theta+n-\alpha)_{(m)}}{(\theta+n)_{(m)}(\theta+n-\alpha)}\left\\{\alpha
n+\theta+\alpha(m-1)+\frac{(\theta-\alpha)(m-1)(1-\alpha)}{\theta+n+1-\alpha}\right\\}$
$\displaystyle\qquad=\frac{m}{n+m}\frac{(\theta+n-\alpha)_{(m)}}{(\theta+n)_{(m)}(\theta+n-\alpha)}\cdot(\theta+n\alpha)\frac{\theta+n+m-\alpha}{\theta+n+1-\alpha}$

where the last equality follows from straightforward calculations. Thus the
expression of (30) follows

$\mathds{E}[W_{1,n+m}\mathds{1}_{A_{1}}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]=\frac{m}{n+m}\cdot\frac{\theta+n\alpha}{\theta+n+1-\alpha}\frac{(\theta+n+1-\alpha)_{(m)}}{(\theta+n)_{(m)}}.$
(35)

We now concentrate on the evaluation of (31), which follows in a similar
manner.

$\displaystyle\mathds{E}[W_{1,n+m}\mathds{1}_{B_{1}}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\sum_{w=m_{1}}^{m_{1}+m}w\frac{\alpha(n+m-w)+\theta
w}{n+m}\cdot\frac{(1-\alpha)_{(w-1)}(\theta+1)_{(n+m-w-1)}}{(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\qquad\qquad\times\frac{n}{\alpha(n-m_{1})+\theta
m_{1}}\binom{m}{w-m_{1}}\frac{(\theta+1)_{(n-1)}}{(1-\alpha)_{(m_{1}-1)}(\theta+1)_{(n-m_{1}-1)}}$
$\displaystyle\qquad=\sum_{w=0}^{m}(w+m_{1})\frac{\alpha(n+m-w-
m_{1})+\theta(w+m_{1})}{n+m}\cdot\frac{(1-\alpha)_{(w+m_{1}-1)}(\theta+1)_{(n+m-w-
m_{1}-1)}}{(\theta+1)_{(n+m-1)}}$
$\displaystyle\qquad\qquad\qquad\times\frac{n}{\alpha(n-m_{1})+\theta
m_{1}}\binom{m}{w}\frac{(\theta+1)_{(n-1)}}{(1-\alpha)_{(m_{1}-1)}(\theta+1)_{(n-m_{1}-1)}}$
$\displaystyle\qquad=\frac{n}{(n+m)(\alpha(n-m_{1})+\theta
m_{1})}\cdot\frac{(\theta+1)_{(n-1)}}{(\theta+1)_{(n+m-1)}(1-\alpha)_{(m_{1}-1)}(\theta+1)_{(n-m_{1}-1)}}$
$\displaystyle\qquad\qquad\qquad\times\frac{1}{\Gamma(1-\alpha)\Gamma(\theta+1)}\sum_{w=0}^{m}[\alpha(n+m-w-
m_{1})+\theta(w+m_{1})](w+m_{1})\binom{m}{w}$
$\displaystyle\qquad\qquad\qquad\qquad\times\Gamma(w+m_{1}-\alpha)\Gamma(\theta+n+m-w-
m_{1}).$

With some simple rearrangements of the terms the sum in the previous
expression can be written as

$\begin{split}&\sum_{w=0}^{m}[\alpha(n+m-w-m_{1})+\theta(w+m_{1})](w+m_{1})\\\
&\qquad\qquad\times\binom{m}{w}\Gamma(w+m_{1}-\alpha)\Gamma(\theta+n+m-w-
m_{1})\\\ &\qquad\qquad\qquad=[\alpha(n+m-m_{1})+\theta m_{1}]m_{1}R_{1}\\\
&\qquad\qquad\qquad\qquad+[m_{1}(\theta-\alpha)+\alpha(n+m-m_{1})+\theta
m_{1}]R_{2}+(\theta-\alpha)R_{3}\end{split}$ (36)

where we have defined

$\displaystyle R_{1}$
$\displaystyle:=\sum_{w=0}^{m}\binom{m}{w}\Gamma(w+m_{1}-\alpha)\Gamma(\theta+n+m-w-
m_{1})$ (37) $\displaystyle R_{2}$
$\displaystyle:=\sum_{w=0}^{m}w\binom{m}{w}\Gamma(w+m_{1}-\alpha)\Gamma(\theta+n+m-w-
m_{1})$ (38) $\displaystyle R_{3}$
$\displaystyle:=\sum_{w=0}^{m}w^{2}\binom{m}{w}\Gamma(w+m_{1}-\alpha)\Gamma(\theta+n+m-w-
m_{1}).$ (39)

We focus on the evaluation of $R_{1}$:

$\displaystyle R_{1}$
$\displaystyle=\Gamma(\theta+n+m-\alpha)\sum_{w=0}^{m}\binom{m}{w}B(w+m_{1}-\alpha,\theta+n+m-m_{1}-w)$
$\displaystyle=\Gamma(\theta+n+m-\alpha)\int_{0}^{1}\sum_{w=0}^{m}\binom{m}{w}x^{w+m_{1}-\alpha-1}(1-x)^{\theta+n+m-m_{1}-w-1}\mathrm{d}x$
$\displaystyle=\Gamma(\theta+n+m-\alpha)\int_{0}^{1}x^{m_{1}-\alpha-1}(1-x)^{\theta+n+m-m_{1}-1}\Big{(}1+\frac{x}{1-x}\Big{)}\mathrm{d}x$
$\displaystyle=\Gamma(\theta+n+m-\alpha)\int_{0}^{1}x^{m_{1}-\alpha-1}(1-x)^{\theta+n-m_{1}-1}\mathrm{d}x$
$\displaystyle=\Gamma(\theta+n+m-\alpha)B(m_{1}-\alpha,\theta+n-m_{1}),$

where we used the integral representation of the Beta function. The two sums
(38)–(39) may be evaluated in a similar fashion to get

$R_{2}=m\Gamma(\theta+n+m-\alpha)B(m_{1}+1-\alpha,\theta+n-m_{1})$

and

$\begin{split}R_{3}&=m(m-1)\Gamma(\theta+n+m-\alpha)B(m_{1}-\alpha+2,\theta+n-m_{1})\\\
&\qquad\qquad+m\Gamma(\theta+n+m-\alpha)B(m_{1}-\alpha+1,\theta+n-m_{1}).\end{split}$

Substituting the previous formulas for $R_{1},R_{2}$ and $R_{3}$ in (36), the
resulting expression can be used to evaluate the posterior expected value of
$W_{1,n+m}\mathds{1}_{B_{1}}$ to get

$\begin{split}&\mathds{E}[W_{1,n+m}\mathds{1}_{B_{1}}|K_{n}=k,\\{M_{i,n}=m_{i}\\}_{i=1}^{k}]\\\
&\qquad=\frac{n}{n+m}\cdot\frac{(\theta+1)_{(n-1)}}{[\alpha(n-m_{1})+\theta
m_{1}](\theta+1)_{(n+m-1)}}\cdot(\theta+n-\alpha)_{(m)}C(n,m,\alpha,\theta,m_{1})\end{split}$
(40)

where

$\begin{split}&C(n,m,\alpha,\theta,m_{1}):=m_{1}[\alpha(n+m-m_{1})+\theta
m_{1}]+m\frac{m_{1}-\alpha}{\theta+n-\alpha}[\alpha(n+m-m_{1})+\theta
m_{1}]\\\
&\qquad+\frac{m_{1}-\alpha}{\theta+n-\alpha}m(\theta-\alpha)(m_{1}+1)+\frac{(m_{1}-\alpha)(m_{1}+1-\alpha)}{(\theta+n-\alpha)(\theta+n+1-\alpha)}(\theta-\alpha)m(m-1),\end{split}$

and the thesis now follows putting together (35) and (40) to evaluate (22).
∎

###  A.7 Auxiliary results

We prove the following corollary of Theorem 1, which provides us with the
posterior probabilities of the events $A_{1}$ and $B_{1}$, defined in the main
paper.

######  Corollary 2.

Let $((T_{1},X_{1}),\ldots,(T_{n},X_{n}))$ be a random sample under the BNP
model (12), such that the sample features $K_{n}=k$ distinct species with
corresponding ordered frequencies $\mathbf{M}_{n}=\mathbf{m}$. Then,

$\text{Pr}(A_{1}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m})=1-\frac{n}{n+m}\cdot\frac{(\theta+n+1-\alpha)_{(m)}}{(\theta+n)_{(m)}}$
(41)

and

$\text{Pr}(B_{1}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m})=\frac{n}{n+m}\cdot\frac{(\theta+n+1-\alpha)_{(m)}}{(\theta+n)_{(m)}}.$
(42)

###### Proof.

The proof of Equation (19) and Equation (20) follows by a direct application
of Theorem 1. We now focus on the evaluation of (21) (42), and obviously (41)
follows from the fact that

$\text{Pr}[A_{1}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]=1-\text{Pr}[B_{1}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}].$

In order to evaluate the probability on the r.h.s. of the previous expression,
we consider (20) and we sum all over the possible values of the random
variable $W_{1,n+m}$:

$\displaystyle\text{Pr}[B_{1}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle\qquad=\sum_{w=0}^{m}\binom{m}{w}\cdot\frac{n}{n+m}\cdot\frac{\alpha(n+m-w-
m_{1})+\theta(w+m_{1})}{\alpha(n-m_{1})+\theta
m_{1}}\cdot\frac{(\theta+n-m_{1})_{(m-w)}(m_{1}-\alpha)_{(w)}}{(\theta+n)_{(m)}}$
$\displaystyle\qquad=\frac{n}{(n+m)(\theta+n)_{(m)}[\alpha(n-m_{1})+\theta
m_{1}]\Gamma(m_{1}-\alpha)\Gamma(\theta+n-m_{1})}(V_{1}+V_{2}),$

where $V_{1}$ and $V_{2}$ are defined as follows

$\displaystyle V_{1}$
$\displaystyle:=\sum_{w=0}^{m}\binom{m}{w}[\alpha(n+m)+m_{1}(\theta-\alpha)]\Gamma(m_{1}+w-\alpha)\Gamma(\theta+n+m-m_{1}-w)$
(43) $\displaystyle V_{2}$
$\displaystyle:=\sum_{w=0}^{m}\binom{m}{w}w(\theta-\alpha)\Gamma(m_{1}+w-\alpha)\Gamma(\theta+n+m-m_{1}-w).$
(44)

We focus on the evaluation of the two sums (43)–(44). As for $V_{1}$ we get

$\begin{split}V_{1}&=[\alpha(n+m)+m_{1}(\theta-\alpha)]\sum_{w=0}^{m}\binom{m}{w}\Gamma(m_{1}+w-\alpha)\Gamma(\theta+n+m-m_{1}-w)\\\
&=[\alpha(n+m)+m_{1}(\theta-\alpha)]\Gamma(\theta+n+m-\alpha)B(m_{1}-\alpha,\theta+n-m_{1})\end{split}$
(45)

where we have observed that the last sum coincides with $R_{1}$, that has been
evaluated in Section A.6 of the Supplementary Material. The sum (44) equals

$\begin{split}V_{2}&=(\theta-\alpha)\sum_{w=0}^{m}\binom{m}{w}w\Gamma(m_{1}+w-\alpha)\Gamma(\theta+n+m-m_{1}-w)\\\
&=(\theta-\alpha)m\Gamma(\theta+n+m-\alpha)B(m_{1}+1-\alpha,\theta+n-m_{1})\end{split}$
(46)

where we have now used the fact that the sum coincides with $R_{2}$, which has
been determined in Section A.6 of the Supplementary Material. We now
substitute (45)–(46) in the posterior probability of $B_{1}$, and with some
simple calculations we get

$\displaystyle\text{Pr}[B_{1}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle=\frac{n\Gamma(\theta+n+m-\alpha)B(m_{1}-\alpha,\theta+n-m_{1})}{(n+m)(\theta+n)_{m}[\alpha(n-m_{1})+\theta
m_{1}]\Gamma(m_{1}-\alpha)\Gamma(\theta+n-m_{1})}$
$\displaystyle\qquad\times\Big{\\{}\alpha(n+m)+(\theta-\alpha)\Big{[}m_{1}+m\frac{m_{1}-\alpha}{\theta+n-\alpha}\Big{]}\Big{\\}}.$

The last expression in parenthesis equals

$\alpha(n+m)+(\theta-\alpha)\Big{[}m_{1}+m\frac{m_{1}-\alpha}{\theta+n-\alpha}\Big{]}=(\theta+n+m-\alpha)\cdot\frac{\alpha(n-m_{1})+\theta
m_{1}}{\theta+n-\alpha},$

as a consequence the posterior probability of $B_{1}$ boils down to

$\displaystyle\text{Pr}[B_{1}|K_{n}=k,\mathbf{M}_{n}=\mathbf{m}]$
$\displaystyle=\frac{n\Gamma(\theta+n+m-\alpha)B(m_{1}-\alpha,\theta+n-m_{1})}{(n+m)(\theta+n)_{(m)}\Gamma(m_{1}-\alpha)\Gamma(\theta+n-m_{1})}\cdot\frac{\theta+n+m-\alpha}{\theta+n-\alpha}$
$\displaystyle=\frac{n}{n+m}\cdot\frac{(\theta+n+1-\alpha)_{(m)}}{(\theta+n)_{(m)}}$

where the last expression follows by a simple rearrangement of the terms, thus
(42) is now proved. ∎

##  Appendix B Additional empirical analyses

###  B.1 Synthetic data evaluation: additional results

In table 1 we report the performance of each parameter-estimation method, for
each of the four quantities considered: the number of distinct species $K$,
the frequency of the first ordered species $W_{1}$ and the frequency $W_{1}$
conditional on the events that the first ordered species is new ($A_{1}$) and
that the first ordered species is old ($B_{1}$). We report the median
percentage error, computed across 500 synthetic datasets, generated with
different distributions: the ordered PYP model (model), or with a clustering
distribution equal to either DP,PYP or zipf, and with a ordering distribution
equal to either alpha-stable or arrival-weighted. For each one of 100 initial
datasets of size $n=500$, we consider the median percentage error across the
25 additional datasets of size $m=5000$; we then report the median percentage
error across those 100 median values.

| alpha-stable | arrival-weighted |
---|---|---|---
method | DP | PYP | zipf | DP | PYP | zipf | model
K
stdPYP | 0.0007494 | 0.0143536 | 0.1244203 | 0.0007494 | 0.0143536 | 0.3653415 | 0.0163301
ordPYP | 0.3836682 | 0.3354845 | 0.1791527 | 0.0606133 | 0.0628106 | 0.5708139 | 0.0471591
ordDP | 0.2987086 | 0.7088840 | 0.7726780 | 0.2987086 | 0.7088840 | 0.9920688 | 0.7789935
lsX1 | 6.0248711 | 2.5801374 | 3.2309563 | 1.2033394 | 0.8679249 | 1.8934422 | 1.6933744
lsK | 0.0404877 | 0.0196868 | 0.0212205 | 0.0404877 | 0.0196868 | 0.3349258 | 0.0252152
W1
stdPYP | 0.6753089 | 2.1948600 | 0.9967755 | 7.6150681 | 9.7485602 | 2.5863652 | 1.1136050
ordPYP | 0.6172892 | 1.2689868 | 0.9969626 | 4.0246539 | 11.9130053 | 3.3076200 | 1.2318097
ordDP | 0.9622991 | 11.4874120 | 9.1622018 | 7.8403247 | 26.6549590 | 18.9922619 | 8.6353945
lsX1 | 0.7083109 | 0.8559754 | 0.8728699 | 5.4151449 | 12.7892906 | 6.6113843 | 0.6504245
lsK | 0.6414169 | 2.3759295 | 0.9962197 | 5.4395825 | 8.8118656 | 2.8914327 | 1.0972410
W1—A1
stdPYP | 0.9495687 | 0.5108009 | 0.4469288 | 0.9495687 | 0.5108009 | 0.4277088 | 0.5124739
ordPYP | 0.5048399 | 0.5054584 | 0.4370104 | 0.5095248 | 0.4916819 | 0.4838472 | 0.4964052
ordDP | 0.9638395 | 0.8003344 | 0.7204361 | 0.9638395 | 0.8003344 | 0.6854973 | 0.8001080
lsX1 | 0.5697234 | 0.4839626 | 0.1351994 | 0.7825655 | 0.5455795 | 0.3683756 | 0.5265883
lsK | 0.9229834 | 0.5206863 | 0.4439866 | 0.9229834 | 0.5206863 | 0.4214174 | 0.5209053
W1—B1
stdPYP | 0.3211614 | 0.5847100 | 0.7981944 | 3.5811886 | 0.9649008 | 0.9812986 | 0.5267659
ordPYP | 0.5540375 | 0.6886318 | 0.7967054 | 1.6272628 | 0.9702555 | 0.9825783 | 0.4759782
ordDP | 0.3211806 | 0.4122675 | 2.3333596 | 3.7804857 | 2.0640371 | 8.4222704 | 0.2828726
lsX1 | 0.6295685 | 0.7910964 | 0.9438083 | 1.9720726 | 0.9984228 | 0.9988799 | 0.5147092
lsK | 0.3978468 | 0.6342513 | 0.7638855 | 2.0510490 | 0.9671770 | 0.9843115 | 0.5222185
Table 1: Median percentage error for simulations based on synthetic data.

###  B.2 Genetic data evaluation: additional results

We report here the plots showing the curve of $K$ and $W_{1}$ for the EDAR
gene and the corresponding predictions and uncertainty quantification. Figure
7, which mirrors Figure 6, reports the curves (as black lines) of $K$ and
$W_{1}$ for the training-testing split that represents the median error for a
given parameter-estimating method, together with the corresponding predicted
curves (depicted in red) and the confidence intervals (red bands). The
confidence bands are created using bootstrap, computing the quantiles of the
curve estimates across the training-testing splits.

Figure 7: Prediction for the curve for the gene EDAR, for $K$ (top panels) and
$W_{1}$ (bottom panels).

For the analysis of genetic variation in gene EDAR, Figure 7 shows that the
best performance in predicting the number of distinct variants is given by
stdPYP and ordPYP, while ordDP and lsX1 better predict the frequency of the
oldest variant.

###  B.3 Analysis of citation data

While ordered SSP were originally motivated by problems arising in population
genetics, many other contexts can give rise to problems species sampling
problems where species are associated to an order. One such context is the
study of citation networks. In particular, given a collection of scientific
articles, each citation reported in any of those articles can be seen as an
individual observation, which can be classified into species based on which
article is cited. Moreover, cited articles can be ordered based on their
publication date.

In the following, we analyze citation data, using the dataset collected and
cleaned by Ji and Jin (2016), which focuses on papers published in the four
top statistics journals from 2003 to 2012. The data is publicly available at
https://www.stat.uga.edu/directory/people/pengsheng-ji, and detailed citation
information (bibtex files, used to compute the ordering of the papers) where
shared by the authors. Note that only the citations from articles in the
dataset to articles in the dataset are considered.

##### Statistics citation data

The Statistics Journals citation data provides the citation network of 3248
articles, of which 1798 contains citations to other papers in the dataset, and
1555 are cited by papers in the datasets, for a total of 5722 citations.
Within each journals, we ordered the papers by considering the year, the
journal issue number, and the pages. To order the articles across journals we
used the months for each issue, and when two journals published an issue in
the same months, the order given by the pages was used. The order is
increasing from the oldest paper to the most recent one.

We divide the original dataset into a training set and a testing set, so that
the size of the former is one tenth of the original sample size. We repeat our
analysis for 100 different randomly sampled training and testing sets. We
focus our analysis on the prediction of the number of distinct species $K$ and
the frequency of the first ordered species $W_{1}$. We compare the performance
of the different parameter-estimation methods described in Sections 3.2.

Figure 8: Prediction performance across several training-testing sets for the
statistics citation data. Note that results in the right panel are displayed
using a log-transformed $y$ axis.

Figure 8 reports the percentage absolute error for the prediction of $K$ and
$W_{1}$ in the statistics citation data, computed across the several training-
testing sets. We note that for the prediction of the number of distinct
species $K$ the best performance is achieved by stdPYP with a median
percentage error of $15\%$, followed by the ordPYP with a median percentage
error of $30\%$. For the prediction of the frequency of the oldest cited
paper, the best performance is achieved by ordDP, with a median percentage
error of $1\%$. We also note that the percentage errors are quite smaller
compared to those observed in Section 4.3, and this is due to the fact that in
this dataset the event $A_{1}$, i.e. that the oldest citation is observed in
the additional sample, is much more likely to occur.

Figure 9: Prediction for the curve for the statistics citation data, for $K$
(top panels) and $W_{1}$ (bottom panels).

Similar behaviors are observed in Figure 9, which displays the predictions as
functions of the number of samples, for the curve of the number of distinct
species $K$ (top panels) and the frequency of the oldest species $W_{1}$
(bottom panels). For simplicity of visualization, we depict in black the
actual curve of $K$ (or respectively $W_{1}$) observed in the training-testing
split that most closely represent the median error achieved by each method.
The prediction curve (represented in red) is also the one corresponding to the
“median” training-testing split, while the red bands represent the $95\%$
confidence bands. We report the empirical confidence bands, computed using the
empirical quantiles of the curve estimates across the various training-testing
splits. From the top panels of Figure 9 it is clear that stdPYP provides a
better fit to the whole curve of $K$, while the bottom panels show that ordDP
has the best prediction for the curve of $W_{1}$. We also note that the curve
of $W_{1}$ is not necessarily monotone increasing, as shown in the second and
fourth bottom panels, because a new older cluster can be observed and reduce
the frequency $W_{1}$ to one.

Overall, the parameters estimated with standard Pitman Yor Process seems to
provide a better fit for the behavior of the number of distinct species
(number of unique cited papers), but the parameters estimated with the ordered
Dirichlet Process seems to better predict the frequency of the oldest cited
paper $W_{1}$. This contrasting behavior is most likely explained by the fact
that the ordered PYP model is somewhat misspecified for these data, which
display power-law behavior for the number of species (consistent with the
standard PYP), and a first ordered species frequency similar to that induced
by the ordered DP.