# A Gibbs Sampling Scheme for a Generalised Poisson-Kingman Class

Robert C. Griffiths  Ross A. Maller  Soudabeh Shemehsavar Email:
<EMAIL_ADDRESS>School of Mathematics, Murdoch University, Perth,
Western Australia; School of Mathematics, Statistics & Computer Sciences,
University of Tehran.
Email<EMAIL_ADDRESS>(corresponding author)

###### Abstract

A Bayesian nonparametric method of James, Lijoi & Prunster (2009) used to
predict future values of observations from normalized random measures with
independent increments is modified to a class of models based on negative
binomial processes for which the increments are not independent, but are
independent conditional on an underlying gamma variable. Like in James et al.,
the new algorithm is formulated in terms of two variables, one a function of
the past observations, and the other an updating by means of a new
observation. We outline an application of the procedure to population
genetics, for the construction of realisations of genealogical trees and
coalescents from samples of alleles.

Keywords:
2010 Mathematics Subject Classification: Primary 60G51, 60G52, 60G55.

##  1 Introduction

(James, Lijoi & Pr$\ddot{\rm u}$nster2009) provided a comprehensive Bayesian
nonparametric analysis of random probabilities, considered as generalizations
of the Dirichlet process, obtained by normalizing random measures with
independent increments (NRMIs). They formulated the posterior distribution of
an NRMI as a mixture with respect to the distribution of a specific latent
variable, showed that the posterior distribution is again an NRMI, and derived
from the ensuing representation a rule for predicting future values of
observations, given those already observed. This takes the form of a mixture
of two variables, one a function of the past observations, and the other an
updating by means of the new observation. From this they developed a Gibbs
sampling algorithm generalising the Blackwell-MacQueen urn scheme which
provides a mechanism for simulating sample observations from the NRMI,
producing explicit and tractable expressions suitable for practical
implementation.

This procedure is closely related to algorithms due to Griffiths & Tavaré
(1994) for computing probability distributions of sample configurations of
alleles in population genetics. Griffiths & Tavaré (1994) derive a recursion
formula satisfied by the sampling probabilities and use it to construct a
Markov chain framework for simulating the likelihood of the data, given a set
of parameter values. This approach corresponds to an importance sampling
procedure (Felsensteinetal (1999), Stephens & Donnelly (2000)). The recursive
nature of the method results from looking back in time in the coalescent tree
to the most recent event in the past (either a coalescence or a mutation) to
see what type of tree configuration there must have been at that time to
produce the current tree configuration.

In the present paper we modify the (James, Lijoi & Pr$\ddot{\rm u}$nster2009)
method (also referring to (James2002, James2005)) to handle a wide class of
models based on negative binomial processes for which the increments are not
independent, but are independent conditional on an underlying gamma variable.
The formulae we derive allow adaptation of the Griffiths & Tavaré (1994)
methods for the construction of realisations of trees and coalescents from
observations, either simulated or actual, on the more general processes.

###  1.1 The setup

Let $(\xi_{i})_{i\in\mathbb{N}}$, $\mathbb{N}=\\{1,2,\ldots\\}$, be a sequence
of i.i.d. random variables on $\mathbb{R}$ with a specified initial “base”
distribution $G_{0}$. Let $(P_{i})_{i\in\mathbb{N}}$ with
$P_{1}>P_{2}>\cdots>0$ be a random weight sequence on the infinite simplex
$\nabla_{\infty}=\\{z_{i}\in[0,1]:z_{1}\geq z_{2}\geq\cdots\geq
0,\,\sum_{i=1}^{\infty}z_{i}=1\\}$, and $G$ a random discrete distribution
generated from $(P_{i})_{i\in\mathbb{N}}$ and $(\xi_{i})_{i\in\mathbb{N}}$.
With $\delta_{x}$ denoting a point mass at $x\in\mathbb{R}$, $G$ can be
written as

$G=\sum_{i\geq 1}P_{i}\delta_{\xi_{i}}.$ (1.1)

The data sample ${\bf X}=\\{X_{1},\ldots,X_{n}\\}$ is obtained by drawing $n$
values from $G$. Thus, conditionally given $G$, the $X_{i}$ are i.i.d. with
distribution $G$. In the context of population genetics, the $X_{i}$ may
represent the allele types of a particular locus on the chromosome.

Due to the assumed discreteness of $G$, the observed data ${\bf X}$ in general
contains ties, and can be reduced to a vector of unique values, ${\bf
Y}=\\{Y_{1},\ldots,Y_{K_{n}}\\}$, and their corresponding multiplicities,
${\bf N}_{[n]}=(N_{[n]}^{1},N_{[n]}^{2},\ldots,N_{[n]}^{Kn})$. Here $K_{n}$ is
a random variable increasing in $n$ representing the number of unique values –
clusters, or blocks – in ${\bf X}$. Given ${\bf X}$, hence $K_{n}$, $K_{n}=k$,
say, ${\bf N}_{[n]}$ takes values among $k$-vectors of positive integers ${\bf
n}=\\{n_{1},\ldots,n_{k}\\}$ with $\sum_{i=1}^{k}n_{i}=n$. We consider the
homogenous case where the $X_{i}$ are independent of the $P_{i}$.

###  1.2 Distributions on the infinite simplex: $\textbf{\rm PK}(\rho)$ and
$\textbf{\rm PK}^{(r)}(\rho)$

In the Poisson case, the $P_{i}$ in (1.1) are generated from a Poisson point
process $\mathbb{B}$ whose intensity measure has density $\rho(s)$ satisfying

$\lim_{s\downarrow 0}\rho(s)=\infty,\quad\rho(s)<\infty\text{ for all
}s>0,\quad\text{ and }\int_{0}^{1}s\rho(s){\rm d}s<\infty.$ (1.2)

We can write $\mathbb{B}$ in the form $\mathbb{B}=\sum_{i\geq
1}\delta_{S_{i}}$ with points $S_{1}>S_{2}>\cdots$. The $P_{i}$ for the
process $\mathbb{B}$ are defined as the ratios of the ranked points of
$\mathbb{B}$ normalised by their sum:

$(P_{i})_{i\in\mathbb{N}}=\Big{(}\frac{S_{i}}{\sum_{\ell\in\mathbb{N}}S_{\ell}}\Big{)}_{i\in\mathbb{N}}.$
(1.3)

This kind of construction was first given by Kingman (1975). In his case the
$S_{i}$ were the ranked jumps $\Delta S_{1}^{(1)}>\Delta S_{1}^{(2)}>\cdots$
of a subordinator $(S_{t})_{t\geq 0}$ with Lévy density $\rho(s)$ taken at
time $1$.

The distributions of the $(P_{i})_{i\in\mathbb{N}}$ thus constructed are said
to be in $\textbf{\rm PK}(\rho)$, the Poisson-Kingman class. These processes
have been intensively studied for certain specifications of $\rho$. We are
interested in a generalised version which we denote as $\textbf{\rm
PK}^{(r)}(\rho)$. To construct it, we modify (1.3) as follows.

######  Definition 1.1.

[Definition of $\mathcal{BN}(r,\rho)$] Take a $\rho$ satisfying (1.2) and a
parameter $r>0$. Following Gregoire (1984), introduce a negative binomial
process, $\mathbb{B}^{(r)}$, defined in terms of its Laplace functional:

$\mathbb{E}\big{(}e^{\mathbb{B}^{(r)}(f)}\big{)}=\Big{(}1+\int_{\mathbb{R}_{+}}(1-e^{-f(x)})\rho(x){\rm
d}x\Big{)}^{-r},$ (1.4)

for any bounded nonnegative measurable function $f$ on $\mathbb{R}_{+}$. Such
processes are said to be in $\mathcal{BN}(r,\rho)$.

Basic properties of a negative binomial process are that numbers of points in
non-overlapping Borel sets are negative binomially (but not independently)
distributed.

(1.4) shows that we can regard a process in $\mathcal{BN}(r,\rho)$ to be a
Poisson point process with random intensity measure $\Gamma_{r}\rho$, where
$\Gamma_{r}:=\Gamma_{r,1}$, a Gamma random variable111We use the notation
$\Gamma_{\beta,\gamma}$ for a Gamma random variable with parameter
$(\beta,\gamma)$ having density $(\gamma x)e^{-\gamma x}(\gamma
x)^{\beta-1}/\Gamma(\beta)$, $\beta,\gamma>0$, where $\Gamma(\cdot)$ is the
gamma function.  with parameter $(r,1)$. Using such a $\mathbb{B}^{(r)}$, we
construct $\textbf{\rm PK}^{(r)}(\rho)$ as follows.

######  Definition 1.2.

[Definition of $\textbf{\rm PK}^{(r)}(\rho)$] Write
$\mathbb{B}^{(r)}=\sum_{i\in\mathbb{N}}J_{i}$ in terms of its ranked points
$J_{1}>J_{2}>J_{3}>\cdots$. Analogous to (1.3), construct a random vector from
the ranked points in $\mathbb{B}^{(r)}$ normalised by their sum. The resulting
random sequence on the infinite simplex is said to have distribution
$\textbf{\rm PK}^{(r)}(\rho)$; thus,

$\Big{(}\frac{J_{i}}{\sum_{\ell\in\mathbb{N}}{J_{\ell}}}\Big{)}_{i\in\mathbb{N}}\sim\textbf{\rm
PK}^{(r)}(\rho),\quad r>0.$ (1.5)

Lévy densities of particular interest are

$\rho_{\alpha}(x)=\frac{\alpha}{\Gamma(1-\alpha)}x^{-\alpha-1},\quad\rho_{\theta}(x)=\theta
x^{-1}e^{-x},\quad\rho_{G}(x)=\frac{\alpha}{\Gamma(1-\alpha)}x^{-\alpha-1}e^{-x},x>0,$
(1.6)

where $\theta>0$ and $0<\alpha<1$. We can identify important members of the PK
class in the $\textbf{\rm PK}^{(r)}$ framework. For example, the two-parameter
Poisson-Dirichlet process $\textbf{\rm PD}(\alpha,\theta)$ (Pitman & Yor
(1997)) is obtained as $\textbf{\rm PD}(\alpha,\theta)=\textbf{\rm
PK}^{(\theta/\alpha)}(\rho_{G})$. Further, $\textbf{\rm
PD}(\alpha,0)=\textbf{\rm PK}^{(r)}(\rho_{\alpha})$ for any $r>0$, and
$\textbf{\rm PD}(0,\theta)$ can be derived as a limiting case: $\textbf{\rm
PD}(0,\theta)=\lim_{\alpha\to 0}\textbf{\rm
PK}^{(\theta/\alpha)}(\rho_{G})=\lim_{r\to\infty}\textbf{\rm
PK}^{(r)}(\rho_{\theta}/r)$ (Pitman & Yor (1997), Ipsen & Maller (2017)).

The choice $\rho_{R}(x)=\alpha x^{-\alpha-1}{\bf 1}_{0<x\leq 1}$ in (1.2)
gives the class $\textbf{\rm PD}_{\alpha}^{(r)}$ studied in Ipsen, Maller &
Shemehsavar (2020a, 2021) and Maller & Shemehsavar (2023b) and related to the
negative binomial process in (1.4). Chegini & Zarepour (2023) and Labadi &
Zarepour (2014) give further properties and applications of $\textbf{\rm
PD}_{\alpha}^{(r)}$. The distribution of $\textbf{\rm PD}_{\alpha}^{(r)}$
converges to that of $\textbf{\rm PD}(\alpha,0)$ as $r\downarrow 0$ (Maller &
Shemehsavar (2023b)).

##  2 Posterior Distributions for $\textbf{\rm PK}^{(r)}(\rho)$

(James2002, James2005) developed a framework to compute the likelihood and
posterior distributions of partitioned data derived from an underlying Poisson
point process. In particular, (James, Lijoi & Pr$\ddot{\rm u}$nster2009)
computed the posterior analysis for $G$ when $(P_{i},i\geq 1)$ are in the
Poisson-Kingman class. We use the same framework to derive the posterior
distribution of $G$ when the weight sequence is distributed as $\textbf{\rm
PK}^{(r)}(\rho)$. Like James et al. we make use of certain auxiliary random
variables, denoted $(V_{i})_{i\in\mathbb{N}}$ in our notation, to be specified
later.

By analogy with their approach, we define the following negative binomial
point process augmented by the inclusion of an independent sequence $\xi_{i}$:

$\mathbb{M}^{(r)}(\cdot,\cdot):=\sum_{i}\delta_{(J_{i},\xi_{i})}\sim\mathcal{BN}(r,\rho,G_{0}),$
(2.1)

where $(J_{i})_{i\in\mathbb{N}}$ are the ranked points of a
$\mathcal{BN}(r,\rho)$ process $\mathbb{B}^{(r)}$, and
$(\xi_{i})_{i\in\mathbb{N}}$ are i.i.d. with a distribution $G_{0}$,
independent of $(J_{i})_{i\in\mathbb{N}}$. Note that $\mathbb{M}^{(r)}$ is no
longer a Poisson random measure, but, conditional on a $\Gamma_{r}$ random
variable, $\mathbb{M}^{(r)}$ is a Poisson point process with intensity measure
$\Gamma_{r}\rho(\cdot)G_{0}(\cdot)$. We can relate $\mathbb{M}^{(r)}$ to the
species sampling model $G$ in (1.1) when $(P_{i})_{i\in\mathbb{N}}$ is
distributed as $\textbf{\rm PK}^{(r)}(\rho)$ by setting

$G(\cdot)=\frac{\mu(\cdot)}{T^{(r)}},\ \text{where}\
\mu(\cdot)=\int_{s>0}s\,\mathbb{M}^{(r)}({\rm d}s,\cdot)\ \text{and}\
T^{(r)}=\mu(\Omega)=\sum_{i\in\mathbb{N}}J_{i}.$ (2.2)

We aim to use the methods of (James, Lijoi & Pr$\ddot{\rm u}$nster2009) to
derive the posterior distribution of the species sampling model $G$ in (2.2)
given observations ${\bf X}$. In the analysis, the $n$ observations in ${\bf
X}$ appear in terms of their $k$ unique values ${\bf
Y}=\\{Y_{1},\ldots,Y_{k}\\}$ and their respective multiplicities ${\bf
N}_{[n]}$.

Our first task is to compute the posterior distribution of the underlying
augmented measure $\mathbb{M}^{(r)}$ in (2.1), given ${\bf X}$. In the
following analysis, we fix $r>0$ and base distribution $G_{0}$. Define

$\psi(v)=1+\int_{x>0}(1-e^{-vx})\rho(x){\rm d}x\ \text{and}\
\pi_{n}(v)=(-1)^{n-1}\psi^{(n)}(v)=\int_{x>0}x^{n}\rho(x)e^{-vx}{\rm d}x,\
v>0.$ (2.3)

Conditional on ${\bf X}$, hence with ${\bf Y}$, $k$ and ${\bf
n}=\\{n_{1},\ldots,n_{k}\\}$ given, define an auxiliary variable
$V_{n}:=V_{n}(r)$ with conditional density proportional to

$g_{r}(v,{\bf
n}):=\frac{r^{[k]}}{(\psi(v))^{r+k}}\cdot\,\frac{v^{n-1}}{\Gamma(n)}\prod_{i=1}^{k}\pi_{n_{i}}(v),\
v>0,$ (2.4)

where $r^{[k]}=\Gamma(r+k)/\Gamma(r)$, $r>0$. Our main results stem from:

######  Theorem 2.1.

Conditional on $V_{n}$ and ${\bf X}$, we have the following decomposition:

$\mathbb{M}_{\bf
X}^{(r,V_{n})}=\widetilde{\mathbb{M}}^{(r+k,V_{n})}+\sum_{i=1}^{k}\delta_{J_{i}^{V_{n}},Y_{i}},$
(2.5)

where, for each $v>0$,

$\widetilde{\mathbb{M}}^{(r+k,v)}\sim\mathcal{BN}\Big{(}r+k,\frac{e^{-vs}}{\psi(v)}\rho(s),G_{0}\Big{)},\
{\rm independent\ of}\ \mathbb{M}^{(r)},$ (2.6)

and the $J_{i}^{v}$, $1\leq i\leq k$, are independent with densities

$P\big{(}J_{i}^{v}\in{\rm d}s\big{)}=\frac{s^{n_{i}}e^{-vs}\rho(s){\rm
d}s}{\pi_{n_{i}}(v)},\ s>0.$ (2.7)

######  Remark 2.1.

The process $\widetilde{\mathbb{M}}$ in (2.6) is constructed from a negative
binomial point process whereas the analogous process in (James, Lijoi &
Pr$\ddot{\rm u}$nster2009) is from a Poisson point process. A significant
distinction is the appearance of the denominator term $\psi(v)$ in (2.6) which
of course changes the dynamics of the process appreciably.

The introduction of the auxiliary random variable $V_{n}$ conveniently
decomposes the posterior structure of $\mathbb{M}^{(r)}$ into an infinite
measure part and a finite dimensional part. Conditional on $V_{n}$, the
posterior distribution of a negative binomial process is again a negative
binomial point process $\mathcal{BN}(r+k,\rho)$, noting in particular the
change in parameter from $r$ to $r+k$, where $k$ is the current number of
distinct values in ${\bf X}$.

Using the posterior distribution for the underlying negative binomial process
in Theorem 2.1, it is immediate to give the posterior distributions of $\mu$
and $G$ in (2.2). The following modifies Theorems 1 and 2 of (James, Lijoi &
Pr$\ddot{\rm u}$nster2009):

######  Theorem 2.2.

(i) Conditional on ${\bf X}$ and $V_{n}$, the posterior distribution of $\mu$
is

$\mu_{\bf X}^{V_{n}}=\mu^{V_{n}}+\sum_{i=1}^{k}J_{i}^{V_{n}}\delta_{Y_{i}}$
(2.8)

where $\mu^{V_{n}}({\rm d}x)=\int_{s>0}s\mathbb{M}^{(r+k,V_{n})}({\rm d}s,{\rm
d}x)$ and $(J_{i}^{V_{n}},{Y_{i}})_{i=1,\ldots,k}$ are defined as in Theorem
2.1.

(ii) Conditional on ${\bf X}$ and $V_{n}$, the posterior distribution of $G$
is, with $T^{V_{n}}=\mu^{V_{n}}(\Omega)$,

$G_{\bf X}^{V_{n}}({\rm d}x)=\frac{\mu^{V_{n}}({\rm
d}x)}{T^{V_{n}}+\sum_{i=1}^{k}J_{i}^{V_{n}}}+\sum_{i=1}^{k}\frac{J_{i}^{V_{n}}}{T^{V_{n}}+\sum_{i=1}^{k}J_{i}^{V_{n}}}\delta_{Y_{i}}({\rm
d}x).$ (2.9)

We obtain the marginal distribution of ${\bf
N}_{[n]}=(N_{[n]}^{1},N_{[n]}^{2},\ldots,N_{[n]}^{K_{n}})$ by integrating out
$V_{n}$ in (2.4). This produces the _exchangeable partition probability
function (EPPF)_ :

$p({\bf
n})=p(n_{1},n_{2},\ldots,n_{k})=\int_{v}\Big{[}\prod_{i=1}^{k}\pi_{n_{i}}(v)\Big{]}\times\frac{v^{n-1}}{\Gamma(n)}r^{[k]}\times(\psi(v))^{-(r+k)}{\rm
d}v=\int g_{r}(v,{\bf n}){\rm d}v.$ (2.10)

##  3 Prediction Rules

Here we state our main result, a formula for the predictive distribution of
$X_{n+1}$ given the data ${\bf X}$. As we show in Section 4, this can be used
to construct an urn scheme for sampling analogous to the Chinese restaurant
process commonly used to visualise these kinds of procedure. The following
formula (3.1) is exactly analogous to formula (10) in Proposition 2 of (James,
Lijoi & Pr$\ddot{\rm u}$nster2009), but again note the appearance of the
denominator term $\psi(v)$ in (3.2) as a major distinction. For proof of
Theorem 3.1 see the Appendix.

######  Theorem 3.1.

Let $G$ be a species sampling model constructed from a negative binomial
process with parameter $r$, intensity measure $\rho$ and base distribution
$G_{0}$ as in (2.2). Then the predictive distribution for $X_{n+1}$ given
${\bf X}$ has the same law as

$\mathbb{P}(X_{n+1}\in{\rm d}x|{\bf X})=\omega_{0}^{(n)}G_{0}({\rm
d}x)+\frac{1}{n}\sum_{i=1}^{k}\omega_{i}^{(n)}\delta_{Y_{i}}({\rm d}x)$ (3.1)

where

$\omega_{0}^{(n)}=\frac{r+k}{n}\int_{0}^{\infty}v\frac{\pi_{1}(v)}{\psi(v)}\,g_{r}(v,{\bf
n}){\rm d}v,$ (3.2)

and for $i=1,\ldots,k$,

$\omega_{i}^{(n)}=\int_{0}^{\infty}v\frac{\pi_{n_{i}+1}(v)}{\pi_{n_{i}}(v)}\,g_{r}(v,{\bf
n}){\rm d}v.$ (3.3)

######  Remark 3.1.

(i) The probability of forming a new cluster given ${\bf X}$ is
$\omega_{0}^{(n)}$, which can be written as

$\omega_{0}^{(n)}=\frac{r}{n}\int_{0}^{\infty}\frac{v\pi_{1}(v)(r+1)^{[k]}}{(\psi(v))^{r+1+k}}\frac{v^{n-1}}{\Gamma(n)}\prod_{i=1}^{k}\pi_{n_{i}}(v){\rm
d}v=\frac{r}{n}\int_{0}^{\infty}v\pi_{1}(v)g_{r+1}(v,{\bf n}){\rm d}v.$ (3.4)

(ii) The allele frequency spectrum (AFS) is the $n$-vector ${\bf
M}_{n}=(M_{1n},\ldots,M_{nn})$, where $M_{jn}=\sum_{i=1}^{k}{\bf
1}_{\\{N_{[n]}^{i}=j\\}}$, $1\leq j\leq n$. We have $\sum_{j=1}^{n}jM_{jn}=n$
and $\sum_{j=1}^{n}M_{jn}=k$. ${\bf M}_{n}$ takes values among $n$-vectors of
nonnegative integers ${\bf m}=(m_{1},\ldots,m_{n})$ satisfying

$A_{nk}:=\Big{\\{}{\bf m}=(m_{1},\ldots,m_{n}):m_{j}\geq
0,\,\sum_{j=1}^{n}jm_{j}=n,\,\sum_{j=1}^{n}m_{j}=k\Big{\\}}.$ (3.5)

(iii) (Posterior Distribution of $K_{n}$) Let $(D_{i},i=1,\ldots,n)$ be
independent Bernoulli($\omega_{0}^{(i)}$) random variables such that $D_{i}=1$
when the $i$th individual sampled takes a different value from the previous
$i$ labels. In the Chinese Restaurant process analogy, the $i$th customer sits
at a different table with probability $\omega_{0}^{(i)}$. The posterior
distribution of $K_{n}$ is $K_{n}^{\bf X}=1+\sum_{i=1}^{n}D_{i}$.

From (2.4) we can write $g_{r}(v,{\bf n})$ alternatively in terms of the
$m_{j}$ as:

$g_{r}(v,{\bf
m}):=\frac{r^{[k]}}{(\psi(v))^{r+k}}\cdot\,\frac{v^{n-1}}{\Gamma(n)}\prod_{j=1}^{n}(\pi_{j}(v))^{m_{j}},\
v>0.$ (3.6)

The next proposition as a check verifies that the probabilities in (3.1) and
(3.2) add to the EPPF, and a renormalisation of the EPPF gives the
corresponding partition distribution.

######  Proposition 3.1.

Using the notation in Theorem 3.1, we have

$\omega_{0}^{(n)}+\frac{1}{n}\sum_{i=1}^{k}\omega_{i}^{(n)}=\int_{0}^{\infty}g_{r}(v,{\bf
n}){\rm d}v;$ (3.7)

and

$\sum_{k=1}^{n}\sum_{{\bf m}\in
A_{nk}}\frac{n!}{\prod_{j=1}^{n}j!^{m_{j}}m_{j}!}\int_{0}^{\infty}g_{r}(v,{\bf
m}){\rm d}v=1.$ (3.8)

##  4 A Gibbs sampling scheme for $\textbf{\rm PD}_{\alpha}^{(r)}$

The following Gibbs sampling scheme generalises the Blackwell-MacQueen urn
scheme in (James, Lijoi & Pr$\ddot{\rm u}$nster2009), Section 3.3. Choose a
value of ${\bf X}=(X_{1},\ldots,X_{n})$ as follows.

  1. (i)

sample $V_{0}$ from density (2.4) with $n=k=1$, which is proportional to
$\frac{r\pi_{1}(v)}{\psi(v)}$;

  2. (ii)

sample $X_{1}$ from $\mathbb{P}({\rm d}X_{1}|V_{0})$, which is proportional to
$\frac{r\pi_{1}(V_{0})}{\psi(V_{0})}G_{0}({\rm d}x)$;

  3. (iii)

for $\ell\geq 1$,

    1. (a)

sample $V_{\ell}$ from $\mathbb{P}({\rm d}V_{\ell}|X_{1},\ldots,X_{\ell})$,
which is proportional to (2.4) with $k$ replaced by $\ell$;

    2. (b)

sample $X_{\ell+1}$ from $\mathbb{P}({\rm
d}X_{\ell+1}|X_{1},\ldots,X_{\ell},V_{\ell})$ using the conditional prediction
rule $\mathbb{P}(X_{n+1}\in{\rm d}x|{\bf X},V_{n}=v)$ proportional to

$\displaystyle\frac{(r+\ell)\pi_{1}(v)}{\psi(v)}G_{0}({\rm
d}x)+\sum_{i=1}^{\ell}\frac{\pi_{n_{i}+1}(v)}{\pi_{n_{i}}(v)}\delta_{Y_{i}}({\rm
d}x);$ (4.1)

  4. (iv)

repeat (iii) until $X_{n}$ is reached.

For analysis, collect the $X_{\ell}$ into vector ${\bf
X}=(X_{1},\ldots,X_{n})$, and reduce ${\bf X}$ to a vector of unique values,
${\bf Y}=\\{Y_{1},\ldots,Y_{K_{n}}\\}$, and their corresponding
multiplicities, ${\bf n}=\\{n_{1},\ldots,n_{K_{n}}\\}$, with
$\sum_{i=1}^{K_{n}}n_{i}=n$, as in Subsection 1.1. Note that the EPPF formula
(2.10) recovers Eq. (5.6) in Ipsen, Maller & Shemehsavar (2020a) when
$\rho(x){\rm d}x=\alpha x^{-\alpha-1}{\rm d}x{\bf 1}_{0<x\leq 1}$, so the
algorithm in Subsection 4 does indeed produce samples from $\textbf{\rm
PD}_{\alpha}^{(r)}$.

##  5 Constructing Trees and Coalescents from a $\textbf{\rm
PD}_{\alpha}^{(r)}$ Sample

A genealogical tree is a diagram that depicts the ancestral history and
relationships between genes in a sample from within a species. The tree’s
branches represent lines of descent and the branches’ tips are labelled with
the allele types. The pattern of branching reflects relationships between
alleles; they are more closely related when there are fewer nodes on the path
from one to another through the tree. A fundamental objective of population
genetics analysis is to construct genealogical trees from DNA sequence data
with the aim of developing probability models to describe the ancestral
processes involved.

Good estimates of tree configurations are crucial for quantifying statistical
uncertainty and estimating parameters such as effective population size,
mutation rate, allele age, etc. See Hudson (1991)) for a general overview of
the applications. Griffiths & Tavaré (1994) describe a sequential importance
sampling algorithm for calculating the probability of a sample of size $n$
under mutation and genetic drift in an unstructured population. De Iorio &
Griffiths (2004) developed an improved method of constructing sequential
importance sampling proposal distributions on coalescent histories. In this
section we outline an extension of the Griffiths & Tavaré (1994) algorithm to
the $\textbf{\rm PK}^{(r)}$ class, deriving a backwards recursion, depending
on ancestral histories, for the sampling distribution. A time scale is
introduced for the backwards processes. Such kinds of recursion are important
in sequential importance sampling and duality with diffusion processes as
illustrated for the one and two parameter Poisson-Dirichlet processes in
Griffiths, Ruggiero, Spano, and Zhou (2023) and Petrov (2009).

Consider a population consisting of a number of genes, each of which has an
infinite number of nucleotide sites, such that at each site there are two or
more possible nucleotides, the original one and mutant ones. In the
infinitely-many-sites model of DNA sequences the “type” of a gene is given by
a sequence of ls and 0s indicating at which sites a mutation has occurred in
the gene or one of its ancestors. The reproductive mechanism assumes in effect
that each of the genes in a generation under consideration selects a parent
gene at random (with replacement) from the previous generation. With a
designated probability the offspring gene is of the same type as the parent
gene; with the complementary probability a mutation occurs, in which case a 0
changes to a 1 at a site where no mutation has occurred in any gene. Which
site is chosen for the mutation does not matter, as sites are unordered. At an
allele level sequences are recorded as being distinct if their DNA types are
different, but the fine details of site mutations are not taken into account.
See Ethier & Griffiths (1987) for more details.

In a reverse procedure, for a given set $D$ (descendents) of input sequences
in the generation under consideration, a genealogy for $D$ can be constructed
backward in time using the evolutionary events of mutation and coalescence.
Types are thought of as mutations and removing a singleton is like removing a
mutation. An ancestral configuration can be regarded as the multiset of all
sequences present at a particular point in time in a possible genealogy for
$D$. It may be thought of as a tree (reading forward in time) or a coalescent
(reading backward in time).

The random vector ${\bf N}_{[n]}$ of Section 2, which takes values among
partitions of $[n]$ is now interpreted as counting the numbers of distinct
alleles in a sample of size $n$ arranged in their order of appearance in the
sample. Suppose we currently have an observation ${\bf
n}=(n_{1},\ldots,n_{k})$ on ${\bf N}_{[n]}$, where $n_{i}$ is the number of
allele type $i$ among the $n$ alleles, of which there are $k$ distinct types.
Let ${\bf e}_{i}=(0,0,\ldots,0,1_{i},0,\ldots,0,0)$ be a unit $n$-vector with
a 1 in the $i$-th position. Recalling the notation in Theorem 3.1, we argue
that to obtain the configuration ${\bf n}$ there was a birth from ${\bf
n}-{\bf e}_{i}$ with probability

$\frac{n_{i}}{n}\frac{1}{n-1}\omega_{i}^{(n-1)}({\bf n}-{\bf e}_{i}),$ (5.1)

when $n_{i}>1$, or, for a singleton ($n_{i}=1$), with probability

$\frac{1}{n}\omega_{0}^{(n-1)}({\bf n}-{\bf e}_{i}).$ (5.2)

All the terms when $n_{i}=1$ are identical. The factor of $n_{i}/n$ appears
because each individual in the current sample (of size $n$) has an equal
probability of being the new individual from the immediate ancestral sample
(of size $n-1$).

A full investigation of the coalescent structure and biological relevance of
(5.1) and (5.2) for particular models is beyond the scope of this paper and
yet to be determined. But we can give some preliminary observations as
follows. From (5.1) and (5.2) we get a fundamental equation for the EPPF,
$\text{p}(\bm{n})$, in (2.10):

$\displaystyle\text{p}({\bf n})$
$\displaystyle=\sum_{n_{i}>1}\frac{n_{i}}{n}\cdot\frac{1}{n-1}\omega_{i}^{(n-1)}({\bf
n}-{\bf e}_{i})+\sum_{\\{i:n_{i}=1\\}}\frac{1}{n}\omega_{0}^{(n-1)}({\bf
n}-{\bf e}_{i}).$ (5.3)

The second term counts the number of singletons in the sample, multiplied by
the probability of a singleton, $\omega_{0}^{(n)}$. All the terms in the
second sum are identical.

Based on (5.3) we define a continuous time Markov process

$\displaystyle\frac{{\rm d}}{{\rm d}t}\text{p}({\bf n},t)$
$\displaystyle=-\varphi({\bf n})\text{p}({\bf n},t)$
$\displaystyle+\varphi({\bf
n})\sum_{n_{i}>1}\frac{n_{i}}{n(n-1)}\cdot\frac{\omega_{i}^{(n-1)}({\bf
n}-{\bf e}_{i})}{\text{p}({\bf n}-{\bf e}_{i})}\cdot\text{p}({\bf n}-{\bf
e}_{i},t)$ $\displaystyle+\varphi({\bf
n})\sum_{\\{i:n_{i}=1\\}}\frac{1}{n}\frac{\omega_{0}^{(n-1)}({\bf n}-{\bf
e}_{i})}{\text{p}({\bf n}-{\bf e}_{i})}\cdot\text{p}({\bf n}-{\bf e}_{i},t),$
(5.4)

where $\varphi({\bf n})$ is a total rate for the configuration ${\bf n}$.
Rescale $\text{p}({\bf n},t)$ as

$\displaystyle H({\bf n},t)$ $\displaystyle=\frac{\text{p}({\bf
n},t)}{\text{p}({\bf n})},\ t\geq 0,$

so that

$\frac{{\rm d}}{{\rm d}t}H({\bf n},t)=-\varphi({\bf n})H({\bf
n},t)+\varphi({\bf n})\sum_{n_{i}>1}\frac{n_{i}}{n}H({\bf n}-{\bf
e}_{i},t)+\varphi({\bf n})\sum_{\\{i:n_{i}=1\\}}\frac{1}{n}H({\bf n}-{\bf
e}_{i},t).$ (5.5)

The rates in (5.5) are $\varphi({\bf n})n_{i}/n$ when both $n_{i}>1$ or
$n_{i}=1$, and of course we have

$\sum_{i=1}^{k}\varphi({\bf n})\frac{n_{i}}{n}=\varphi({\bf n}).$

$H({\bf n},t)$ is a function of a process $({\bf n}(t),t\geq 0)$, which begins
at ${\bf n}(0)={\bf n}$ and decreases back in time by coalescing within a
group where $n_{i}>1$ or removing a singleton (“mutation”) when $n_{i}=1$. The
rates ${\bf n}\to{\bf n}-{\bf e}_{i}$ are $\varphi(\bm{n})n_{i}/n$. The choice
of $H({\bf n},0)$ is important. In the one and two parameter Poisson-Dirichlet
processes the construction of $H({\bf n},t)$ leads to a coalescent process
dual to the diffusions.

It is easy to convert p$({\bf n})$ into a partition distribution by taking

${\cal P}({\bf n})={n\choose{\bf n}}\frac{1}{m_{1}!m_{2}!\cdots
m_{n}!}\text{p}({\bf n})$

where $m_{\ell}$ is the number of $n_{j}$’s equal to $\ell$, but we restrict
ourselves mainly to the function p for simplicity.

The ratios after calculation for $n_{i}>1$ and $n_{i}=1$ are

$\displaystyle\frac{n_{i}}{n}\cdot\frac{1}{n-1}\frac{\omega_{i}^{(n-1)}({\bf
n}-{\bf e}_{i})}{\text{p}({\bf n}-{\bf e}_{i})}$
$\displaystyle=\frac{n_{i}}{n(n-1)}\cdot\frac{\int_{0}^{\infty}\frac{v^{n-1}}{(\psi(v))^{r+k}}\prod_{j=1}^{k}\pi_{n_{j}}(v)dv}{\int_{0}^{\infty}\frac{v^{n-2}}{(\psi(v))^{r+k}}\prod_{j=1}^{k}\pi_{n_{j}-\delta_{ij}}(v)dv}$
(5.6)

and

$\displaystyle\frac{1}{n}\frac{\omega_{0}^{(n-1)}({\bf n}-{\bf
e}_{i})}{\text{p}({\bf n}-{\bf e}_{i})}$
$\displaystyle=\frac{r+k-1}{n(n-1)}\cdot\frac{\int_{0}^{\infty}\frac{v^{n-1}}{(\psi(v))^{r+k}}\prod_{j=1}^{k}\pi_{n_{j}}(v)dv}{\int_{0}^{\infty}\frac{v^{n-2}}{(\psi(v))^{r+k-1}}\prod_{j=1,j\neq
i}^{k}\pi_{n_{j}}(v)dv}.$ (5.7)

Note that the numerator in (5.6) is the same in both the ratios, so can be
factored out of the sums in (5) and connected with $\varphi({\bf n})$. The two
ratios in (5.6) and (5.7) can be written as

$\displaystyle\frac{n_{i}}{n}\cdot\frac{\int_{0}^{\infty}g(v;{\bf
n})dv}{\int_{0}^{\infty}g(v;{\bf n}-{\bf e}_{j})dv}\quad{\rm
and}\quad\frac{1}{n}\cdot\frac{\int_{0}^{\infty}g(v;{\bf
n})dv}{\int_{0}^{\infty}g(v;{\bf n}-{\bf e}_{j})dv}.$

The two ratios are the same, however the meaning is different because a
singleton is removed in the second ratio. For comparison, in the two-parameter
Pitman-Yor model the two terms in (5.6) and (5.7) are

$\frac{n_{i}(n_{i}-\alpha-1)}{n(\theta+n-1)}\text{~{}and~{}}\frac{\theta+(k-1)\alpha}{n(\theta+n-1)}$

respectively. These have an important meaning in duality with the two-
parameter diffusion.

Acknowledgement We are very grateful to Yuguang Ipsen for her help in the
initial stages of this investigation.

##  6 Appendix: Proofs

It is well known (e.g. (James2005)) that the conditional distribution of a
Poisson point measure $\mathbb{B}$ with intensity $\rho$ given that a point
$S$ is observed equals $\mathbb{B}_{S}=\widetilde{\mathbb{B}}+\delta_{S}$,
where $\widetilde{\mathbb{B}}$ is again a Poisson point process with the same
distribution as, but independent of, $\mathbb{B}$. We denote the law of
$\mathbb{B}$ with intensity $\rho$ by $\mathcal{P}({\rm d}\mathbb{B}|\rho)$
and the Palm measure of $\mathbb{B}$, conditional on observing $S$, as
$\mathcal{P}_{S}({\rm d}\mathbb{B}|\rho)$, which is the law of
$\mathbb{B}_{S}$.

The Palm measure for the negative binomial point process,
$\mathcal{BN}(r,\rho)$, also has a neat form. Given that a point $S$ is
observed in $\mathbb{B}^{(r)}$, the conditional distribution of
$\mathbb{B}^{(r)}$, denoted by $\mathbb{B}^{(r)}_{S}$, can be written as
$\mathbb{B}^{(r)}_{S}=\widetilde{\mathbb{B}}^{(r+1)}+\delta_{S}$, where
$\widetilde{\mathbb{B}}^{(r+1)}$ is distributed as $\mathcal{BN}(r+1,\rho)$,
independent of $\mathbb{B}^{(r)}$ and $S$. See Prop. 4.3 in Gregoire (1984)
and Prop. 3.1 in Ipsen, Maller & Shemehsavar (2020a).

###### Proof of Theorem 2.1.

We make use of the characterisation that a negative binomial process
$\mathcal{BN}(r,\rho,G_{0})$ is equal in distribution to a Poisson point
process with intensity measure $\Gamma_{r}\rho\times G_{0}$, where
$\Gamma_{r}$ is an independent Gamma(r,1) rv. Conditional on $\Gamma_{r}=w$,
let $N=\sum_{i\geq 1}\delta_{J_{i},Y_{i}}$ be a Poisson point process with
parameter $w\rho(s)G_{0}(x)$. We will use the notation $\mathcal{P}({\rm
d}N|\Gamma_{r}\rho\cdot G_{0})$ for the law of $N$ and $\mathcal{P}({\rm
d}N|w\rho\cdot G_{0})$ for the law of $N$ conditional on $\Gamma=w$.

The observed data $({\bf n},{\bf Y})$ are samples from $\mathbb{M}^{(r)}$,
which can be written with latent variables $(N,{\bf J},{\bf Y},\Gamma_{r})$.
We start with Eq. (23) of (James, Lijoi & Pr$\ddot{\rm u}$nster2009) (which we
write with our variable $v$ rather than their $u$ and our initial distribution
$G_{0}$ rather than their $H$) to get the joint distribution of $(N,{\bf
J},{\bf Y},V_{n})$ as proportional to:

$\prod_{i=1}^{k}s_{i}^{n_{i}}e^{-vs_{i}}\rho({\rm d}s_{i}|Y_{i})G_{0}({\rm
d}Y_{i})v^{n-1}e^{-\psi(v)}\mathcal{P}_{\nu_{v}}({\rm d}N).$ (6.1)

In our case we take this conditional on $\Gamma_{r}=w$, let $\rho({\rm
d}s|Y)=w\rho(s){\rm d}s$ (since we assume the homogeneous case) and note that,
with $f(s,v)=vs$,

$\mathbb{E}(e^{-N(f)}|w\rho)=\exp\Big{(}\int_{s>0}(1-e^{-f(s,v)})w\rho(s){\rm
d}s\Big{)}=e^{-\int_{s>0}(1-e^{-vs})w\rho(s){\rm d}s}=e^{-w\psi(v)}$

(so we substitute $e^{-w\psi(v)}$ for the $e^{-\psi(v)}$ in (6.1)). Eq. (2) of
(James2005) gives a formula for a disintegration of a Poisson process $N$ with
point $J$ in terms of the Palm distribution: $N({\rm d}J)\mathcal{P}({\rm
d}N|\nu)=\mathcal{P}({\rm d}N|\nu,J)\nu({\rm d}J)$. In our case
$\mathcal{P}_{\nu_{v}}({\rm d}N)$ in (6.1) takes the form $\mathcal{P}_{{\bf
J},{\bf Y}}({\rm d}N|we^{-vs}\rho(s)\cdot G_{0})$. Using this we can write the
joint likelihood of $\\{{\bf N}_{[n]},\,{\bf J},{\bf Y},\Gamma_{r}\\}$ as
proportional to

$\prod_{i=1}^{k}\big{[}s_{i}^{n_{i}}e^{-vs_{i}}w\rho(s_{i}){\rm
d}s_{i}G_{0}({\rm
d}Y_{i})\big{]}\times\frac{v^{n-1}}{\Gamma(n)}e^{-w\psi(v)}\times\mathcal{P}_{{\bf
J},{\bf Y}}({\rm d}N|we^{-vs}\rho(s)\cdot
G_{0})\times\frac{w^{r-1}}{\Gamma(r)}e^{-w}.$ (6.2)

After some rearrangement, integrate with respect to $\\{{\bf
J},N,\Gamma_{r}\\}$, to get the joint distribution of $\\{{\bf N}_{[n]},{\bf
Y}\\}$ as

$\displaystyle\Big{[}\prod_{i=1}^{k}\int_{s_{i}}s_{i}^{n_{i}}e^{-vs_{i}}\rho(s_{i}){\rm
d}s_{i}\Big{]}\times\prod_{i=1}^{k}G_{0}(dY_{i})\times\frac{v^{n-1}}{\Gamma(n)}\frac{\Gamma(r+k)}{\Gamma(r)}$
(6.3) $\displaystyle\hskip 56.9055pt\times\int_{{\cal
M}}\int_{w}\mathcal{P}_{{\bf J},{\bf Y}}({\rm d}N|we^{-vs}\rho(s)\cdot
G_{0})\times\frac{w^{r+k-1}}{\Gamma(r+k)}e^{-w\psi(v)}e^{-w}{\rm d}w.$ (6.4)

The integral in the second line of (6.3) equals

$(\psi(v))^{-(r+k)}\int_{{\cal M}}\mathcal{P}_{{\bf J},{\bf Y}}({\rm
d}N|\Gamma_{r+k,\psi(v)}e^{-vs}\rho(s)\cdot G_{0})$

where $\Gamma_{r+k,\psi(v)}$ is an independent Gamma$(r+k,\psi(v))$ random
variable.

Note that a Poisson point process with intensity measure
$\Gamma_{r,\beta}\rho$ is in $\mathcal{BN}(r,\rho/\beta)$, as we can see by
checking that the Laplace functionals are the same; namely,

$\displaystyle
E\exp\Big{(}-\int_{0}^{\infty}(1-e^{-f(x)})\Gamma_{r,\beta}\rho(x){\rm
d}x\Big{)}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\exp\Big{(}-\int_{0}^{\infty}(1-e^{-f(x)})y\rho(x){\rm
d}x\Big{)}\beta e^{-\beta y}\frac{(\beta y)^{r-1}}{\Gamma(r)}{\rm d}y$
$\displaystyle=$
$\displaystyle\beta^{r}\int_{0}^{\infty}\exp\big{(}-y\int_{0}^{\infty}(1-e^{-f(x)})\rho(x){\rm
d}x+\beta\big{)}\frac{y^{r-1}}{\Gamma(r)}{\rm d}y$ $\displaystyle=$
$\displaystyle\beta^{r}\Big{(}\beta+\int_{0}^{\infty}(1-e^{-f(x)})\rho(x){\rm
d}x\Big{)}^{-r}=\Big{(}1+\frac{1}{\beta}\int_{0}^{\infty}(1-e^{-f(x)})\rho(x){\rm
d}x\Big{)}^{-r},$

with $f$ any bounded nonnegative measurable function on $\mathbb{R}_{+}$. This
is of the form in (1.4). Thus, we can write

$\mathcal{P}({\rm d}N|\Gamma_{r+k,\psi(v)}e^{-vs}\rho(s)\cdot
G_{0})=\mathcal{P}\Big{(}{\rm
d}N|\Gamma_{r+k}\frac{e^{-vs}\rho(s)}{\psi(v)}\cdot G_{0}\Big{)}.$ (6.5)

Since the integral over the space of point measures ${\cal M}$ is $1$, we get
from (6.3) the joint distribution of $V_{n}$ and ${\bf X}$ as

$\Big{[}\prod_{i=1}^{k}\pi_{n_{i}}(v)\Big{]}\times\prod_{i=1}^{k}G_{0}({\rm
d}Y_{i})\times\frac{v^{n-1}}{\Gamma(n)}r^{[k]}\times(\psi(v))^{-(r+k)}=\prod_{i=1}^{k}G_{0}({\rm
d}Y_{i})\times g_{r}(v,{\bf n}),$ (6.6)

where $g_{r}(v,{\bf n})$ is defined in (2.4). Dividing this into (6.2) gives
the joint density of $\\{{\bf J},N,\Gamma_{r}\\}$ conditional on $\\{{\bf
X},V_{n}\\}$ as

$\frac{(\psi(v))^{r+k}}{[r]^{k}}\prod_{i=1}^{k}\frac{s_{i}^{n_{i}}e^{-vs_{i}}\rho(s_{i})}{\pi_{n_{i}}(v)}\frac{w^{r+k-1}}{\Gamma(r)}e^{-w(\psi(v))}\times\mathcal{P}_{{\bf
J},{\bf Y}}({\rm d}N|we^{-vs}\rho(s)\cdot G_{0}).$ (6.7)

Integrate out $w$ from (6.7) to obtain the joint density of $\\{{\bf J},N\\}$
given ${\bf X}$ and $V_{n}$ as

$\prod_{i=1}^{k}\frac{s_{i}^{n_{i}}e^{-vs_{i}}\rho(s_{i})}{\pi_{n_{i}}(v)}\times\mathcal{P}_{{\bf
J},{\bf Y}}\Big{(}{\rm d}N|\Gamma_{r+k}\frac{e^{-vs}\rho(s)}{\psi(v)}\cdot
G_{0}\Big{)}.$

Then note that the Palm measure of a Poisson measure with law
$\mathcal{P}_{{\bf J},{\bf Y}}({\rm d}N|w\frac{e^{-vs}\rho(s)}{\psi(v)}\cdot
G_{0})$ has the same distribution as

$N^{v,w}+\sum_{i=1}^{k}\delta_{J_{i}^{v},Y_{i}},\quad\text{where $N^{v,w}$ has
law $\mathcal{P}\Big{(}{\rm d}N|w\frac{e^{-vs}\rho(s)}{\psi(v)}\cdot
G_{0}\Big{)}$},$

and the $J_{i}^{v}$ have the density in (2.7). From here we can read off the
posterior distribution of $\\{N,{\bf J}\\}$ given ${\bf X}$ and $V_{n}$. We
obtain (2.5) by observing that $\mathcal{P}\Big{(}{\rm
d}N|\Gamma_{r+k}\frac{e^{-vs}\rho(s)}{\psi(v)}\cdot G_{0}\Big{)}$ is the law
of $\mathcal{BN}(r+k,\frac{e^{-vs}\rho(s)}{\psi(v)},G_{0})$. ∎

###### Proof of Theorem 3.1..

Let $G$ be the species sampling model derived from the negative binomial
process as in (2.1). Let ${\bf X}$ be an independent sample of size $n$ from
$G$, then, following the proof of Prop. 2, p.96, of (James, Lijoi &
Pr$\ddot{\rm u}$nster2009),

$\displaystyle\mathbb{P}(X_{n+1}\in{\rm d}x|{\bf X})$
$\displaystyle=\int\mathbb{P}(X_{n+1}\in{\rm d}x|G)\mathbb{P}(G|{\bf X})$
$\displaystyle=\int G({\rm d}x)\mathbb{P}(G|{\bf X})$
$\displaystyle=\mathbb{E}(G_{{\bf X}}({\rm d}x))=\int_{v}\mathbb{E}(G_{{\bf
X}}^{v}({\rm d}x))\,g_{r}(v,{\bf n}){\rm d}v.$

Here $G_{\bf X}$ is the posterior distribution of $G$ given ${\bf X}$, $G_{\bf
X}^{V_{n}}$ is specified in Theorem 2.2 and $g_{r}(v,{\bf n})$ is defined in
(2.4). Conditioning on $V_{n}=v$, by (2.9) we have

$\mathbb{E}(G_{{\bf X}}^{v}({\rm d}x))=\mathbb{E}\Bigg{(}\frac{\mu^{v}({\rm
d}x)}{T^{v}+\sum_{i=1}^{k}J_{i}^{v}}\Bigg{)}+\mathbb{E}\Bigg{(}\frac{\sum_{i=1}^{k}J_{i}^{v}\delta_{Y_{i}}({\rm
d}x)}{T^{v}+\sum_{i=1}^{k}J_{i}^{v}}\Bigg{)}:=I_{1}(v,{\rm d}x)+I_{2}(v,{\rm
d}x),$ (6.8)

with $\mu^{v}$ and $T^{v}$ as in Theorem 2.2 and $(J_{i}^{v},Y_{i})$ as in
Proposition 2.1. We evaluate each component in (6.8) separately.

Write as a shorthand $R:=\Gamma_{r+k,\psi(v)}$. Then by Theorem 2.2 and using
(6.5) for the distribution of $\mathbb{M}^{(r+k,v)}$, we get

$\displaystyle I_{1}(v,{\rm d}x)=$
$\displaystyle\mathbb{E}\Bigg{(}\frac{\mu^{v}({\rm
d}x)}{T^{v}+\sum_{i=1}^{k}J_{i}^{v}}\Bigg{)}$ $\displaystyle=$
$\displaystyle\mathbb{E}\int_{\lambda}e^{-\lambda
T^{v}}\int_{s}s\mathbb{M}^{(r+k,v)}({\rm d}s,{\rm
d}x)e^{-\lambda\sum_{i}J_{i}^{v}}{\rm d}\lambda$ $\displaystyle=$
$\displaystyle\int_{w}\int_{\cal M}\int_{\lambda}\int_{s}sN({\rm d}s,{\rm
d}x)\cdot e^{-\lambda T^{v}}\mathcal{P}({\rm d}N|we^{-vs}\rho(s)\cdot
G_{0})\cdot\mathbb{E}(e^{-\lambda\sum_{i}J_{i}^{v}}){\rm
d}\lambda\mathbb{P}(R\in{\rm d}w)$ (6.9)

Here note that by Eq. (2) of (James2005),

$e^{-\lambda T^{v}}\mathcal{P}({\rm d}N|we^{-vs}\rho(s)\cdot
G_{0})=\mathbb{E}(e^{-\lambda T^{v}}|we^{-vs}\rho(s)\cdot
G_{0})\mathcal{P}({\rm d}N|we^{-(\lambda+v)s}\rho(s)\cdot G_{0}),$

where $\mathbb{E}(e^{-\lambda T^{v}}|we^{-vs}\rho(s)\cdot
G_{0})=\exp(-w\psi^{(v)}(\lambda))$, and

$\displaystyle\psi^{(v)}(\lambda)$ $\displaystyle=$
$\displaystyle\int(1-e^{-\lambda x})e^{-vx}\rho(x){\rm
d}x=\int(1-e^{-(\lambda+v)x})\rho(x){\rm d}x-\int(1-e^{-vx})\rho(x){\rm d}x$
$\displaystyle=$ $\displaystyle\psi(\lambda+v)-\psi(v).$

Also

$\mathbb{E}(e^{-\lambda\sum_{i}J_{i}^{v}})=\prod_{i=1}^{k}\frac{1}{\pi_{n_{i}}(v)}\int_{s}e^{-\lambda
s}e^{-vs}s^{n_{i}}\rho(s){\rm
d}s=\prod_{i=1}^{k}\frac{\pi_{n_{i}}(\lambda+v)}{\pi_{n_{i}}(v)}.$

Thus, (6) equals

$\displaystyle\int_{w}\int_{\lambda}\int_{s}s\int_{\cal M}N({\rm d}s,{\rm
d}x)\mathcal{P}({\rm d}N|we^{-(\lambda+v)s}\rho(s)G_{0})\cdot
e^{-w(\psi(\lambda+v)-\psi(v))}$ (6.10) $\displaystyle\hskip
199.16928pt\times\prod_{i=1}^{k}\frac{\pi_{n_{i}}(\lambda+v)}{\pi_{n_{i}}(v)}{\rm
d}\lambda\mathbb{P}(R\in{\rm d}w).$ (6.11)

Also note that

$\displaystyle\int_{s}s\int_{\cal M}N({\rm d}s,{\rm d}x)\mathcal{P}({\rm
d}N|we^{-(\lambda+v)s}\rho(s)\cdot G_{0})$ $\displaystyle=$ $\displaystyle
G_{0}({\rm d}x)\cdot w\int_{s}e^{-(\lambda+v)s}s\rho({\rm d}s)$
$\displaystyle=$ $\displaystyle w\pi_{1}(v+\lambda)G_{0}({\rm d}x).$

Recalling the density $g_{r}(v,{\bf n})$ in (2.4) and that
$R=\Gamma_{r+k,\psi(v)}$, we see that

$\displaystyle\int_{v}I_{1}(v,{\rm d}x)g_{r}(v,{\bf n}){\rm
d}v=\int_{v}\mathbb{E}\Bigg{(}\frac{\mu^{v}({\rm
d}x)}{T^{v}+\sum_{i=1}^{k}J_{i}^{v}}\Bigg{)}g_{r}(v,{\bf n}){\rm d}v$
$\displaystyle=G_{0}({\rm
d}x)\int_{v}\int_{\lambda}\prod_{i=1}^{k}\frac{\pi_{n_{i}}(v+\lambda)}{\pi_{n_{i}}(v)}\int_{w}w\pi_{1}(v+\lambda)e^{-w(\psi(\lambda+v)-\psi(v))}\frac{w^{r+k-1}}{\Gamma(r+k)}e^{-w\psi(v)}{\rm
d}w\,{\rm d}\lambda$ $\displaystyle\hskip
85.35826pt\times\prod_{i=1}^{k}\pi_{n_{i}}(v)\frac{v^{n-1}}{\Gamma(n)}r^{[k]}(\psi(v))^{r+k}(\psi(v))^{-(r+k)}{\rm
d}v.$ (6.12)

Collecting the $w$ terms together, this equals

$\displaystyle G_{0}({\rm
d}x)\int_{v}\int_{\lambda}\pi_{1}(v+\lambda)\prod_{i=1}^{k}\pi_{n_{i}}(v+\lambda)\int_{w}e^{-w(\psi(\lambda+v))}\frac{w^{r+k}}{\Gamma(r+k+1)}{\rm
d}w\,{\rm d}\lambda$ (6.13) $\displaystyle\hskip 199.16928pt\times
v^{n-1}\frac{\Gamma(r+k+1)}{\Gamma(r)\Gamma(n)}{\rm d}v$
$\displaystyle=G_{0}({\rm
d}x)\int_{v}\int_{\lambda}\pi_{1}(v+\lambda)\prod_{i=1}^{k}\pi_{n_{i}}(v+\lambda)(\psi(\lambda+v))^{-(r+k)-1}\,{\rm
d}\lambda\,v^{n-1}\frac{\Gamma(r+k+1)}{\Gamma(r)\Gamma(n)}{\rm d}v.$ (6.14)

Changing variable to $z=\lambda+v$ and recalling (2.4), we get the RHS of
(6.13) equal to

$\displaystyle G_{0}({\rm
d}x)\int_{z=0}^{\infty}\Big{(}\int_{v=0}^{z}v^{n-1}{\rm
d}v\Big{)}\pi_{1}(z)\prod_{i=1}^{k}\pi_{n_{i}}(z)(\psi(z))^{-(r+k)-1}\,\,\frac{\Gamma(r+k+1)}{\Gamma(r)\Gamma(n)}{\rm
d}z$ $\displaystyle=\frac{G_{0}({\rm
d}x)}{n}\int_{z=0}^{\infty}\frac{z\pi_{1}(z)(r+k)}{\psi(z)}\Big{\\{}\frac{z^{n-1}}{\Gamma(n)}\prod_{i=1}^{k}\pi_{n_{i}}(z)(\psi(z))^{-(r+k)}r^{[k]}\,\Big{\\}}{\rm
d}z.$

Comparing the last expression to (3.2) we see that

$\int_{v}I_{1}(v,{\rm d}x)g_{r}(v,{\bf n}){\rm d}v=\omega_{0}^{(n)}G_{0}({\rm
d}x).$

To get (3.3), we work with $I_{2}(v,{\rm d}x)$ in (6.8) in a similar way:

$\displaystyle I_{2}(v,{\rm d}x)=$
$\displaystyle\mathbb{E}\Bigg{(}\frac{\sum_{i=1}^{k}J_{i}^{v}\delta_{Y_{i}}({\rm
d}x)}{T^{v}+\sum_{i=1}^{K}J_{i}^{v}}\Bigg{)}$ $\displaystyle=$
$\displaystyle\int_{\lambda}\mathbb{E}(e^{-\lambda
T^{v}|\Gamma_{r+k,\psi(v)}e^{-v}\rho})\mathbb{E}\Big{(}e^{-\lambda\sum_{i=1}^{k}J_{i}^{v}}\sum_{i=1}^{k}J_{i}^{v}\delta_{Y_{i}}({\rm
d}x)\Big{)}{\rm d}\lambda.$ (6.16)

Integrating out the $\Gamma_{r+k,\psi(v)}$ in a similar way as for (6), we get
the first expectation in (6) as

$\displaystyle\mathbb{E}(e^{-\lambda T^{v}|\Gamma_{r+k,\psi(v)}e^{-v}\rho})$
$\displaystyle=$
$\displaystyle\int_{w}\exp\Big{(}-\int_{0}^{\infty}(1-e^{-\lambda
y})we^{-vy}\rho(y){\rm
d}y\Big{)}e^{-w\psi(v)}\frac{w^{r+k-1}(\psi(v))^{r+k}}{\Gamma(r+k)}{\rm d}w$
$\displaystyle=$
$\displaystyle\Big{(}\frac{\psi(v)}{\psi(\lambda+v)}\Big{)}^{r+k}.$

Recalling the density of $(J_{i}^{v})$ in (2.7), the second expectation in (6)
can be simplified as

$\displaystyle\mathbb{E}\Big{(}e^{-\lambda\sum_{i=1}^{k}J_{i}^{v}}\sum_{i=1}^{k}J_{i}^{v}\delta_{Y_{i}}({\rm
d}x)\Big{)}$ $\displaystyle=$
$\displaystyle\mathbb{E}\Big{(}\sum_{i=1}^{k}e^{-\lambda\sum_{j}J_{j}^{v}}J_{i}^{v}\delta
Y_{i}({\rm d}x)\Big{)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{k}\int_{s_{1}}\cdots\int_{s_{k}}\prod_{j=1}^{k}e^{-\lambda
s_{j}}s_{i}\delta_{Y_{i}}({\rm
d}x)\prod_{j=1}^{k}\frac{e^{-vs_{j}}s_{j}^{n_{j}}\rho(s_{j})}{\pi_{n_{j}}(v)}{\rm
d}s_{1}\cdots{\rm d}s_{k}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{k}\int_{s_{1}}\cdots\int_{s_{k}}\Big{(}\prod_{j=1}^{k}\frac{e^{-(v+\lambda)s_{j}}s_{j}^{n_{j}}\rho(s_{j})}{\pi_{n_{j}}(v)}\Big{)}s_{i}\delta_{Y_{i}}({\rm
d}x)\,{\rm d}s_{1}\cdots{\rm d}s_{k}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{k}\prod_{j\neq
i}\frac{\pi_{n_{j}}(v+\lambda)}{\pi_{n_{j}}(v)}\int_{s_{i}}s_{i}\delta_{Y_{i}}({\rm
d}x)\,\frac{e^{-(v+\lambda)s_{i}}s_{i}^{n_{i}}\rho(s_{i})}{\pi_{n_{i}}}{\rm
d}s_{i}$ $\displaystyle=$ $\displaystyle\sum_{i=1}^{k}\prod_{j\neq
i}\frac{\pi_{n_{j}}(v+\lambda)}{\pi_{n_{j}}(v)}\frac{\pi_{n_{i}+1}(\lambda+v)}{\pi_{n_{i}(v)}}\,\delta_{Y_{i}}({\rm
d}x)\,$ $\displaystyle=$
$\displaystyle\prod_{j=1}^{k}\frac{1}{\pi_{n_{j}}(v)}\sum_{i=1}^{k}\Big{(}\pi_{n_{i}+1}(\lambda+v)\prod_{j\neq
i}\pi_{n_{j}}(v+\lambda)\Big{)}\,\delta_{Y_{i}}({\rm d}x).$

Putting these together and recalling $g_{r}(v,{\bf n})$ in (2.4), we get

$\displaystyle\int_{v}I_{2}(v,{\rm d}x)g_{r}(v,{\bf n}){\rm d}v$
$\displaystyle=\int_{v}\int_{\lambda}\Big{(}\frac{\psi(v)}{\psi(\lambda+v)}\Big{)}^{r+k}\prod_{j=1}^{k}\frac{1}{\pi_{n_{j}}(v)}\sum_{i=1}^{k}\Big{(}\pi_{n_{i}+1}(\lambda+v)\prod_{j\neq
i}\pi_{n_{j}}(v+\lambda)\Big{)}\,\delta_{Y_{i}}({\rm d}x){\rm d}\lambda$
$\displaystyle\hskip
144.54pt\times\prod_{i=1}^{k}\pi_{n_{i}}(v)\frac{v^{n-1}}{\Gamma(n)}r^{[k]}(\psi(v))^{-(r+k)}{\rm
d}v$
$\displaystyle=\int_{v}\int_{\lambda}\sum_{i=1}^{k}\Big{(}\pi_{n_{i}+1}(\lambda+v)\prod_{j\neq
i}\pi_{n_{j}}(v+\lambda)\Big{)}\,\delta_{Y_{i}}({\rm d}x){\rm
d}\lambda\times\frac{v^{n-1}}{\Gamma(n)}r^{[k]}(\psi(\lambda+v))^{-(r+k)}{\rm
d}v$ $\displaystyle=\int_{z=0}^{\infty}\Big{(}\int_{v=0}^{z}v^{n-1}{\rm
d}v\Big{)}\sum_{i=1}^{k}\Big{(}\pi_{n_{i}+1}(z)\prod_{j\neq
i}\pi_{n_{j}}(z)\Big{)}\,\delta_{Y_{i}}({\rm
d}x)\times\frac{1}{\Gamma(n)}r^{[k]}(\psi(z))^{-(r+k)}{\rm d}z$
$\displaystyle=\int_{z=0}^{\infty}\frac{z}{n}\Big{(}\prod_{j=1}^{k}\pi_{n_{j}}(z)\frac{z^{n-1}}{\Gamma(n)}r^{[k]}(\psi(z))^{-(r+k)}\Big{)}\sum_{i=1}^{k}\Big{(}\frac{\pi_{n_{i}+1}(z)}{\pi_{n_{i}}(z)}\Big{)}\,\delta_{Y_{i}}({\rm
d}x){\rm d}z$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{k}\int_{z}z\frac{\pi_{n_{i}+1}(z)}{\pi_{n_{i}}(z)}g_{r}(z,{\bf
n}){\rm d}z\,\delta_{Y_{i}}({\rm d}x).$ (6.17)

We recognise (6) as the second term on the RHS of (3.1). This completes the
proof. ∎

###### Proof of Proposition 3.1.

Using (2.4) the RHS of (3.7) is

$\displaystyle\int_{0}^{\infty}v\Big{(}\frac{(r+k)}{n}\frac{\pi_{1}(v)}{\psi(v)}+\frac{1}{n}\sum_{i=1}^{k}\frac{\pi_{n_{i}+1}(v)}{\pi_{n_{i}}(v)}\Big{)}g_{r}(v,{\bf
n}){\rm d}v$
$\displaystyle=\frac{r^{[k]}}{n\Gamma(n)}\int_{0}^{\infty}\Big{(}\frac{(r+k)\pi_{1}(v)}{\psi(v)}+\sum_{i=1}^{k}\frac{\pi_{n_{i}+1}(v)}{\pi_{n_{i}}(v)}\Big{)}\frac{v^{n}\prod_{i=1}^{k}\pi_{n_{i}}(v)}{(\psi(v))^{r+k}}.$

Leaving aside the factor of $r^{[k]}/n\Gamma(n)$, the integral term equals

$\displaystyle\int_{0}^{\infty}\Big{(}\frac{(r+k){\rm d}\psi(v)/{\rm
d}v}{\psi(v)}+\sum_{i=1}^{k}\frac{1}{\pi_{n_{i}}(v)}\big{(}-\frac{{\rm
d}\pi_{n_{i}}}{{\rm
d}v}\Big{)}\frac{v^{n}\prod_{i=1}^{k}\pi_{n_{i}}(v)}{(\psi(v))^{r+k}}{\rm d}v$
$\displaystyle=\int_{0}^{\infty}\Big{(}(r+k)\frac{{\rm d}\log\psi(v)}{{\rm
d}v}-\sum_{i=1}^{k}\frac{{\rm d}\log\pi_{n_{i}}}{{\rm
d}v}\Big{)}\frac{v^{n}\prod_{i=1}^{k}\pi_{n_{i}}(v)}{(\psi(v))^{r+k}}{\rm d}v$
$\displaystyle=\int_{0}^{\infty}\frac{{\rm d}}{{\rm
d}v}\log\Big{(}\frac{(\psi(v))^{r+k}}{\prod_{i=1}^{k}\pi_{n_{i}}(v)}\Big{)}\frac{v^{n}\prod_{i=1}^{k}\pi_{n_{i}}(v)}{(\psi(v))^{r+k}}{\rm
d}v.$

Reintroducing the factor of $r^{[k]}/n\Gamma(n)$ and defining

$f(v):=\log\Big{(}\frac{(\psi(v))^{r+k}}{\prod_{i=1}^{k}\pi_{n_{i}}(v)}\Big{)},$

the last expression can be written as

$\frac{r^{[k]}}{n\Gamma(n)}\int_{0}^{\infty}\frac{{\rm d}f(v)}{{\rm
d}v}v^{n}e^{-f(v)}{\rm
d}v=\frac{r^{[k]}}{\Gamma(n)}\int_{0}^{\infty}v^{n-1}e^{-f(v)}{\rm d}v,$

with an integration by parts at the last step. Substituting for $f(v)$, the
RHS is seen to equal $\int_{0}^{\infty}g_{r}(v,{\bf n}){\rm d}v$.

(ii) Using (2.4) and (3.6) we can write

$\frac{n!}{\prod_{j=1}^{n}j!^{m_{j}}m_{j}!}g_{r}(v,{\bf
m})=\frac{n!r^{[k]}v^{n-1}}{\Gamma(n)(\psi(v))^{r+k}}\prod_{j=1}^{n}\frac{(\pi_{j}(v))^{m_{j}}}{j!^{m_{j}}m_{j}!}.$

Let $f(x)=x^{-r}$, so that $f^{(k)}(x)=r^{[k]}(-1)^{k}x^{-r-k}$, and let
$\phi(v)=\psi(v)$. Then $\pi_{j}(v)=(-1)^{j-1}\phi^{(j)}(v)$ and

$\displaystyle\sum_{k=1}^{n}\sum_{{\bf m}\in
A_{nk}}\frac{n!}{\prod_{j=1}^{n}j!^{m_{j}}m_{j}!}g_{r}(v,{\bf m})$
$\displaystyle=\frac{(-1)^{n}v^{n-1}}{\Gamma(n)}\sum_{k=1}^{n}\sum_{{\bf m}\in
A_{nk}}n!f^{(k)}(\phi(v))\prod_{j=1}^{n}\frac{1}{m_{j}!}\Big{(}\frac{(-1)^{j}\phi^{(j)}(v)}{j!}\Big{)}^{m_{j}}.$

Recall from $A_{nk}$ in (3.5) that $\sum_{j=1}^{n}jm_{j}=n$ and
$\sum_{j=1}^{n}m_{j}=k$. Then, after cancelling the factors of $-1$, we
recognise in the sum and product Faà di Bruno’s formula for
$f^{(n)}(\phi(v))$, so the integral over $v$ of the expression equals

$\frac{1}{\Gamma(n)}\int_{0}^{\infty}v^{n-1}f^{(n)}(\phi(v)){\rm d}v,$

and this is seen to equal 1 after $n-1$ integrations by parts. Hence (3.8). ∎

## References

  * Chegini & Zarepour (2023) Chegini, S. & Zarepour, M. (2023)  Random discrete probability measures based on negative binomial process.  arxiv.org/pdf/2307.00176.
  * De Iorio & Griffiths (2004) De Iorio, M. & Griffiths, R.C. (2004)  Importance sampling on coalescent histories. I  Adv. Appl. Probab., 36, 434–454.
  * Ethier & Griffiths (1987) Ethier, S.N. & Griffiths. R.C. (1987)  The infinitely-many-sites model as a measure-valued diffusion.  Ann. Probab., 15, 515–545.
  * Ewens (1972) Ewens, W. (1972)  The sampling theory of selectively neutral alleles.  Theoret. Pop. Biol., 3, 87–112.
  * Felsensteinetal (1999) Felsenstein, J., Kuhner, M.K., Yamato, J. & Beerli, P (1999)  Likelihoods on coalescents: A Monte Carlo sampling approach to inferring parameters from population samples of molecular data.  In: Statistics in Molecular Biology and Genetics, Seillier-Moiseiwitsch, F. (Ed.), Inst. Math. Statistics & American Math. Soc., Haywood, CA 163–185.
  * Gregoire (1984) Gregoire, G. (1984)  Negative binomial distributions for point processes.  Stoch. Proc. Appl., 16, 179–188.
  * Griffiths (1979) Griffiths, R.C. (1979)  On the distribution of allele frequencies in a diffusion model.  Theoret. Pop. Biol., 15, 140–158.
  * Griffiths & Tavaré (1994) Griffiths, R.C. & Tavaré, S. (1994)  Simulating probability distributions in the coalescent.  Theoret. Popn. Biol., 46, 131–159.
  * Griffiths, Ruggiero, Spano, and Zhou (2023) Griffiths, R.C., Ruggiero, M., Spano, D. and Zhou, Y., (2023) Dual process in the two-parameter Poisson-Dirichlet diffusion. arXiv preprint arXiv:2102.08520.
  * Hobolth, Uyenoyama & Wiuf (2008) Hobolth, A., Uyenoyama, M.K., Wiuf, C. (2008)  Importance sampling for the infinite sites model.  Statist. Appl. Genetics Molecular Biology, Vol. 7: Issue 1, Article 32.
  * Hoppe (1984) Hoppe, F.M. (1984)  Pólya-like urns and the Ewens’ sampling formula.  J. Math. Biol. 20, 91–94.
  * Hoppe (1987) Hoppe, F.M. (1987)  The sampling theory of neutral alleles and an urn model in population genetics.  J. Math. Biol. 25, 123–159.
  * Hudson (1991) Hudson, R.R. (1991)  Gene genealogies and the coalescent process.  Oxford Surveys in Evolutionary Biology 7, D. Futuyma and J. Antonovics (eds.), pp. 1-44.
  * Ipsen & Maller (2017) Ipsen, Y.F. & Maller, R.A. (2017)  Negative binomial construction of random discrete distributions on the infinite simplex.  Theor. Stoch. Proc., 22, 34–46.
  * Ipsen, Shemehsavar & Maller (2018) Ipsen, Y.F., Shemehsavar, S. & Maller, R.A. (2018)  Species sampling models generated by negative binomial processes.  arXiv: 1904.13046.
  * Ipsen, Maller & Shemehsavar (2020a) Ipsen, Y.F., Maller, R.A. & Shemehsavar (2020a)  Limiting distributions of generalised Poisson-Dirichlet distributions based on negative binomial processes.  J. Theoret. Prob., 33, 1974–2000.
  * Ipsen, Maller & Shemehsavar (2020b) Ipsen, Y.F., Maller, R.A. & Shemehsavar, S. (2020b)  Size biased sampling from the Dickman subordinator.  Stoch. Proc. Appl., 130, 6880–6900.
  * Ipsen, Maller & Shemehsavar (2021) Ipsen, Y.F., Maller, R.A. & Shemehsavar, S. (2021)  A generalised Dickman distribution and the number of species in a negative binomial process model.  Adv. Appl. Prob., 53, 370–399.
  * Karlin & McGregor (1972) Karlin, S. and McGregor, J. (1972)  Addendum to a paper of W. Ewens.  Theoret. Popn. Biol., 3, 113–116.
  * (20) James, Lancelot F. (2002)  Poisson process partition calculus with applications to exchangeable models and Bayesian nonparametrics.  arXiv:0205093.
  * (21) James, Lancelot F. (2005)  Bayesian Poisson process partition calculus with an application to Bayesian Lévy moving averages.  The Annals of Statistics, 33, 1771–1799.
  * (22) James, Lancelot F. (2008)  Large sample asymptotics for the two-parameter Poisson–Dirichlet process.  Pushing the Limits of Contemporary Statistics: Contributions in Honor of Jayanta K. Ghosh, Vol. 3 (2008) 187–199, Inst. Math. Statistics, 2008.
  * (23) James, L. F., Lijoi, A. & Pr$\ddot{\rm u}$nster, I. (2009). Posterior analysis for normalized random measures with independent increments.  Scand. J. Statistics, 36, 76–97.
  * Kingman (1975) Kingman, J.F.C. (1975)  Random discrete distributions.  J. Roy. Statist. Soc. B, 37, 1–22.
  * Labadi & Zarepour (2014) Labadi, L.A. & Zarepour, M. (2014)  On simulations from the two-parameter Poisson-Dirichlet process and the normalized Inverse-Gaussian process.  Sankhyã, 76-A, Part 1, 158–176.
  * Maller & Shemehsavar (2023a) Maller, R.A. & Shemehsavar, S. (2023a)  Generalized Poisson–Dirichlet distributions based on the Dickman subordinator.  Theor. Prob. Appl., 67, 593-612.
  * Maller & Shemehsavar (2023b) Maller, R.A. & Shemehsavar, S. (2023b)  Limiting distributions of Poisson-Dirichlet and generalized Poisson–Dirichlet distributions.  arXiv:2309.00815 [math.ST] https://doi.org/10.48550/arXiv.2309.00815
  * Petrov (2009) Petrov, L. (2009) Two-parameter family of diffusion processes in the Kingman simplex. _Funct. Anal. Appl._ 43 279–296.
  * Pitman (2003) Pitman, J. (2003)  Poisson-Kingman partitions.  Lect. Notes Monogr. Ser., 40, 1–34.
  * Pitman (2006) Pitman, J. (2006)  Combinatorial Stochastic Processes.  Springer-Verlag, Berlin.
  * Pitman & Yor (1997) Pitman, J. & Yor, M. (1997)  The two-parameter Poisson–Dirichlet distribution derived from a stable subordinator.  Ann. Probab., 25, 855–900.
  * Stephens & Donnelly (2000) Stephens, M. & Donnelly, P. (2000)  Inference in molecular population genetics.  J. R. Stat. Soc., B 62, 605–635.