#  Is infinity that far?
A Bayesian nonparametric perspective of finite mixture models

Raffaele Argiento
ESOMAS Department, University of Torino and Collegio Carlo Alberto, Torino,
Italy
and
Maria De Iorio
Yale-NUS College, Singapore and Dept. of Statistical Science, University
College London

###### Abstract

Mixture models are one of the most widely used statistical tools when dealing
with data from heterogeneous populations. This paper considers the long-
standing debate over finite mixture and infinite mixtures and brings the two
modelling strategies together, by showing that a finite mixture is simply a
realization of a point process. Following a Bayesian nonparametric
perspective, we introduce a new class of prior: the Normalized Independent
Point Processes. We investigate the probabilistic properties of this new
class. Moreover, we design a conditional algorithm for finite mixture models
with a random number of components overcoming the challenges associated with
the Reversible Jump scheme and the recently proposed marginal algorithms. We
illustrate our model on real data and discuss an important application in
population genetics.

Keywords: Bayesian Clustering, Bayesian Mixture Models, Conditional
Algorithms, Dirichlet process, Markov Chain Monte Carlo Methods.

##  1 Introduction

Mixture models are a very powerful and natural statistical tool to model data
from heterogeneous populations. In a mixture model, observations are assumed
to have arisen from one of $M$ (finite or infinite) groups, each group being
suitably modelled by a density typically from a parametric family. The density
of each group is referred to as a component of the mixture, and is weighted by
the relative frequency (weight) of the group in the population. This model
offers a conceptually simple way of relaxing distributional assumptions and a
convenient and flexible way to approximate distributions that cannot be
modelled satisfactorily by a standard parametric family. Moreover, it provides
a framework by which observations may be clustered together into groups for
discrimination or classification. For a comprehensive review of mixture models
and their applications see McLachlan et al. (2000); Frühwirth-Schnatter (2006)
and Fruhwirth-Schnatter et al. (2019). More in details, let
$Y\in\mathcal{Y}\subset{\mathbb{R}}^{r}$ be the population variable, each
observation is assumed to have arisen from one of $0<M\leq\infty$ groups:

$f_{Y}(y\mid
P)=\int_{\Theta}f(y\mid\theta)P(d\theta)=\sum_{j=1}^{M}w_{j}f(y\mid\tau_{j})$
(1)

where $\\{f(\cdot\mid\theta),\theta\in\Theta\subset{\mathbb{R}}^{d}\\}$ is a
parametric family of densities on $\mathcal{Y}$, while $P(\cdot)$ is an almost
sure discrete measure on $\Theta$, and it is referred to as _mixing measure_.
Here $\\{\tau_{j},j=1,\ldots,M\\}$ is a collection of points in $\Theta$, that
defines the support of $P$. For each $j=1,\dots,M$, the density
$f(y\mid\tau_{j})$ is the kernel of the mixture, and is weighted by $w_{j}$,
the relative frequency of the group in the population. Model (1) defines a
framework by which observations may be clustered together into groups, so that
conditionally, data are independent and identically distributed within the
groups and independent between groups. To avoid confusion in terminology, in
what follows $M$ will denote the number of components in a mixture, i.e. of
possible clusters/sub-populations, while by number of clusters, $k$, we mean
the number of allocated components, i.e. components to which at least one
observation has been assigned. The latter quantity can only be estimated a
posteriori.

We believe that in the context of mixture modelling the words cluster and
component are often misused in terminology, i.e. the distinction between
number of components and number of clusters has generally been overlooked in
the parametric world, leading to the criticism that if we fix a priori the
number $M$, we cannot estimate the number of clusters. What needs to be
highlighted (see Rousseau and Mengersen, 2011) is that when in a finite
mixture model we fix $M$, we are specifying the number of components (i.e.
possible clusters) that corresponds to the data generating process, but still
we need to estimate the actual number of clusters in the sample (allocated
components). Already Nobile et al. (2004) had pointed out this difference,
noticing that the posterior distribution of the number of components $M$ might
assigns considerable probabilities to values greater than the number of
allocated components. Similar observations have been by Richardson and Green
(1997), who specify a prior on the number of components $M$, highlighting the
fact that some of the components might be empty as not all the components
might be represented in a finite sample and the data are non-informative on
unallocated components. This leads to an identifiability problem for $M$ and,
as a consequence, fully non-informative priors cannot be elicited in a mixture
context. Nevertheless, Richardson and Green (1997) still focus their inference
problem on $M$ and do not investigate the relationship between $M$ and $k$.
More recently, Malsiner-Walli et al. (2016) introduce sparse finite mixture
models as an alternative to infinite (nonparametric) mixtures, and impose
sparsity to estimate the number of non-empty components in a deliberately
over-fitting mixture model where $M$ is fixed relatively large. On the other
hand, in Bayesian nonparametrics $M$ is set equal to infinity (i.e.,
$M=\infty$) and the focus of inference is only $k$. In this work we stress the
importance of the distinction between $M$ and $k$ as it will allow us to
collocate nonparametric and parametric mixtures in exactly the same framework.

In Bayesian Nonparametrics (i.e., $M=\infty$), the Dirichlet process mixture
model (DPM, Lo, 1984; Neal, 2000) – i.e. the Model in (1) where the mixing
measure is indeed the Dirichlet process – plays a pivotal role. DPM popularity
is mainly due to its high flexibility and mathematical and computational
tractability both in density estimation as well as in clustering problems.
However, in some statistical applications, the use of the Dirichlet process as
a clustering mechanism may be restrictive (see, for instance, Lau and Green,
2007; Miller and Harrison, 2013): the clustering results often depend on the
choice of a particular kernel, partitions will typically be dominated by few
large clusters (the rich-get-richer property), the number of clusters
increases as the number of observations $n$ increases (as $\log n$), often
leading to the creation of too many non-interpretable singleton clusters. To
overcome these drawbacks many alternative mixing measures have been proposed
(e.g. Ishwaran and James, 2003; Dey et al., 2012). In particular, Lijoi et al.
(2007) replace the Dirichlet process with a large and flexible class of random
probability measures obtained by _normalization_ of random (infinite
dimensional) measures with independent increments. Once again, all these
approaches assume $M=\infty$ and focus on estimating $k$.

On the other hand, in a Bayesian parametric context (i.e., $M<\infty$ almost
surely) the most popular approaches are (i) fix $M$ and then focus mainly on
density estimation (ii) treat $M$ as a random parameter and make it the focus
of inference. Then, conditionally on $M=m$, the mixture weights
$(w_{1},\dots,w_{m})$ are chosen according to a $m-1$ dimensional
$Dirichlet_{m}$ distribution. We refer to the latter model as _finite
Dirichlet mixture model_ (FDMM). Refer, among the others, to Nobile (1994);
Richardson and Green (1997); Stephens (2000) and Miller and Harrison (2018)
for more details. The literature is rich of proposals on how to estimate the
number of components $M$, but there is no consensus on the best method.
Likelihood based inference typically relies on model choice criteria, such as
BIC or the approximate weight of evidence (see Biernacki et al., 2000, for a
review). Although in the Bayesian paradigm there are approaches based on model
choice criteria, such as DIC, it is usually preferable to perform full
posterior inference on $M$ as well, eliciting an appropriate prior. A fully
Bayesian approach in FDMM is often based on the reversible jump Markov chain
Monte Carlo (Richardson and Green, 1997; Dellaportas and Papageorgiou, 2006)
or, alternatively, on the marginal likelihood $p(\bm{y}\mid M)$. Both methods
present significant computational challenges.

Although mainly for computational purposes, the connection between finite and
infinite mixture models has been present in the literature for at least two
decades since the work of Muliere and Tardella (1998). Since then extensive
research effort has been devoted to find approximate representation of the
Dirichlet process (e.g. Ishwaran and Zarepour, 2002). Moreover, algorithms for
posterior inference of infinite mixture models often truncate the infinite
measure and approximate it with a finite mixture with $L$ components, where
$L$ is sufficiently large or random (see for instance Ishwaran and James,
2001; Argiento et al., 2016) and, in practice, the inferential problem
translates into estimating the number of allocated components (clusters) and
the cluster-specific parameters. The focus of this work is to provide a
probabilistic treatment of mixture modelling, that reconciles the two
approaches: $M=\infty$ and $M<\infty$. This connection has received some
attention from a theoretical point of view (see e.g. Miller and Harrison,
2018; Frühwirth-Schnatter and Malsiner-Walli, 2018), but has never been
investigated thoroughly and fully resolved.

###  1.1 Contribution of this work

In this work, instead of approximating an infinite mixture with a parametric
one, we show that a finite mixture model is simply a realization of a
stochastic process whose dimension is random and has an infinite dimensional
support. To this end, we introduce a new class of random measures obtained by
normalization of a point process and use it as mixing measures in Model (1).
We refer to this new class as _Normalized Independent Finite Point Processes_
and we derive the family of prior distributions induced on the data partition
by providing a general formula for exchangeable partition probability
functions (Pitman, 1996). Finally, we characterize the posterior distribution
of the Normalized Independent Finite Point Process. Our construction is
exactly in the spirit of Bayesian nonparametrics, as it is based on the
normalization of a point process, leading to an almost surely discrete
measure. Indeed, there is a fundamental and simple idea behind the
construction of almost surely discrete random measures: they can be obtained
by normalization of stochastic processes. Already Ferguson (1973) in his
seminal work derived the Dirichlet process as normalization of a Gamma
process. More recently, Regazzini et al. (2003) propose a new class of
nonparametric priors, called Normalized Random Measures with Independent
Increments, obtained through the normalization of a Lévy process. This latter
work has opened the door to one of the most active lines of current research
in Bayesian statistics as well as in machine learning. On one hand it has led
to the development of nonparametric priors beyond the Dirichlet process (e.g.
Lijoi et al., 2010) and on the other the same techniques are widely used for
clustering in the machine learning community under the name of Normalised
Completely Random Measures  (e.g. Jordan, 2010).

The class we propose is rich and includes as a particular case the popular
finite Dirichlet mixture model. Several inference methods have been proposed
for the finite mixture models, of which the most commonly-used are the
Reversible Jump Markov chain Monte Carlo (Green, 1995; Richardson and Green,
1997) and the recently proposed marginal algorithm (Miller and Harrison,
2018). The Reversible jump algorithm is a very general technique and has been
successfully applied in many contexts, but it can be difficult to implement
since it requires designing problem-specific moves, which is often a
nontrivial task particularly in high-dimensions. The algorithm proposed by
Miller and Harrison (2018), although more efficient and in some ways more
automatic, restrains the class of prior distributions for the weights and does
not allow inference on the hyper-parameters of the process, which could
constitute a serious limitation in complex set-ups. In Miller and Harrison
(2018), by integrating out the mixing measure, inference is limited to the
number of allocated components and the mean of linear functionals of the
posterior distribution of the mixture model. See Gelfand and Kottas (2002) for
a discussion of these issues.

Among the main achievements of this work, there is the construction of a Gibbs
sampler scheme to simulate from the posterior distribution of the Normalized
Finite Independent Point Process, in particular a conditional MCMC algorithm
based on the posterior characterization of such process. This algorithm, in
the particular case of a Dirichlet prior on the mixture weights, leads to
conjugate updating with a substantial gain in computational efficiency over
current algorithms. The key result (associated to the nonparametric
construction of the process) is to be able to propose transdimensional moves
which are automatic and naturally implied by the prior process. We illustrate
the proposed prior process through the benchmark example offered by the the
Galaxy data (Roeder, 1990) and an important application in population
genetics.

The manuscript is structured as follows. In Section 2 we introduce the finite
mixture model framework, highlighting the connection between parametric and
nonparametric constructions. Section 3 reviews necessary theory from Finite
Point Processes. In Section 4 we introduce the prior process, the Normalised
Independent Finite Point Process, and discuss its clustering properties, while
in Section 5 we characterise its posterior distribution. In Section 6 we show
how the new construction leads to efficient conditional and marginal
algorithms. In Section 7 we briefly show how the new prior can be used as a
component in more complex hierarchies. In Section 8 we discuss possible
choices for the prior on the number of mixture components, while Section 9
presents special cases of the process. We conclude the paper with two
examples: (i) the Galaxy data example, which provides an opportunity to
benchmark our method in Section 10; (ii) a real data genetic application aimed
to identify population structure from microsatellites loci genotyped in a
sample of thrushes in Section 11. We conclude the paper in Section 12.

##  2 Finte Mixture Models

In this section we introduce the finite mixture model (FMM) and show how it
can be written in three equivalent ways. The first two representations are
widely used in the parametric literature, while the last one uses notions
typical of random mixing measure (nonparametric) set-ups. Our starting point
is the general finite mixture model. Let $Y_{1},\dots,Y_{n}$ be a set of
observations taking values in an Euclidean space $\mathcal{Y}$. We consider
the following mixture model:

$\begin{split}&Y_{i}|\tau_{1},\ldots,\tau_{M},\mathbf{w},M\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\sum_{m=1}^{M}w_{m}f(y\mid\mathbb{\tau}_{m})\\\
&\tau_{i}\mid
M\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}P_{0}(d\tau),\qquad
i=1,\ldots,M\\\ &\mathbf{w}\mid M\sim P_{W}(\mathbf{w})\\\ &M\sim
q_{M}\end{split}$ (2)

where $f(\cdot;\tau_{m})$ is a parametric density on $\mathcal{Y}$, which
depends on a vector of parameters $\tau_{m}$. The vector of parameters assume
values in $\Theta\subset{\mathbb{R}}^{d}$ and is assigned a non-atomic prior
density $p_{0}$ corresponding to the probability measure $P_{0}$ on $\Theta$.
The number of components is an important parameter of the mixture model and in
a fully Bayesian approach it is given a prior $q_{M}$. Conditionally on $M$,
the vector of weights $\mathbf{w}=(w_{1},\ldots,w_{M})$, which represents the
probability of belonging to each mixture component, is given a prior
probability $P_{W}$ on the simplex of dimension $M-1$. The model in Eq. (2)
can be rewritten in terms of latent variables, since this representation
allows for simpler computations. To this end, we need to introduce a latent
allocation vector $\mathbf{c}=(c_{1},\ldots,c_{n})$, whose element $c_{i}$
denotes to which component observation $Y_{i}$ is assigned,
$c_{i}\in\\{1,\ldots,M\\}$. Then the model in (2) is equivalent to:

$\begin{split}&Y_{i}|\theta_{i}\stackrel{{\scriptstyle\mathrm{ind}}}{{\sim}}f(y\mid\theta_{i}),\qquad
i=1,\dots,n\\\
&\theta_{i}|c_{i}\stackrel{{\scriptstyle\mathrm{ind}}}{{\sim}}\delta_{\tau_{c_{i}}}(d\theta_{i})\\\
&\tau_{m}\mid
M\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}P_{0}(d\tau),\qquad
m=1,\ldots,M\\\ &\mathbf{w}\mid M\sim P_{W}(\bm{w}\mid M)\\\ &c_{i}\mid
M,\mathbf{w}\sim\text{Multinomial}_{M}(1,w_{1},\ldots,w_{M})\\\ &M\sim
q_{M}\end{split}$ (3)

where $\delta_{\tau}$ is the Dirac measure assigning unit mass at location
$\tau$. Usually $P_{W}$ is assumed to be a
Dirichlet${}_{M}(\gamma,\ldots,\gamma)$ distribution, while typical choices
for $q_{M}$ include a discrete uniform on some finite space, a Negative
Binomial or a Poisson distribution. Note that prior information about the
relative sizes of the mixing weights $w_{1},\ldots,w_{M}$ can be introduced
through $\gamma$ – roughly speaking, small $\gamma$ favours lower entropy
$w$’s, while large $\gamma$ favours higher entropy $w$’s. In general, the
hyperparameter $\gamma$ is either set equal to a constant (e.g. equal 1 or
$1/M$) or is assigned a Gamma hyperprior. Frühwirth-Schnatter and Malsiner-
Walli (2018) propose a sparsity prior on $\mathbf{w}$, which allows the number
of non-empty components to be much smaller than $M$, where $M$ is a non-random
constant. Alternatively, Miller and Harrison (2018) do not make strong
assumptions on the prior on $M$ but, by showing the connection between FMM and
exchangeable partition probability functions (eppf), manage to apply the well-
developed inferential methods for DPMs to FMMs with significant gains in
computational efficiency. The strategy proposed by Miller and Harrison (2018)
is limited to the Dirichlet prior on $\mathbf{w}$ and employs a marginal-type
of algorithm to perform posterior inference. This approach, often used in
DPMs, marginalises over the weights of the mixture and it is most appropriate
when the main object of scientific interest is $k$. In this work we propose a
richer construction, where the prior on $\mathbf{w}$ is obtained normalising a
finite point process. Advantages of the proposed approach include: (i)
extension of the family of prior distributions for the weights; (ii) full
Bayesian inference on all the unknowns (in particular $M$ and $\mathbf{w}$);
(iii) possibility of inducing sparsity through appropriate choice of hyper-
parameters; (iv) ease of interpretation; (v) possibility of extending the
construction to covariate dependent weights and (vi) extension to more general
processes.

In a nutshell, we build a general class of finite mixture models by proposing
a new prior process for $P_{W},q_{M},P_{0}$ which admits the conventional
mixture model described in Eq. (3) as special case. To introduce this new
class of prior distributions, which we refer to as Normalized Independent
Finite Point Processes (Norm-IFPP), we first need to review some background
theory and introduce some notation. Then, in Section 9 we provide examples
which do not require the prior for the weights to be a Dirichlet distribution.

The theoretical developments are based on the key observation that a
realization $M,\bm{w},\bm{\tau}$ from the prior on the mixture model
parameters defined in Eq. (2) in terms of hierarchical parametric distribution
$q_{M},P_{W},P_{0}$ defines an almost surely (a.s.) finite-dimensional random
probability measure on the parameter space $\Theta$:

$P(\theta)=\sum_{m=1}^{M}w_{m}\delta_{\tau_{m}}(d\theta)$ (4)

This implies that the joint probability distribution on $M$, $\mathbf{w}$ and
$\bm{\tau}$ induces a distribution on $P$ defined in Eq. (4), whose support is
the space of the a.s. finite-dimensional random probability measures on
$\Theta$. Moreover, it is straightforward to prove (see Argiento et al., 2019)
that by letting $\theta_{i}=\tau_{c_{i}}$, as in Eq. (3), the variables
$\theta_{1},\dots,\theta_{n}$ can be considered as a sample from $P$, i.e.
$\theta_{1},\dots,\theta_{n}|P\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}P$.
From this observation, the link between infinite (nonparametric) and finite
mixture models becomes evident as the model in Eq. (2) can be easily rewritten
as

$\begin{split}&Y_{1},\dots,Y_{n}|\theta_{1},\dots,\theta_{n}\stackrel{{\scriptstyle
ind}}{{\sim}}f(y;\theta_{i})\\\
&\theta_{1},\dots,\theta_{n}|P\stackrel{{\scriptstyle iid}}{{\sim}}P\\\
&P\sim\mathcal{P}\end{split}$ (5)

where $P$ is defined in Eq. (4) and $\mathcal{P}$ is the law of $P$ defined
via $q_{M},P_{W},P_{0}$. The main theoretical contribution of this work is to
give a constructive definition of $\mathcal{P}$, introducing a class of FMM
for which the weights $w_{m}$ represent the normalised jumps of a _finite
point process_ and the parameters $\tau_{m}$ are defined in terms of
realisations of the same point process. As in any mixture, $\theta_{i}$s in
Model (5) are equal to one of the $\tau_{m}$ in Eq. (4), depending on which
component the $i$th observation is assigned to. The link between finite
mixture models and point processes is not unknown, as pointed out in the
introduction. In particular, Stephens (2000) highlights this connection, but
defines the point process on the complex space of normalized weights (i.e. the
union of infinite simplexes). In this work through normalization not only we
are able to work on a simpler space, but also to build a new general class of
distribution $P_{W}$, i.e. a new class of prior for the weights of the mixing
measure $\bm{w}$.

##  3 Finite Point Processes

In this section we review some concepts from point process theory which are
necessary to construct the Norm-IFPP process. We refer to the books of Daley
and Vere-Jones (2007) and Møller and Waagepetersen (2003) for a complete
treatment of finite point processes.

Let $\mathcal{X}$ be a complete separable metric space, a _finite point
process_ $X$ is a random countable subset of $\mathcal{X}$. In this paper we
restrict our attention to processes whose realizations are a finite subset of
$\mathcal{X}$. For any realization of the process, $x\subseteq\mathcal{X}$,
let $\\#(x)$ denote the cardinality of $x$. The realizations of $X$ are
constrained on $N_{f}=\\{x\subseteq\mathcal{X}:\\#(x)<\infty\\}.$ Elements of
$N_{f}$ are called _finite point configurations_. The law of a finite point
process is identified by the following quantities:

  1. 1.

a discrete probability density $q_{M},M=0,1,\dots$ determining the law of the
total number $M$ (i.e. $\\#(x)$) of points of the process,

  2. 2.

for each integer $M\geq 1$, a probability distribution $\Pi_{M}(\cdot)$ on the
Borel sets of $\mathcal{X}^{M}$, that determines the joint distribution of the
positions of the points of the process, given that their total number is $M$.

In particular, $q_{M}$ and $\Pi_{M}$ provide a constructive definition of the
process which is very useful in simulations. First generate a random integer
$M$ from $q_{M}$ and then, given $M\neq 0$, generate a random set
$X=\\{\xi_{1},\dots,\xi_{M}\\}$ which is a sample from $\Pi_{M}(\cdot)$. If
$M=0$, the random generation stops and $X$ coincides with the empty set.

Note that a point process $X=\\{\xi_{1},\dots,\xi_{M}\\}$ is a set of
unordered points. To this end, the distributions $\Pi_{M}(\cdot)$ needs to
give equal weight to all $M!$ permutations of the elements in the vector
$(\xi_{1},\dots,\xi_{M})$, i.e. $\Pi_{M}(\cdot)$ must be symmetric. A
convenient way to specify the law of $X$ is based on the _Janossy measure_
(Daley and Vere-Jones, 2007):

$J(A_{1}\times\dots\times A_{M})=M!q_{M}\Pi_{M}(A_{1}\times\dots\times A_{M})$

where the $A_{m}$s are Borel-sets of $\mathcal{X}$, with $\xi_{m}\in A_{m}$.
The _Janossy measure_ is unnormalised and plays a fundamental role in the
study of finite point processes and spatial point patterns. It has a simple
interpretation which makes it easy to work with. Let
$\mathcal{X}={\mathbb{R}}^{d}$ and let $j(\xi_{1},\dots,\xi_{M})$ denote the
density of $J(\cdot)$ with respect to the Lebesgue measure with
$\xi_{m}\neq\xi_{\tilde{m}}$ for $m\neq\tilde{m}$. Then
$j(\xi_{1},\dots,\xi_{M})d\xi_{1}\dots d\xi_{M}$ is the probability that there
are exactly $M$ points in the process, one in each of the distinct
infinitesimal regions $(\xi_{m},\xi_{m}+d\xi_{m})$. Here, we will use the
Janossy measure to characterize the prior and the posterior distribution of
the new class of finite discrete random probability measures, i.e. the family
of _Normalised Independent Finite Point Processes_ (Norm-IFPP). We now
introduce a simplified version of this process which assumes that the points
$\xi_{m}$ are conditionally independent and identically distributed.

######  Definition 1.

Let $\nu(\cdot)$ and $\\{q_{M},M=0,1,\dots\\}$ be a density on $\mathcal{X}$
and a probability mass function respectively. $X$ is an independent finite
point process, $X\sim IFPP(\nu,q_{M})$, if its Janossy density can be written
as

$j(\xi_{1},\dots,\xi_{M})=M!q_{M}\prod_{m=1}^{M}\nu(\xi_{m})$ (6)

In what follows, our construction is based on Eq. (6).

##  4 Normalized Independent Finite Point Processes

Let $\Theta\subset{\mathbb{R}}^{d}$, for some positive integer $d$ and let
$\mathcal{X}$ be ${\mathbb{R}}^{+}\times\Theta$. We denote with $\xi=(s,\tau)$
a point of $\mathcal{X}$. Let $\nu(s,\tau)$ be a density on $\mathcal{X}$ such
that $\nu(s,\tau)=h(s)p_{0}(\tau)$, where $h(\cdot)$ is a density on
${\mathbb{R}}^{+}$ and $p_{0}(\cdot)$ is a density on $\Theta$. Finally, we
consider only $q_{M}$ such that $q_{0}=0$, i.e. the prior probability of $M=0$
is zero. We consider the independent finite point process
$\tilde{P}=\\{(S_{1},\tau_{1}),\dots,(S_{M},\tau_{M})\\}$ with parameters
$\nu$ and $\\{q_{M}\\}$, i.e. $\tilde{P}\sim IFPP(\nu,q_{M})$. In what
follows, it is easier to introduce a slight change of notation and define
$IFPP(h,q_{m},p_{0})=IFPP(\nu,q_{M})$ to highlight the dependence of the
process also on $p_{0}(\cdot)$. Let $\mathcal{M}:=\\{1,\dots,M\\}$ be the set
of indexes corresponding to the points of the process. Since we are assuming
that $q_{0}=0$ the random variable $T:=\sum_{m\in\mathcal{M}}S_{m}$ is almost
surely larger than $0$ so that we can give the following definition:

######  Definition 2.

Let $\tilde{P}=\\{(S_{1},\tau_{1}),\dots,(S_{M},\tau_{M})\\}\sim
IFPP(h,q_{m},p_{0})$, with $q_{0}=0$. A normalized independent finite point
process (Norm-IFPP) with parameters $h(\cdot),\ p_{0}$ and $\\{q_{M}\\}$ is a
discrete probability measure on $\Theta$ defined by

$\begin{split}P(A)&=\sum_{m\in\mathcal{M}}w_{m}\delta_{\tau_{m}}(A)\stackrel{{\scriptstyle
d}}{{=}}\sum_{m\in\mathcal{M}}\frac{S_{m}}{T}\delta_{\tau_{m}}(A)\end{split}$
(7)

where $T=\sum_{m\in\mathcal{M}}S_{m}$ and $A$ denotes a measurable set of
$\Theta$. We refer to the process in Eq. (7) as $P\sim{\rm Norm-
IFPP}(h,q_{m},p_{0})$.

The finite dimensional process defined in Eq. (7) belongs to the wide class of
species sampling models (see Pitman, 1996) and this will allow us to use all
the efficient machinery developed for such models. Let
$\left(\theta_{1},\dots,\theta_{n}\right)$ be a sample from a Norm-IFPP. It is
well known that sampling from a discrete probability measure induces ties
among the $\theta_{i}$s and, therefore, a random partition of the
observations. Let ${\rho}_{n}:=\\{C_{1},\dots,C_{k}\\}$ indicate a partition
of the set $\\{1,\dots,n\\}$ in $k$ subsets, where
$C_{j}=\\{i:\theta_{i}=\theta^{\star}_{j}\\}$ for $j=1,\dots,k\leq n$, and let
$\\{\theta_{1}^{\star},\ldots,\theta_{k}^{\star}\\}$ denote the set of
distinct $\theta_{i}$s associated to each $C_{i}$. The marginal law of
$\left(\theta_{1},\dots,\theta_{n}\right)$ has a unique characterization:

$\mathcal{L}\left(d\theta_{1},\dots,d\theta_{n}\right)=\mathcal{L}({\rho}_{n},d\theta^{\star}_{1},\dots,d\theta^{\star}_{k})=\pi(n_{1},\dots,n_{k})\prod_{j=1}^{k}P_{0}(d\theta^{\star}_{j})$

where $n_{j}=\\#(C_{j})$, $\sum_{j=1}^{k}n_{j}=n$ and $\pi(\cdot)$ is the
exchangeable partition probability function (eppf) associated to the random
probability $P$ (Pitman, 1996). For each $n$, the eppf $\pi$ is a probability
law on the set of the partitions of $\\{1,\dots,n\\}$, which determines the
(random) number of clusters $k$ and the numerosity of each cluster $C_{i}$.
The partition is exchangeable because its law depends only on the number and
size of the clusters, and not on the allocation of the individuals to each
clusters. The eppf is a key tool in Bayesian analysis as mixture models can be
rewritten in terms of random partitions and such equivalence is often
exploited to improve computational efficiency, in particular of marginal
algorithms (Lijoi et al., 2010). The following proposition provides an
expression for the eppf of a Norm-IFPP measure.

######  Theorem 1.

Let $(n_{1},\dots,n_{k})$ be a vector of positive integers such that
$\sum_{j=1}^{k}n_{j}=n$. Then, the eppf associated with a Norm-
IFPP$(h,p_{0},q_{M})$ is

$\pi(n_{1},\dots,n_{k})=\int_{0}^{+\infty}\dfrac{u^{n-1}}{\Gamma(n)}\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}\prod_{j=1}^{k}\kappa(n_{j},u)du$
(8)

where $\psi(u)$ is the Laplace transform of the density $h(s)$, i.e.

$\psi(u):=\int_{0}^{\infty}e^{-us}h(s)ds$ (9)

and

$\kappa(n_{j},u):=\int_{0}^{\infty}s^{n_{j}}e^{-us}h(s)ds=(-1)^{n_{j}}\frac{d}{du^{n_{j}}}\psi(u)$

Proof: See Appendix $\blacksquare$

For what follows it is important to highlight the difference between $M$ and
$k$. The number of components of the finite mixture $M$ is given by a
realisation of the process in Eq. (7). On the other hand, $k$ denotes the
number of non-empty (allocated) components, with $k\leq M$. This difference
has been noted before in the literature (see, for example, Nobile et al.
(2004); Miller and Harrison (2018); Frühwirth-Schnatter and Malsiner-Walli
(2018)). Suppose that a realization from $\mathcal{P}$ is a discrete measure
with $M=4$ atoms in Eq. (4) and $\bm{\tau}=(0.2,2.4,1.5,4.1)$. Furthermore we
have a realization from $P$, $\bm{\theta}=(0.2,0.2,4.1,2.4,0.2,2.4)$, with
$n=6$. Then the allocated components are $k=3$ and the total number of mixture
components is $M=4$. Note that the representation in Eq. (7) implies that the
jumps of the point process, indexed by the elements of
$\mathcal{M}=\\{1,\dots,M\\}$, correspond to the components of the finite
mixture, and their relative size defines the weights. More formally, we denote
by $\mathcal{M}^{(a)}$ the set of indexes of _allocated_ jumps of the Process
(7), i.e. the indexes $m\in\mathcal{M}$ corresponding to some jumps $S_{m}$
such that there exists a location for which $\tau_{m}=\theta^{\star}_{i}$,
$i=1,\dots,k$. The remaining values of $\mathcal{M}$ correspond to the _non-
allocated_ jumps and we denote this set with $\mathcal{M}^{(na)}$. We use the
superscripts $(a)$ and $(na)$ for random variables related to _allocated_ and
_non-allocated_ jumps respectively.

One of the main focus of inference when using finite mixture models is to
determine the clustering allocation of the observations. The eppf gives the
prior distribution on the space of possible partitions. Moreover,
marginalising over the cluster sizes, it is also possible to derive the
implied prior distribution on the number of clusters, $k$, which corresponds
to the number of allocated components.

######  Corollary 1.

Under the assumptions of Theorem 1, the marginal prior probability of sampling
a partition with $k$ clusters is given by

$p^{\star}_{k}=\Pr\\{M^{(a)}=k\\}=\int_{0}^{+\infty}\dfrac{u^{n-1}}{\Gamma(n)}\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}B_{n,k}(\kappa(\cdot,u))$
(10)

where $k=1,\ldots,n$, and $B_{n,k}(\kappa(\cdot,u))$ is the partial Bell
polynomial (Pitman, 2006) over the sequence of coefficients $\\{\kappa(n,u),\
n=1,2,\dots\\}$.

Proof: See Appendix A.2 $\blacksquare$

Moreover, from de Finetti’s theorem it follows that, $k$ converges almost
surely to $M$, as $n\rightarrow\infty$.

##  5 Posterior Carachterization of a Norm-IFPP Process

In this section we characterise the posterior distribution of the process
$P\sim\text{Norm-IFPP}(h,p_{0},q_{M})$. To this end, we introduce the random
variable $U_{n}=\Gamma_{n}/T$, where $\Gamma_{n}\sim\text{Gamma}(n,1)$, with
$\Gamma_{n}$ and $T$ independent, where $T=\sum_{i\in\mathcal{M}}S_{i}$. It is
easy to show (see the Appendix A.4) that if $P\sim\text{Norm-
IFPP}(h,p_{0},q_{M})$ then, for any $n\geq 1$, the marginal density of $U_{n}$
is given by

$f_{U_{n}}(u;n)=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\frac{d}{du^{n}}\mathbb{E}\left(\psi(u)^{M}\right)\\\
$ (11)

where $\psi(u)$ is the Laplace transform of the density $h$, as defined in Eq.
(9). We give the derivation of Eq. (11) in Appendix A.4. The posterior
distribution of $U_{n}$, given $\bm{\theta}=(\theta_{1},\ldots,\theta_{n})$,
is crucial to perform posterior inference and allows us to derive the
posterior distribution of the unnormalised process $\widetilde{P}$. To this
end, we need to show that a posteriori, conditionally to $U_{n}$,
$\widetilde{P}$ is the superposition (union) of two independent process: a
point process and a finite process with fixed locations at
$(\theta^{\star}_{1},\ldots,\theta^{\star}_{k})$. Note that $k$ corresponds to
the number of allocated jumps $M^{(a)}$ and $M$ is equal to the sum of $k$ and
the number $M^{(na)}$ of unallocated jumps, assuming values in
$\mathbb{N}\cup\\{0\\}$. The process of unallocated jumps is a latent variable
which links the parametric part of the model in $P$ to a nonparametric
process. This link is essential for computations as it will become clearer in
Section 6, where we discuss the algorithm. The results below are conditional
on the realizations of the random variable $U_{n}$, which is a typical
strategy in the theory of normalised random measures, since working on the
augmented space allows us to exploit the quasi-conjugacy of the process $P$
(see James et al., 2009). We now present the main theoretical contribution of
this work.

######  Theorem 2.

If $P\sim\text{Norm-IFPP}(h,p_{0},q_{M})$, then the unnormalized process
$\widetilde{P}$, given
${\bm{\theta}}^{\star}=(\theta^{\star}_{1},\dots,\theta^{\star}_{k})$,
$\bm{n}=(n_{1},\dots,n_{k})$ and $U_{n}=u$, is the superposition of two
processes:

$\widetilde{P}\stackrel{{\scriptstyle
d}}{{=}}\widetilde{P}^{(na)}\cup\widetilde{P}^{(a)}$

where

  1. 1.

The process of non-allocated jumps $\widetilde{P}^{(na)}$ is an independent
finite point process with Janossy density given by

$j_{m}((s_{1},\tau_{1}),\dots,(s_{m},\tau_{m}))=m!q_{m}^{\star}\prod_{j=1}^{m}h^{\star}(s_{j})p_{0}(\tau_{j})$

where $h^{\star}_{u}(s)\propto\mathrm{e}^{-us}h(s)$,
$q^{\star}_{m}\propto\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}$, $\psi(u)$ is the
Laplace transform of $h$, and $m$ is a realization of $M^{(na)}$, the number
of unallocated jumps, taking values in $\\{0,1,2,\dots\\}$.

  2. 2.

The process of allocated jumps $\widetilde{P}^{(a)}$ is the unordered set of
points $(S_{1},\tau_{1}),\dots,(S_{k},\tau_{k})$, such that, for
$j=1,\dots,k$, $\tau_{j}=\theta^{\star}_{j}$ and the distribution of $S_{j}$
is proportional to $s^{n_{j}}\mathrm{e}^{-us}h(s)$.

  3. 3.

Conditionally on $\mathcal{M}^{(a)}$ and $U_{n}=u$, $\widetilde{P}^{(a)}$ and
$\widetilde{P}^{(na)}$ are independent.

Moreover, the posterior law of $U_{n}$ given
${\bm{\theta}}=(\theta_{1},\ldots,\theta_{n})$ depends only on the partition
$\rho_{n}$ and has density on the positive reals given by

$f_{U}(u\mid\rho_{n})\propto\dfrac{u^{n-1}}{\Gamma(n)}\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}\prod_{i=1}^{k}\kappa(n_{j},u)$

Proof: See Appendix A.3 $\blacksquare$

The result in Theorem 2 is the finite dimensional counterpart of Theorem 1 in
James et al. (2009) for normalised completely random measure. This theorem
will allow building an efficient block Gibbs sampler for finite mixture
models. Since the order in which the points of a point process arise is not
important, without loss of generality, given a realization of the posterior
process $\widetilde{P}$, we assume that, in
$\widetilde{P}=\\{S_{1},\dots,S_{M}\\}$, $M=k+M^{(na)}$, i.e. the first $k$
points $\\{S_{1},\dots,S_{k}\\}$ correspond to the allocated jumps, while the
last $M^{(na)}$ to the non-allocated ones.

##  6 Posterior inference

To perform posterior inference tailored MCMC algorithms need to be devised.
The two most popular strategies in Bayesian nonparametrics are marginal (Neal,
2000) and conditional algorithms (Ishwaran and James, 2001; Kalli et al.,
2011; Argiento et al., 2016). Our construction allows for straightforward
extension of such startegies to the finite mixture case, offering a convenient
alternative to the often inefficient and labour intensive reversible jump. To
implement marginal algorithms it is desirable (although not necessary, but at
the cost of extra computations) to be able to compute the sum in Eq. (8) to
obtain the probability of a random partition. On the other hand, for
conditional algorithms we need to sample from the posterior distribution of a
Norm-IFPP which requires a closed form expression for the Laplace transform in
Theorem 2. More specifically, it is essential to be able to sample from the
posterior distribution of the number of the non-allocated jumps,
$q^{\star}_{m}$, as well as from the distribution of the allocated and
unallocated jumps, i.e. the densities proportional to $e^{-us}h(s)$
(Exponential tilted) and $s^{n_{j}}e^{-us}h(s)$ (Gamma tilted). Specific
solutions for well known processes will be presented in the following
sections. Here we give a general outline of both algorithms.

###  6.1 Marginal Algorithm

As mentioned before, a sample $\theta_{1},\dots,\theta_{n}$ from $P$ induces a
partition of the set of the data indexes, denoted by
$\rho_{n}=\\{C_{1},\dots,C_{k}\\}$, such that $i\in C_{j}$ implies that datum
$i$ belongs to cluster $j$. Marginal algorithms rely on the fact that, by
integrating out the measure $P$, the only parameters left in Eq. (5) are the
random partition $\rho_{n}$ and the cluster specific parameters
$\theta_{1}^{\star},\dots,\theta_{k}^{\star}$. Posterior sampling strategies
for $\rho_{n}$ are based on the Chinese restaurant process (Aldous, 1985),
which describes the (a priori) predictive generative process for $\rho_{n}$,
and relies on the evaluation of the eppf associated with $P$. Nevertheless,
when $P$ corresponds to the Norm-IFPP model, this evaluation can be
computationally burdensome due to the integral with respect to $u$ in Eq. (8).
To design efficient algorithms we adopt a _disintegration_ technique following
a strategy similar to the one suggested by James et al. (2009) and Favaro and
Teh (2013) for NRMI. In particular, we augment the state space introducing the
latent variable $U$ (see Theorem 2).

We now explain how, conditional to the latent variable $U_{n}:=U$, the Chinese
restaurant process can be adapted to this set-up. Recall that the marginal
distribution of $U_{n}$, defined in Section 7, with $U_{n}\mid
T\sim\text{Gamma}(n,T)$, has been derived in Eq. (11):
$f_{U_{n}}(u;n)=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\frac{d}{du^{n}}\mathbb{E}\left(\psi(u)^{M}\right)$.

The partition (or clustering) $\rho_{n}$ can be generated using the eppf
derived in Theorem 1. It is straightforward to show that

$\begin{split}\mathcal{L}(\rho_{n},U_{n})=\pi(n_{1},\dots,n_{k};u)&=\dfrac{u^{n-1}}{\Gamma(n)}\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}\prod_{j=1}^{k}\kappa(n_{j},u)\\\
&=\dfrac{u^{n-1}}{\Gamma(n)}\Psi(u,k)\prod_{j=1}^{k}\kappa(n_{j},u)\end{split}$
(12)

where
$\Psi(u,k):=\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}$.
This joint distribution allows us to derive the predictive probability
(conditionally on $U_{n}=u$) that observation $n+1$ belongs to a new cluster
$C_{k+1}$ is

$\displaystyle\mathbb{P}(n+1\in C_{k+1}|u,\rho_{n})$
$\displaystyle\propto\frac{\pi(n_{1},\dots,n_{k},1;u)}{\pi(n_{1},\dots,n_{k};u)}=\frac{S(u,k+1)}{S(u,k)}\kappa(1,u)$
(13)

while the predictive probability of belonging to an existing cluster is

$\displaystyle\mathbb{P}(n+1\in C_{j}|u,\rho_{n})$
$\displaystyle\propto\frac{\pi(n_{1},\dots,n_{j}+1,\dots,n_{k},1;u)}{\pi(n_{1},\dots,n_{j},\dots,n_{k};u)}=\frac{\kappa(n_{j}+1,u)}{\kappa(n_{j},u)},\quad
j=1,\dots,k$ (14)

As in a standard Chinese restaurant process (Aldous, 1985), a sequence of
customers (data $i=1,2,\dots$) enter a restaurant with an infinite number of
tables (groups $C_{1},C_{2},\dots$). The first customer sits at the first
table and a random variable $U_{1}$ is drawn. Then each subsequent customer
joins a new table with probability proportional to Eq. (13), or an existing
table with $n_{j}$ customers with probability proportional to Eq. (14). For
each new customer $i$, a variable $U_{i}$ is drawn. After $n$ customers have
entered the restaurant, the seating arrangement of customers around tables
corresponds to a partition $\rho_{n}$ of $\\{1,\ldots,n\\}$ with numerosity
$(n_{1},\ldots,n_{k})$, $n_{j}=\\#C_{j}$. The seating arrangement of the
customers is exchangeable, in the sense that any seating that leads to the
same number of occupied tables and the same number of customers per table has
the same probability. The main difference with the standard Chinese process
consists in updating the cluster allocation conditional on $U_{i}$’s. The
strategy of conditioning on a sequence of auxiliary variables to generalise
the Chinese restaurant process was introduced for infinite dimensional
measures by James et al. (2009). Here, we have derived the finite dimensional
counterpart.

A general scheme to implement a posterior Gibbs sampler for Norm-IFPP mixture
model is the following:

  * i.

Draw $\rho_{n}$ from $\mathcal{L}(\rho_{n}\mid\text{rest})$. This can be done,
for instance, using one of the several algorithms presented in Neal (2000), by
simply substituting the predictive distributions of the Dirichlet process with
the conditional predictive structure of a Norm-IFPP given in Eq.s (13) and
(14).

  * ii.

Draw $U_{n}$ from $\mathcal{L}(U_{n}\mid\text{rest})$. This update requires a
Metropolis step (or any other alternative that ensures that the chain is
invariant) with target distribution proportional to
$\mathcal{L}(\rho_{n},U_{n})$ in Eq. (12)

  * iii.

Draw $\theta_{j}^{\star}$, for each $j=1,\dots,k$, from
$\mathcal{L}(\theta^{\star}_{j}\mid\text{rest})$. In general, this is
straightforward and involves a simple parametric update from

$\prod_{i\in C_{j}}f(y_{i}\mid\theta_{j}^{\star})p_{0}(\theta_{j}^{\star})$

Special cases in which the full conditional distributions of $U_{n}$ and
$\rho$ have a simple expression will be discussed later.

###  6.2 Conditional Algorithm

Conditional algorithms are usually of wide applicability. The most famous
example of this type of strategy is the one proposed by Ishwaran and James
(2001), which consists of a blocked Gibbs sampler based on the stick-breaking
representation of a discrete random measure. Conditional algorithms allow us
to draw from the joint distribution of
$(M,\mathbf{\tau},\mathbf{S},\mathbf{c})$ in Eq. (3), where $w_{i}=S_{i}/T$,
which in turns defines a draw of the random probability measure on $\Theta$:

$P(d\theta)=\sum_{m=1}^{M}w_{m}\delta_{\tau_{m}}(d\theta)$

As the algorithm samples from the posterior distribution of the random
measure, we are able to perform full posterior inference, at least
numerically, on any functional of such distribution. These issues are
discussed in detail in Gelfand and Kottas (2002). Moreover, it is simple to
make inference on the hyper-parameters of the distributions of $M$ and
$\mathbf{S}$. An outline of the MCMC algorithm is given in Figure 1. The
scheme follows directly from Theorem 2, adapted to the mixture case. Note that
in step 2 of the algorithm, the relabelling of the mixture components is
essential so that the non-empty components correspond to the first $k$
components.

Repeat for g in 1…G: 1. Sample ${u}^{(g)}$ from a Gamma$(n,T)$ 2. For i=1,..,n
sample ${c_{i}}^{(g)}$ from a discrete distribution s.t.
$\mathbb{P}(c_{i}=j\mid\text{rest})\propto S_{j}f(y_{i}\mid\tau_{j}),\quad
j=1,\dots,M$ After resampling the vector $\mathbf{c}^{(g)}$, calculate the
number $k^{(g)}$ of unique values of $\mathbf{c}^{(g)}$ and relabel the
mixture components in a way that the first $k^{(g)}$ ones are allocated. 3.a
Sample the hyperparametrs $\eta_{1}^{(g)}$ of the density $h$ from
$\mathbb{P}(\eta_{1}=d\eta_{1}\mid\text{rest})\propto\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}\prod_{j=1}^{k}\kappa(n_{j},u)\pi_{1}(\eta_{1})d\eta_{1}$
where $\pi_{1}(\eta_{1})$ denotes the prior density for $\eta_{1}$. 3.b Sample
the hyperparametrs $\eta_{2}^{(g)}$ of the density $q_{M}$ from
$\mathbb{P}(\eta_{2}=d\eta_{2}\mid\text{rest})\propto\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}\prod_{j=1}^{k}\kappa(n_{j},u)\pi_{2}(\eta_{2})d\eta_{2}$
where $\pi_{2}(\eta_{2})$ denotes the prior density for $\eta_{2}$. 4.a Sample
$M^{(na)(g)}$ from $q^{\star}_{m}\propto\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}$
and set $M^{(g)}=k^{(g)}+M^{(na)(g)}$ 4.b.$i$ Allocated jumps: for
$m=1,\ldots,k^{(g)}$, sample $S_{m}^{(g)}$ independently from
$\mathbb{P}(S_{m}=ds\mid\text{rest})\propto$ $s^{n_{m}}\mathrm{e}^{-su}h(s)ds$
4.b.$ii$ Non-allocated jumps: for $m=k^{(g)}+1,\ldots,M^{(g)}$, sample
$S_{m}^{(g)}$ independently from $\mathbb{P}(S_{m}=ds\mid\text{rest})\propto
e^{-us}h(s)ds$ 4.c.$.i$ Allocated points of support: sample $\tau_{m}^{(g)}$
independently from
$\mathbb{P}(\tau_{m}=d\tau_{m}\mid\text{rest})\propto\left\\{\prod_{i\in
C_{m}}f(y_{i}\mid\tau_{m})\right\\}p_{0}(\tau_{m})d\tau_{m}$ 4.c$.ii$ Non-
allocated points of support: sample $\tau_{m}^{(g)}$ independently from the
prior, i.e.
$\mathbb{P}(\tau_{m}=d\tau_{m}\mid\text{rest})=p_{0}(\tau_{m})d\tau_{m}$
Figure 1: Blocked Gibbs sampler scheme; the conditioning arguments of all full
conditionals have been omitted to simplify notation.

##  7 Norm-IFFP hierarchical mixture models

Most real world applications of discrete random measures involve an additional
layer in the model hierarchy and convolve the random measure with a continuous
kernel leading to nonparametric mixture models. In this context, data are
assumed to be generated from a parametric distribution indexed by some
parameter $\theta$, with $\theta\sim P$. Usually $P$ is assigned a
nonparametric prior, in our case a Norm-IFFP. This leads to models of the form

$\begin{split}&Y_{1},\dots,Y_{n}|\theta_{1},\dots,\theta_{n}\stackrel{{\scriptstyle
ind}}{{\sim}}f(y\mid\theta_{i})\\\
&\theta_{1},\dots,\theta_{n}|P\stackrel{{\scriptstyle iid}}{{\sim}}P\\\ &P\sim
Norm-IFFP(h,p_{0},q_{m})\end{split}$ (15)

where $f(\cdot\mid\theta_{i})$ is a parametric density on $\mathcal{Y}$, for
all $\theta\in\Theta\subset{\mathbb{R}}^{d}$. We point out that $p_{0}$ is the
density of a non-atomic probability measure $P_{0}$ on $\Theta$, such that
$\mathbb{E}(P(A))=P_{0}(A)$ for all $A\in\mathcal{B}(\Theta)$. Model (15) will
be addressed here as a Norm-IFFP hierarchical mixture model. The model can be
extended by specifying appropriate hyperpriors. It is well known that this
model is equivalent to assuming that the $Y_{i}$’s, conditional on $P$, are
independently distributed according to the random density (1). We point out
that Model (15) admits as a special case the popular finite Dirichlet mixture
model (see Nobile (1994); Richardson and Green (1997); Stephens (2000); Miller
and Harrison (2018)) discussed in more details in Section 9.1.1. The posterior
characterization given in Theorem 2, as well as the analytical expression for
the eppf given in Theorem 1, allow us device conditional or marginal
algorithms to perform inference under Model (15) as discussed in Section 6.

##  8 Special Choices of $q$ in Norm-IFFP

The exact evaluation of the eppf in Eq. (8) presents two challenges: an
integral and an infinite sum. Numerical solution of the integral is handled
within the MCMC via the augmentation trick, while here we discuss more in
detail the infinite sum, defined as

$\Psi(u,k):=\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}$

for each real $u>0$ and each integer $k\geq 1$. As it is shown in the proof of
Theorem 2, $\Psi(u,k)\leq\frac{k!}{(1-\psi(u))^{k+1}}$, i.e. the sum always
converges.

The analytical solution of the latter depends on the particular choice of
prior distribution $q$ for $M$. Since $\psi(u)$ is less than 1, $\Psi(u,k)$ is
related to a binomial series. This implies that if $q_{m}$ is a Poisson or a
Negative Binomial, we can derive conjugate updates for $M^{(na)}$ and we can
find a closed form solution for $\Psi(u,k)$. In particular, if
$q_{m}=\mathcal{P}_{1}(m,\Lambda)$, corresponding to the density of a random
variable shifted on $\\{1,2,\dots,\\}$, then we obtain

$\Psi(u,k)=\Lambda^{k-1}(\Lambda\psi(u)+k)\exp\\{\Lambda(\psi(u)-1)\\}$

Moreover, the full conditional distribution of $M^{(na)}$, i.e.
$q^{\star}_{m}$ in item (a) of Theorem 2, is

$q^{\star}_{m}\propto
k\psi(u)\mathcal{P}_{0}(m,\Lambda\psi(u))+\Lambda\psi(u)\mathcal{P}_{1}(m,\Lambda\psi(u))$

where $\mathcal{P}_{0}$ and $\mathcal{P}_{1}$ are the probability mass
function of a Poisson and of a shifted Poisson respectively.

Finally, it is worth to mention that the shifted Poisson choice for $M$
implies that in Eq. (11), we have

$\displaystyle
f_{U_{n}}(u;n)=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\frac{d}{du^{n}}\psi(u)\,{{\mathrm{e}}^{\Lambda\,\left(\psi(u)-1\right)}}$

Note that it is also possible to use a Truncated Poisson distribution for $M$,
with a slight difference in results.

On the other hand if we choose $q_{m}=\text{NegBin}(m;p,r)$, a Negative
Binomial density with parameters $0\leq p\leq 1$ and $r>0$ and support on
$\\{1,2,\dots\\}$, i.e.

$q_{m}=\frac{\Gamma(r+m-1)}{(m-1)!\Gamma(r)}p^{m-1}(1-p)^{r}\mathbb{I}_{\\{1,2,\dots\\}}(m)$

then, it is easy to show that

$\Psi(u,k)=\frac{\Gamma\left(r+k-1\right)}{\Gamma\left(r\right)}{p}^{k-1}\left(1-p\right)^{r}\frac{p\psi(u)(r-1)+k}{\left(1-p\,\psi(u)\right)^{k+r}}\quad
u>0\text{ and }\quad k\geq 1$

In this case we obtain that the full conditional for the number of non-
allocated components has support in $m=0,1,\dots$ with probability mass
function

$q^{\star}_{m}\propto(r+k)p\,\psi(u)\text{NegBin}(m;p\,\psi(u),r+k)+k\left(1-p\,\psi(u)\right)\text{NegBin}(m+1;p\,\psi(u),r+k-1)$

Moreover,

$\displaystyle f_{U_{n}}(u;n)$
$\displaystyle=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\frac{d}{du^{n}}{\frac{\psi(u)\,\left(1-p\right)^{r}}{\left(1-p\,\psi(u)\right)^{r}}}$

Finally, in applications, we might want to fix the number of mixture
components, i.e. the number of points of the point process, leading to the
standard finite mixture setup. In this case, if $M$ is set very large, we
recover the sparse mixture framework of Frühwirth-Schnatter and Malsiner-Walli
(2018). Let $M=\widetilde{M}\geq 1$ with probability 1, we obtain

$\Psi(u,k)=\left\\{\begin{array}[h]{ll}\frac{\widetilde{M}!}{(\widetilde{M}-k)!}\psi(u)^{\widetilde{M}-k}&\text{if
}k\leq\widetilde{M}\\\ 0&\text{if }k>\widetilde{M}\end{array}\right.$

This prior specification for $M$ implies that the support for $k$ is bounded,
$k\leq\widetilde{M}$, $q^{\star}_{m}$ assigns probability mass one to
$\widetilde{M}-k$ and

$\displaystyle
f_{U_{n}}(u)=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\widetilde{M}\psi(u)^{\widetilde{M}-1}\psi^{\prime}(u)$

##  9 Important Examples

The ${\rm Norm-IFPP}(h,p_{0},q_{M})$ depends on three densities. The prior on
$\theta_{i}$ is $p_{0}$ and, in applications, a conjugate prior is usually
preferred. The choice of $q_{M}$ has been discussed in Section 8. We now focus
on the choice of $h$. Once again the particular choice of $h$ influences the
induced clustering in Eq. (8) as well as efficiency of computations. There are
two possible alternatives: either to choose as $h$ a parametric density or to
select the Laplace transform of $h$, $\psi(u)$.

###  9.1 Choice of $h$

####  9.1.1 Finite Dirichlet process

Let $h$ be the $\text{Gamma}(\gamma,1)$ density. Under this choice of $h$ the
Norm-IFPP is a finite Dirichlet measure, that is an almost surely discrete
probability measure as in Eq. (7), where, conditionally on $M>0$, the jump
sizes $(w_{1},\dots,w_{M})$ of $P$ are a sample from the $M$-dimensional
Dirichlet${}_{M}(\gamma,\dots,\gamma)$ distribution. Therefore, this is
equivalent to a conventional finite mixture model as described in Section 2.
Recall that the Laplace transform and its derivatives for a
$\text{Gamma}(\gamma,1)$ density are given by
$\psi(u)=\frac{1}{(u+1)^{\gamma}},\ \text{and}\
\kappa(n_{j},u)=\frac{1}{(u+1)^{n_{j}+\gamma}}\frac{\Gamma(\gamma+n_{j})}{\Gamma(\gamma)},\quad
u>0,\ n_{j}=1,2,\dots$ Then, applying Theorem 1, we obtain that the eppf of
this model is

$\displaystyle p(n_{1},\dots,n_{k})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{u^{n-1}}{\Gamma(n)}\left(\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\frac{1}{(u+1)^{m\gamma}}q_{m+k}\right)\prod_{j=1}^{k}\frac{\Gamma(\gamma+n_{j})}{\Gamma(\gamma)}\frac{1}{(u+1)^{n_{j}+\gamma}}du$
$\displaystyle=$
$\displaystyle\left\\{\frac{1}{\Gamma(n)}\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}q_{m+k}\int_{0}^{\infty}\frac{u^{n-1}}{(u+1)^{m\gamma+n+kn}}du\right\\}\prod_{j=1}^{k}\frac{\Gamma(\gamma+n_{j})}{\Gamma(\gamma)}$
$\displaystyle=$
$\displaystyle\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}q_{m+k}\frac{\Gamma((k+m)\gamma)}{\Gamma((k+m)\gamma+n)}\right\\}\prod_{j}^{k}\frac{\Gamma(\gamma+n_{j})}{\Gamma(\gamma)}$
$\displaystyle=$ $\displaystyle
V(n,k)\prod_{j=1}^{k}\frac{\Gamma(\gamma+n_{j})}{\Gamma(\gamma)}$

See also Chapter 2 in Pitman (2006) and Miller and Harrison (2018). In
particular, if we choose _semi-conjugate_ priors for the number of components
as we discussed in Section 8, we can obtain integral representations for
$V(n,k)$. When $q_{m}$ is shifted Poisson distribution with parameter
$\Lambda$,

$V(n,k)=\Lambda^{k-1}\int_{0}^{\infty}\frac{u^{n-1}}{\Gamma(n)}\frac{\Lambda+k(u+1)^{\gamma}}{(u+1)^{n+\gamma(k+1)}}\exp\left\\{-\Lambda\frac{(u+1)^{\gamma}-1}{(u+1)^{\gamma}}\right\\}du$

while when $q_{m}$ is a Negative Binomial with with parameters $p$ and $r$

$V(n,k)=\frac{\Gamma(r+k-1)}{\Gamma(r)}{p^{k-1}(1-p)^{r}}\int_{0}^{\infty}\frac{u^{n-1}}{\Gamma(n)}\frac{p(r-1)+k(u+1)^{\gamma}}{(u+1)^{n+\gamma(r-1)}}\frac{1}{\left((u+1)^{\gamma}-p\right)^{r+k}}du$

Finally when $q_{m}$ assigns probability one to $\widetilde{M}$

$V(n,k)=\frac{\widetilde{M}!}{(\widetilde{M}-k)!}\frac{\Gamma(\gamma\widetilde{M})}{\Gamma(n+\gamma\widetilde{M})}\mathbb{I}_{\\{1,\dots,\widetilde{M}\\}}(k)$

The finite Dirichlet process case has been extensively discussed by Miller and
Harrison (2018). They propose a marginal algorithm which requires evaluating
the sum in Eq. (9.1.1). This restriction implies that an approximation of the
infinite sum needs to be evaluated at every step of the algorithm, slowing
down computations for large $n$ and making it difficult to specify a prior on
$\gamma$ and on the number of allocated components. We note that these
difficulties are easily overcome, for convenient choices of $q_{m}$, by
employing the disintegration trick and implementing a conditional MCMC scheme,
as described in Section 6.2. We highlight that when $q_{m}$ is a point mass,
the marginal algorithm becomes even more straightforward as we can obtain a
closed form expression for $V(n,k)$.

Note that, when $q_{m}$ is a shifted Poisson, if $\gamma=\alpha/\Lambda$, for
$\alpha>0$, and $\Lambda\rightarrow\infty$, then $P$ converges in distribution
to the Dirichlet process with mass parameter $\alpha$ (see Appendix A.6 for a
proof). Similarly, we recover the Dirichlet process when $q_{m}$ assigns mass
one to $\widetilde{M}$, $\gamma=\alpha/\widetilde{M}$ and $\widetilde{M}$ goes
infinity. This case has been extensively investigated in the Bayesian
nonparametric literature from both computational and methodological
perspective (see Ishwaran and Zarepour, 2002, for a thourough discussion).

Furthermore, Eq. (9.1.1) implies that the finite Dirichlet process is a member
of the family of Gibbs partition distributions (Pitman, 1996; Lijoi et al.,
2010). The Gibbs type structure allows us to simplify the prior for number of
occupied components given in Eq. (17), which in the FDMM case becomes

$p^{\star}_{k}=\Pr\\{M^{(a)}=k\\}=V(n,k)\gamma^{k}S_{n,k}^{-1,\gamma}$ (17)

where $S_{n,k}^{-1,\gamma}$ is the _Generalized Stirling number_ computed for
$k$ compositions of $n$ with parameters $-1$ and $\gamma$, i.e.

$S_{n,k}^{-1,\gamma}=\sum_{j=k}^{n}(-1)^{n-j}s_{n,j}S_{j,k}\gamma^{j-k}$

where $s_{n,j}$ is a Stirling number of the first kind and $S_{j,k}$ is a
Stirling number of second kind as defined in Eq. (1.16) and (1.13) of Pitman
(2006) respectively.

Finally, since the finite Dirichlet process is widely used in applications, we
give in Section B of the Appendix a detailed description of the conditional
algorithm when $q_{M}$ is the density of a shifted Poisson distribution and
appropriate hyperpriors are specified on $\gamma$ and $\Lambda$.

####  9.1.2 Uniform weights

Let

$h(s)=\mathbb{I}_{(0,1)}(s)$

i.e. the un-normalized jumps are uniformly distributed. To implement the
conditional algorithm and compute the eppf we need to evaluate the Laplace
transform as well as its derivatives of degree $n$ for each $n>1$. To this
end, we need to solve, for each $n\geq 0$, the following integral

$\displaystyle\int_{0}^{\infty}s^{n}\mathrm{e}^{-us}h(s)ds=\int_{0}^{1}s^{n}\mathrm{e}^{-us}ds=\frac{1}{u^{n+1}}\int_{0}^{u}z^{n}e^{-z}dz=\frac{\gamma(n+1,u)}{u^{n+1}}$
(18)

where $\gamma(\alpha,u)=\int_{0}^{u}z^{\alpha-1}\mathrm{e}^{-z}dz$ is the
upper-incomplete gamma function (Gradshteyn and Ryzhik, 2007). Moreover, for
$n\geq 0$, the upper incomplete gamma function simplifies to

$\gamma(n+1,u)=n!\left[1-\mathrm{e}^{-u}\sum_{m=0}^{n}\frac{u^{m}}{m!}\right]$

Exploiting the above result leads to

$\begin{array}[h]{ccc}\psi(u)=\frac{1-e^{u}}{u},&\text{and}&\kappa(n_{j},u)=\frac{\gamma(n_{j}+1,u)}{u^{n_{j}+1}}\end{array}$

The evaluation of $\psi$ is essential to implement the conditional alogorithm
in the case of uniform weights. Moreover, an efficient implementation of the
conditional algorithm requires us to be able to sample from the tilted Gamma
version of $h$ in Eq. (18), which in this case is simply a truncated Gamma
distribution on $(0,1)$.

Applying Eq. (8), we can also easily obtain the eppf of the process with
uniform jumps

$\pi(n_{1},\dots,n_{k})=\int_{0}^{\infty}\frac{u^{n-1}}{\Gamma(n)}\frac{\Psi(n,k)}{u^{n+k}}\prod_{j=1}^{k}\gamma(n_{j}+1,u)$

where $\Psi(n,k)$ is defined in Section 6.1.

####  9.1.3 Gamma approximation

Any absolutely continuous density on ${\mathbb{R}}^{+}$ can be approximated by
a mixture of Gamma densities. Indeed, DeVore and Lorentz (1993, p.14) showed
that if $h(s)$ defined on $(0,\infty)$ has limit zero as $s\rightarrow\infty$,
then $\mathscr{S}_{\varepsilon}(s)$, defined as

$\mathscr{S}_{\varepsilon}(s)=\mathrm{e}^{-s/\varepsilon}\sum_{l=0}^{\infty}\frac{s^{l}}{\varepsilon^{l}l!}h\left(\varepsilon
l\right)=\sum_{l=0}^{\infty}\varepsilon h(\varepsilon
l)\text{Gamma}(s;l+1,\frac{1}{\varepsilon})$ (19)

admits as limit

$\lim_{\varepsilon\rightarrow 0}\mathscr{S}_{\varepsilon}(s)=h(s)$

uniformly for $0<s<\infty$. Therefore, $h$ can be approximated by a mixture of
Gamma densities, i.e Gamma$(s;l+1,1/\varepsilon)$, with unnormalized weights
$\varepsilon h(\varepsilon l)$, for $l=0,1,\dots$ . See Wiper et al. (2001)
for an extensive discussion of mixtures of Gamma distributions and their
convergence properties. This is a powerful result as we can approximate any
$h$ with a mixture of Gamma densities and allows us to consider a large class
of weight distributions at the cost of computational complexity. In practice,
to approximate $h$ we need to set a tolerance level $\varepsilon$ sufficiently
small. Let $\mathscr{C}=\left(\sum_{l=0}^{\infty}\varepsilon h(\varepsilon
l)\right)^{-1}$. It is easy to show that

$\psi(u)=\mathscr{C}\sum_{l=0}^{\infty}\left(\frac{1}{1+\varepsilon
u}\right)^{l+1}\varepsilon h(\varepsilon l)$

Moreover, it is obvious that $\mathrm{e}^{-us}\mathscr{S}_{\varepsilon}(s)$
and $s^{n_{j}}\mathrm{e}^{-us}\mathscr{S}_{\varepsilon}(s)$ are both infinite
mixtures of Gamma densities. The normalising constant of the first one is
$\psi(u)$, while the normalizing constant of the second one is given by the
function $\kappa$:

$\kappa(n_{j},u)=\mathscr{C}\left(\frac{\varepsilon}{1+u\varepsilon}\right)^{n_{j}}\sum_{l=0}^{\infty}\frac{\Gamma(n_{j}+l+1)}{\Gamma(l+1)}\left(\frac{1}{1+\varepsilon
u}\right)^{l+1}\varepsilon h(\varepsilon l)$

###  9.2 Choice of $\psi$: Point processes with infinite divisible jumps

It is well known (it follows from the Levy-Khinchine formula in Jacod and
Shiryaev, 2013), that the Laplace transform of an infinite divisible random
variable $S$ has the form

$\psi(u)=\exp\left\\{-\int_{0}^{\infty}(\mathrm{e}^{uz}-1)\omega(z)dz\right\\}$

where the _Lévy intensity_ $\omega(z)$ is the intensity of a measure on the
positive reals satisfying the regularity condition
$\int_{0}^{\infty}\min(1,z)\omega(dz)<\infty$. Moreover, if the distribution
of $S$ is absolutely continuous (with respect to the Lebesgue measure) with
strictly positive and continuous density $h(s)$ on ${\mathbb{R}}^{+}$, then
$\int_{0}^{\infty}\omega(z)dz=\infty$. As an alternative to specifying $h$ we
can choose $\omega$ that uniquely identifies $h$, for instance via the
integral equation:

$h(s)=\int_{0}^{s}\omega(z)h(s-z)\frac{z}{s}dz$

The theory of positive infinite divisible random variables has been very
useful for the study of Normalised Random Measures with Independent Increments
(Regazzini et al., 2003), because $S$ can be written as an infinite sum of
positive random variables, i.e. $0<S=\sum_{j=1}^{\infty}S_{j}<\infty$. For the
Norm-IPPF, we define the law of the unnormalised jumps $S_{j}$, which then
defines the distribution of $S$. This is different from the work of Regazzini
et al. (2003) where the distribution of the weights of the a.s. discrete
random measure is derived by first specifying the law of the normalising
constant $S$. In detail, we want to assign the density $h$ of the unnormalized
weights $S_{j}$ in Section 4 such that the Laplace transform has a closed form
and posterior inference is computationally manageable. We point out how, once
we have a closed form of the Laplace transform $\psi(u)$ and of its
derivatives $k(n_{j},u)$, we can easily compute the eppf using Eq. (8) so that
marginal algorithms are straightforward to implement as discussed in Section
6. On the other hand, if we are able to build a sampler to draw from the
tilted density $\mathrm{e}^{-us}h(s)$ and from the Gamma tilted density
$s^{n_{j}}\mathrm{e}^{-us}h(s)$ the implementation of a conditional algorithm
is straightforward. In the following we present three relevant examples for
which computations are feasible.

####  9.2.1 Gamma Process

For the particular choice of

$\displaystyle\omega(z;\gamma)=\gamma z^{-1}e^{-z}$ (20)

where $\gamma>0$, we obtain that the density $h$ of $S_{j}$ coincides with a
Gamma density with parameters $(\gamma,1)$ density. This is the exact same
situation of 9.1.1.

####  9.2.2 $\sigma$-Stable Process

Consider the Levy density

$\omega(z;\sigma)=\frac{\sigma z^{-\sigma-1}}{\Gamma(1-\sigma)}$ (21)

with $0<\sigma<1$. It is straightforward to show that

$\displaystyle-\log(\psi(u))=\int_{0}^{\infty}(1-\mathrm{e}^{uz})\omega(z,\sigma)dz=u^{\sigma}$

and the Laplace transform is

$\displaystyle\psi(u)=\exp(-u^{\sigma})$ (22)

Pollard (1946) shows that the density of $S_{j}$ with Laplace transform in Eq.
(22) can be represented as follows:

$\displaystyle
h(s;\sigma)=\frac{-1}{\pi}\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\sin(\pi\sigma
k)\frac{\Gamma(\alpha k+1)}{s^{\sigma k+1}}$

Although the density $h$ is computationally intractable, since a closed form
expression for the Laplace transform is available, it is possible to implement
a marginal algorithm by calculating the derivatives of $\psi$. Exploiting Eq
(13) in Favaro et al. (2015) we obtain that

$k(n_{j},u)=\frac{\mathrm{e}^{-u^{\sigma}}}{\sigma^{n_{j}}}\sum_{k=1}^{n_{j}}u^{\sigma
k}\mathcal{C}(n_{j},k;\alpha,0)$

where, for any non-negative integer $n\geq 0$, $0\leq k\leq n$, real numbers
$\alpha$, $\beta$, $\mathcal{C}(n,k;\alpha,\beta)$ denotes the noncentral
generalized factorial coefficient (see Charalambides (2005) for details). Here
we mention that these indices can be easily computed when $\beta=0$ using the
recursive formula

$\mathcal{C}(n,k,\alpha,0)=\alpha\mathcal{C}(n-1,k-1;\alpha,0)+(n-1-k\alpha)\mathcal{C}(n-1,k;\alpha,0)$

with $\mathcal{C}(1,1,\alpha,0)=\alpha$.

To implement a conditional algorithm we need to be able to sample from an
Exponential tilted stable density $\mathrm{e}^{-us}h(s)$ (also known as
generalized Gamma) as well as from a Gamma tilted density
$s^{-n_{j}}\mathrm{e}^{-us}h(s)$. Strategies to sample from an Exponential
tilted density are presented in Devroye (2009) and Hofert (2011) while a
method to sample from the Gamma tilted $\sigma$-stable density is discussed in
Section 3.1 of Favaro et al. (2015).

####  9.2.3 Bessel Process

Consider the intensity

$\omega(z;\alpha,\beta)=\frac{\alpha}{z}\mathrm{e}^{-\beta z}I_{0}(z),\quad
z>0$

where $\beta\geq 1$, $\alpha>0$ and let

$I_{\alpha}(z)=\sum_{l=0}^{+\infty}\frac{(z/2)^{2l+\alpha}}{l!\Gamma(\alpha+l+1)}$

be the modified Bessel function of order $\alpha\geq 0$ (see Erdélyi et al.,
1953, Section 7.2.2). Then, for $s>0$,

$\omega(z;\alpha,\beta)=\frac{\alpha}{z}\mathrm{e}^{-\beta
z}+\sum_{l=1}^{+\infty}\frac{\alpha}{2^{2l}(l!)^{2}}z^{2l-1}\mathrm{e}^{-\beta
z}$ (23)

so that $\omega$ is the sum of the Lévy intensity of a Gamma process and of
the Lévy intensities

$\omega_{l}(z;\alpha,\beta)=\frac{\alpha}{2^{2l}(l!)^{2}}z^{2l-1}\mathrm{e}^{-\beta
z},\quad z>0,\qquad l=1,2,\ldots$ (24)

corresponding to finite activity Poisson processes (see Argiento et al.,
2016).

######  Proposition 1.

When $\alpha>0$, we have

  1. (a)

the density $h$ corresponding to the Lévy intensity $\omega(z;\alpha,\beta)$
is

$h(s)=\alpha\left(\beta+\sqrt{\beta^{2}-1}\right)^{\alpha}\dfrac{\mathrm{e}^{-\beta
s}}{s}I_{\alpha}(s),\ \ s>0$

  2. (b)

the Laplace transform of $S_{m}\sim h(s)$ is

$\psi(u)=\left(\frac{\beta+\sqrt{\beta^{2}-1}}{\beta+u+\sqrt{(\beta+u)^{2}-1}}\right)^{\alpha}$

  3. (c)

the function $\kappa(n_{j},u)$ has the following expression

$\kappa(n_{j},u)=\frac{\alpha\left(\beta+\sqrt{\beta^{2}-1}\right)^{\alpha}}{2^{\alpha}(u+\beta)^{n_{j}+\alpha}}\frac{\Gamma(\alpha+n_{j})}{\Gamma(\alpha+1)}\
_{2}F_{1}\left(\frac{n_{j}+\alpha}{2},\frac{n_{j}+\alpha+1}{2};\alpha+1;\frac{1}{(u+\beta)^{2}}\right)$

where

$\
{}_{2}F_{1}(\alpha_{1},\alpha_{2};\gamma;z):=\sum_{l=0}^{\infty}\frac{\left(\alpha_{1}\right)_{l}\left(\alpha_{2}\right)_{l}}{\left(\gamma\right)_{l}}\frac{1}{l!}\left(z\right)^{l},\quad\textrm{
with }(\alpha)_{l}:=\frac{\Gamma(\alpha+l}{\Gamma(\alpha)}$

is the hypergeometric series (see Gradshteyn and Ryzhik, 2007, Equation
(9.100)).

Proof: See Appendix $\blacksquare$

Using Eq. (8) we can derive the eppf for the Bessel case:

$\displaystyle\pi(n_{1},\dots,n_{k})$ $\displaystyle=$
$\displaystyle\int_{0}^{+\infty}\dfrac{u^{n-1}}{\Gamma(n)}\Psi(u,k)\prod_{j=1}^{k}\kappa(n_{j},u)du$
$\displaystyle=$
$\displaystyle\int_{0}^{+\infty}\dfrac{u^{n-1}}{\Gamma(n)(u+\beta)^{n+k\alpha}}\left(\frac{\alpha(\beta+\sqrt{\beta^{2}-1})^{\alpha}}{2^{\alpha}}\right)^{k}\Psi(u,k)$
$\displaystyle\cdot$
$\displaystyle\prod_{j=1}^{k}\frac{\Gamma(\alpha+n_{j})}{\Gamma(\alpha)}\
_{2}F_{1}\left(\frac{n_{j}+\alpha}{2},\frac{n_{j}+\alpha+1}{2};\alpha+1;\frac{1}{(u+\beta)^{2}}\right)du$

Therefore, the difficulty in implementing a marginal algorithm is simply in
the evaluation of the $\ {}_{2}F_{1}$ function.

To implement a conditional algorithm, also in this case we need to sample from
an Exponential tilted density $\mathrm{e}^{-us}h(s)$ and a Gamma tilted
$s^{n_{j}}\mathrm{e}^{-us}h(s)$ density for each real $u>0$ and integer
$n_{j}>1$. From the proof of Proposition 1, it is clear that these densities
are both a mixture of Gamma distributions. In particular, we have that

$\mathrm{e}^{-us}h(s)\propto\sum_{l=1}^{\infty}o_{l}(u)\text{Gamma}(s;2l+\alpha,u+\beta)$

and

$s^{n_{j}}\mathrm{e}^{-us}h(s)\propto\sum_{l=1}^{\infty}o_{l}(u,n_{j})\text{Gamma}(s;2l+n_{j}+\alpha,u+\beta)$

where the mixture weights are given by

$\displaystyle
o_{l}(u)=\frac{\alpha}{2^{\alpha}}\left(\frac{\beta+u+\sqrt{(\beta+u)^{2}-1}}{\beta+u}\right)^{\alpha}\frac{1}{l!\Gamma(\alpha+l+1)2^{2l}}\frac{\Gamma(2l+\alpha)}{(\beta+u)^{2l}}$
$\displaystyle o_{l}(u,n_{j})=\frac{1}{\
{}_{2}F_{1}\left(\frac{n_{j}+\alpha}{2},\frac{n_{j}+\alpha+1}{2};\alpha+1;\frac{1}{(u+\beta)^{2}}\right)}\frac{\Gamma(\alpha+1)}{\Gamma(\alpha+n_{j})}\frac{1}{l!\Gamma(\alpha+l+1)2^{2l}}\frac{\Gamma(2l+\alpha+n_{j})}{\left(u+\beta\right)^{2l}}$

To sample from these two infinite mixtures we can use rejection samplings as
the weights go quickly to zero.

##  10 Galaxy data

We illustrate our model using the Galaxy dataset (Roeder, 1990), which offers
a standard benchmark for mixture models. It contains $n=82$ measurements on
velocities of different galaxies from six well-separated conic sections of
space. Values are expressed in Km/s, scaled by a factor of $10^{-3}$. We fit
Model (15), using a Gaussian density $\mathcal{N}(\mu,\sigma^{2})$ on
$\mathbb{R}$ as $f(y\mid\tau)$, $\tau=(\mu,\sigma^{2})$. We specify the
following prior
$p_{0}(\mu,\sigma^{2})=\mathcal{N}(\mu;m_{0},\sigma^{2}/\kappa_{0})\times\text{Inv-
gamma}(\sigma^{2};\nu_{0}/2,\nu_{0}/2\sigma^{2}_{0})$. Here $\text{Inv-
Gamma}(a,b)$ denotes the Inverse-Gamma distribution with mean $b/(a-1)$ (if
$a>1$). We set $m_{0}=\bar{x}_{n}=20.8315$, $\kappa_{0}=0.01$, $\nu_{0}=4$,
$\sigma_{0}^{2}=4$. Finally we assume a shifted $\text{Poisson}(\Lambda)$ as
prior on $M$ and a Gamma$(\gamma,1)$ as a prior for $S_{m}$ (i.e. a finite
Dirichlet process as mixing distribution). We implement the conditional
algorithm described in Appendix B to perform posterior inference. In
particular, we focus on density estimation and inference on the number of
mixture components and clusters.

First of all, we fit the model with $\Lambda$ and $\gamma$ fixed, with the aim
of comparing the performance of our algorithm with the reversible jump sampler
of Richardson and Green (1997) as implemented in the mixAK R-package (Komárek,
2009). Implementation of our algorithm has been done in C++ using the Rcpp
library (Eddelbuettel and François, 2011), while post processing of the MCMC
results in R. For each MCMC run, we have discarded the first 5000 iterations
as burn-in and thinned every 10, obtaining a final sample size of $5000$. We
have considered different scenarios, and in Figure 2 we show the predictive
density with 95% credible bounds for one of them.

First of all, we fix the hyperparameters $\gamma$ and $\Lambda$ in Eq. (3) in
such a way that the prior mean for the number of clusters is (A)
$\mathbb{E}(k)=1$; (B) $\mathbb{E}(k)=5$; (C) $\mathbb{E}(k)=10$.

Figure 2: Density estimation for $m_{0}=\bar{x}_{n}=20.8315$,
$\kappa_{0}=0.01$, $\nu_{0}=4$, $\sigma_{0}^{2}=0.5$. The hyperparameter
settings of the mixing distribution are specified in simulation scenario in
$D.1$ corresponding to the optimal value of the LPML index.

In order to compare the conditional algorithm with the Reversible Jump, we
compute the integrated autocorrelation time (IAC) and the effective sample
size (ESS) for the number $M$ of components for the all combinations of hyper-
parameters. The IAC (Sokal, 1997) index provides a measure of the efficiency
of the sampling algorithm in terms of accuracy of the estimates (see, e.g.,
Kalli et al., 2011). A small absolute value of the integrated autocorrelation
time (near 0) implies good mixing and hence an efficient method. The Effective
Sample Size (ESS, Kong, 1992) provides an estimate of the number of
independent draws from the posterior distribution of a parameter of interest
and small values indicate high autocorrelation between draws, implying that
the estimate of the posterior distribution of that parameter will be poor.
Posterior results are summarised in Table 1: it is evident that our algorithm
outperforms the reversible jump in terms of both the IAC and ESS.

Table 1:  Running times, posterior mean of $M$ and integrated autocorrelation
times $\hat{\rho}$ for the Gibbs sampler (GS) in Appendix B and the Reversible
Jump (RJ) MCMC implemented in the R-package mixAK.

GS RJ $(\Lambda,\gamma)$ $\mathbb{E}(M|data)$ ESS IAC $\mathbb{E}(M|data)$ ESS
$M$ IAC $M$ $(100,2e^{-4})$ 103.78 5474.59 1.49 70.48 3.07 602.26 A
$(10,2e^{-3})$ 13.51 5000.00 1.51 11.43 16.38 149.08 $(1,10e^{-2})$ 4.22
1854.67 1.03 3.71 298.23 7.78 $(100,1e^{-2})$ 103.67 4231.39 0.58 95.84 8.66
206.10 B $(10,0.143)$ 13.51 1551.56 1.38 10.18 467.97 5.77 $(5,0.5)$ 8.69
1178.10 2.14 7.09 1168.48 1.99 $(1000,2.8e^{-3})$ 1001.01 5393.03 1.50 819.49
1.94 874.78 C $(100,3.2e^{-2})$ 101.62 3846.80 0.57 86.23 22.55 86.35
$(10,1.8)$ 13.71 1595.40 1.64 8.50 1271.37 2.03

Moreover, through an appropriate choice of $(\Lambda,\gamma)$, we are able to
introduce in the model a desired level of sparsity. In Figure 3 we report the
posterior distribution of the number of clusters (allocated components) for
the same combinations of hyper-parameters in Table 1. It is clear that the
posterior distribution of the number of clusters is robust to the choice of
hyper-parameters within each scenario (A, B and C), since the prior mean on
the number of allocated components is constant.

Figure 3: Posterior distribution of $k$ for the three scenarios.

To gain more insight, in Figure 4 we show the posterior distribution of
$M^{(na)}$, the number of non-allocated components. We highlight: (i) these
posteriors are more concentrate on large number for large values of $\Lambda$
(ii) for the same value of $\Lambda$ the level of sparsity increases for small
values of $\gamma$ (see variations within columns). Large values of $\Lambda$
and small values for $\gamma$ favour a posterior distribution for $M^{(na)}$
centred on large values.

Figure 4: Posterior distribution of $M^{(na)}$ for the three scenarios.

We conclude that $\Lambda$ controls the number of unallocated clusters, while
$\gamma$ controls degree of sparsity of the mixture.

Finally, we fit the same model adding an extra layer to the hierarchy, by
specifying prior distributions on both $\gamma$ and $\Lambda$. We consider two
scenarios described in Table 2. In this case we consider the posterior
distribution of the number of allocated components, that is more meaningful
from an inferential point of view, as well as a predictive goodness of fit
criterion, the Logarithm of the Pseudo Marginal Likelihood (LPML – Geisser and
Eddy, 1979).

Table 2:  Results obtained specifying a prior on $\Lambda$ and $\gamma$.

$\Lambda$ $\gamma$ $\mathbb{E}(k|data)$ ESS $k$ IAC $k$ LPML Gamma$(1,0.01)$
Gamma$(1,1)$ 13.41 909.07 3.39 -1632.44 D Gamma$(1,0.1)$ Gamma$(1,1)$ 11.73
723.87 3.65 -1515.55 Gamma$(1,1)$ Gamma$(1,1)$ 7.72 307.77 10.38 -1092.03
Gamma$(1,0.01)$ Gamma$(0.1,0.1)$ 11.79 454.57 5.20 -1471.48 E Gamma$(1,0.1)$
Gamma$(0.1,0.1)$ 11.41 397.18 7.15 -1443.14 Gamma$(1,1)$ Gamma$(0.1,0.1)$ 7.89
199.44 12.79 -1139.50

In Figure C.8 of Appendix C we show the posterior distribution of number of
clusters: we note how adding an extra layer to the hierarchy makes inference
more robust to hyper-prior specifications. In Figure C.9 of Appendix C the
posterior distribution of non-allocated components is shown: it is evident
that $\Lambda$ still influences such distribution, while $\gamma$ determines
the level of sparsity. Moreover, adding this extra level of randomness induces
more parsimonious posteriors: the posterior on the number of non-allocated
components is now shrunk toward zero in both scenarios D and E. Although
treating $\Lambda$ and $\gamma$ as random variables leads to more robust
estimates, it also increases the autocorrelation in the MCMC chains. This is
evident from Table 2 as well as from Figure 5 that shows the joint marginal
posterior of $\Lambda$ and $\gamma$ for scenarios D and E. We highlight the
strong negative correlation between the two hyper-parameters, which is natural
as $\Lambda$ controls the number of on non-allocated components while $\gamma$
the number of allocated ones.

Figure 5: Three-dimensional histogram of the joint marginal posterior of
$\Lambda$ and $\gamma$ for scenarios D and E.

##  11 Population Structure: Taita Thrush Data

In population genetics, population structure refers to the presence of a
systematic differences in genetic markers’ allele frequencies between
subpopulations due to variation in ancestry. This phenomenon arises from the
bio-geographical distribution of species, due to the fact that either natural
populations occupy a vast geographic area and cannot act as randomly mating or
geographical barriers reduce migration between different regions. Consequently
population structure affects the dynamics of alleles in populations and
impacts the type of statical analysis to perform in many applications, for
example in genetic association studies. Broadly speaking, the analysis of
population structure focuses on: (i) detecting population structure in a
sample of chromosomes; (ii) estimating the number of populations in a sample;
(iii) assigning individuals to populations and (iv) defining the number of
ancestral populations in a sample. A variety of statistical approaches have
been proposed to infer population structure. Arguably the most widely method
is the one proposed by Pritchard et al. (2000) based on Bayesian mixture
models and implemented in the software STRUCTURE (Pritchard and Wen, 2003).
Pritchard et al. (2000) assume that individuals come from one of $M$ (fixed)
subpopulations and population membership and population specific allele
frequencies are jointly estimated from the data. Independent priors on the
allelic profile parameters of each population are specified and posterior
inference is performed through MCMC. In Pritchard et al. (2000), the number of
mixture components is fixed and their method clusters individual in one of a
fixed number of populations. Determination of the number of populations in a
sample is achieved using a model selection criteria based on MCMC estimates of
the log marginal probabilities of the data and the Bayesian deviance
information criterion, though it has been noted by Falush et al. (2003) that
such estimates are highly sensitive to prior specifications regarding the
relatedness of the populations. To avoid such model selection, Huelsenbeck and
Andolfatto (2007) propose a method for the analysis of population structure
based on a Dirchlet process mixture model and implemented in the software
Structurama (Huelsenbeck et al., 2011), which does not require the
specification of a fixed and finite $M$.

We now illustrate the performance of our method in a population structure
problem, using an empirical data set of $n=237$ Taita thrushes kindly made
available by Dr P. Galbusera. A previous smaller version of these data
(Galbusera et al., 2000) has been analysed by Pritchard et al. (2000) and
Huelsenbeck and Andolfatto (2007) as benchmark example. We have run an
analysis using our algorithm on this old data and have drawn identical
conclusions. Here we prefer to focus on the new dataset. The Taita Hills in
Kenya represent the northernmost part of the Eastern Arc Mountains
biodiversity hotspot of Kenya and Tanzania. They are isolated from other
highlands by over 80 km of semiarid plains in either direction. During the
last 200 years, indigenous forest cover in the Taita Hills has decreased by
circa 98% and the critically endangered Taita thrush, endemic to the Taita
Hills, is currently restricted to the fragments of Mbololo, Ngangao and Chawia
(Callens et al., 2011). These fragments are separated from each other by
cultivated areas and human settlements. This dataset is ideal to test the
performance of our method as the geographic samples are likely to represent
distinct populations, i.e. mixture components. Each bird was sampled at $L=6$
microsatellite loci. We follow the notation of Huelsenbeck and Andolfatto
(2007). Recall that the Taita thrush is diploid, i.e. has two sets of
chromosomes and for each locus we have genotype data. At locus $l$, we observe
$J_{l}$ unique alleles. The number of copies of allele $j$ at locus $l$ in
individual $i$ is denoted by $Y_{ilj}\in\\{0,1,2\\}$ and the number of copies
of all alleles observed at locus $l$ in individual $i$ is denoted by
$Y_{il}=\sum_{j=1}^{J_{l}}Y_{ilj}$. The allelic information for individual $i$
at locus $l$ is contained in the vector
${Y}_{il}=(Y_{il1},Y_{il2},\dots,Y_{ilJ_{l}})$, with the constrain
$\sum_{j=1}^{J_{l}}Y_{ilj}=2$. Given $M$ possible populations, let
$\tau_{mlj}$ denote the frequency of allele $j$ al locus $l$ in population
$m$, let $\tau_{ml}=(\tau_{ml1},\ldots,\tau_{mJ_{l}})$ be the vector of allele
frequencies at locus $l$ in population $m$ and let
$\tau_{m}=(\tau_{m1},\ldots,\tau_{mL})$. Finally, let
$c_{i}\in\\{1,\ldots,M\\}$ be the allocation variable of bird $i$, i.e.
$c_{i}=m$ if the bird comes from population $m$. Following Huelsenbeck and
Andolfatto (2007) we assume that

$f\left({y}_{il}|{\tau}_{ml}\right)=\mathbb{P}({Y}_{il}={y}_{il}\mid{\tau}_{ml},c_{i}=m)\propto\prod_{j=1}^{J_{l}}\tau_{mlj}^{y_{ilj}}$

We assume independence across loci, so that, if
${Y}_{i}=({Y}_{i1},\dots,{Y}_{iL})$ is the multidimensional array of the
allelic information at the $L$ loci for individual $i$, we have

$f\left({y}_{i}\mid{\tau}_{m},c_{i}=m\right)=\mathbb{P}({Y}_{i}={y}_{i}\mid{\tau}_{m},c_{i}=m)=\prod_{l=1}^{L}f\left({y}_{il}\mid{\tau}_{ml}\right)$
(25)

We fit Model (3), with the sampling model defined in Eq. (25). The mixing
measure is a finite Dirchlet process as in Section 9.1.1, with the following
prior specification: $M$ has a shifted Poisson prior distribution with
parameter $\Lambda$, $P_{0}$ is the convolution of $L$ independent Dirichlet
distributions with parameter 1, $\gamma$ in the finite Dirichlet process has a
Gamma prior with parameter $(0.1,0.1)$, $\Lambda$ has a Gamma prior with
parameter $(3/2,1/2)$. For the parameter $\gamma$ we have specified a vague
prior distribution, while the hyper-parameters in the prior for $\Lambda$ are
chosen so that the prior mean is 3, corresponding to the three geographical
fragments, and the prior variance is large. We employ the conditional
algorithm described in Appendix B to perform posterior inference. The mode of
the posterior distribution for $k$ is at 3 ($\mathbb{E}(k\mid
data)=3,\operatorname{Var}(k\mid data)=0.03$), as well as the one of the
posterior of $M$ ($\mathbb{E}(M\mid data)=3.12,\operatorname{Var}(M\mid
data)=0.42$). From Figure 6 it is evident that the three clusters coincide
with the three geographical fragments, except in three cases where the birds
appear to be out of the obvious clusters. This could be due to rare migration
events (Galbusera et al., 2000).

Figure 6: Posterior estimate of the clustering allocation: each colour
correspond to a cluster. Note that he two green thrushes have been captured in
Ngangao, but have the genetic profile of the Mbololo birds. The opposite is
true for thrush B1441.

An important goal of population structure analysis is not only to uncover the
group structure of the observations, but also to identify variables that best
distinguish the different populations. The results could lead to a better
understanding of the evolutionary patterns of population differentiation. To
this end we would like to identify the microsatellite loci that most influence
the clustering structure. Variable selection for clustering is a challenging
problem since there is no observed response to inform the selection and the
inclusion of unnecessary variables could complicate or mask the recovery of
the clusters (Tadesse et al., 2005; Kim et al., 2006). As such there are few
contributions in the literature. Here we opt for a model choice method
proposed by Goutis and Robert (1998) in the generalised linear model
framework, which we adapt to our context. The approach of Goutis and Robert
(1998) focuses on the predictive properties of a model and, employing the
Kullback-Leibler distance as discrepancy measure, aims to assess the relevance
of some restriction on the parameter $\Theta$ (leading to a simpler model)
with respect to a full model described by a density $f(y\mid\theta)$. More in
details, for each locus $l$, let $Y_{l}=(Y_{1l},\ldots,Y_{nl})$,
$\theta_{il}=\tau_{ml}$ if $c_{i}=m$ and
$\theta_{l}=(\theta_{1l},\ldots,\theta_{nl})$. Let $f(y_{l}\mid\theta_{l})$ be
the full general mixture model:

$f(y_{l}\mid\theta_{l})\propto\prod_{i=1}^{n}\prod_{j=1}^{J_{l}}\theta_{ilj}^{y_{ijl}}$
(26)

We define a model choice hypothesis $H_{0}$ through a restriction on the
parameter space, i.e. $\theta_{l}\in\Theta_{0}\subset\Theta$, where
$\Theta_{0}$ is the subset of the parameter space such that
$\theta_{ilj}=\widetilde{\theta}_{lj}$ for each $i$. In our application
$H_{0}$ represents a fuly parametric model for locus $l$. Goutis and Robert
(1998) define the projection $\theta_{l}^{\perp}$ of $\theta_{l}$ according to
the Kullback-Leibler distance $d$ to be the point in $\Theta_{0}$ that
achieves the infimum

$d\left\\{f(\cdot\mid\theta_{l}),f(\cdot\mid\theta_{l}^{\perp})\right\\}=\inf_{\widetilde{\theta}_{l}\in\Theta_{0}}d\left\\{f(\cdot\mid\theta_{l}),f(\cdot\mid\widetilde{\theta}_{l})\right\\}$

where
$\widetilde{\theta}_{l}=(\widetilde{\theta}_{l1},\ldots,\widetilde{\theta}_{lJ_{l}})$
and $f(\cdot\mid\theta_{l}^{\perp})$ is the projection of
$f(\cdot\mid\theta_{l})$. Obviously small values of $d$ support $H_{0}$. We
opt for this approach because, instead of phrasing the problem in terms of the
classical dichotomy between null and alternative hypothesis, it interprets
model choice in terms of the approximation efficacy of a more parsimonious
model, focusing on whether or not $\theta_{l}$ is far away from the subspace
$\Theta_{0}$. In Figure 7

Figure 7: Posterior distribution of the KL divergence for each microsatellite
locus.

we show the posterior distribution of
$d\left\\{f(\cdot\mid\theta_{l}),f(\cdot\mid\theta_{l}^{\perp})\right\\}$ for
each locus $l$. It is evident that locus PC3 contributes the least to the
clustering structure as the distance is concentrated near zero, implying the
its allele frequencies are similar across Taita thrush populations. The other
loci, in particular PAT43, present allele frequency differences among the
three groups, which in our case well correspond to geographical locations.

##  12 Conclusions

In this work we contribute to the growing understanding of mixture models by
providing an unifying framework which encompasses both finite and infinite
mixtures. A key concept is the distinction between the number of components
and number of clusters, where by components we refer to the number of
subpopulations that are likely to have generated the data, while clusters
indicate the number of non-empty components in a sample. Already Nobile et al.
(2004) had pointed out this difference, noticing that the posterior
distribution of the number of components $M$ (corresponding to the data
generating process) might assign considerable probability to values greater
than $k$, the number of clusters. More recently, the concept of sparse finite
mixtures has been introduced as a first attempt to bridge between finite
mixture models and nonparametric mixtures (see Malsiner-Walli et al., 2016,
2017). In this context, Frühwirth-Schnatter and Malsiner-Walli (2018) fix $M$
very large so that they are close to infinity, justifying this choice from the
asymptotic point of view and then their work focuses on sparse estimation of
the number of clusters. Our construction is based on the normalization of a
point process, which is a standard trick in Bayesian nonparametrics. We
introduce the Norm-IFFP prior process and we provide theoretical results
characterizing the induced prior on the partition of the observations and the
posterior distribution of this process. Our framework allows for efficient
computations (inherited from the nonparametric construction) and for data
driven estimation of both number of clusters and components, as well as of any
functional of interest.

## Acknowledgement

We would like to thank Dr Peter Galbusera at the Royal Zoological Society of
Antwerp for sharing the enriched Taita Thrush Dataset. Dr Argiento is grateful
to Yale-NUS College, Singapore for the funding provided.

##  Appendix A Appendix: Proofs

###  A.1 Proof of Theorem 1

###### Proof.

We have

$\displaystyle\pi(n_{1},\dots,n_{k})=\sum_{m=1}^{+\infty}\pi(n_{1},\dots,n_{k}|M=m)q_{m}$
(A.27)

since $M\sim q_{M}$ and we have assumed $q_{0}=0$. Then, equation (3) in
Pitman (2003) yields

$\pi(n_{1},\dots,n_{k}|M=m)=\mathbb{I}_{\\{k\leq
m\\}}\sum_{(c_{1}^{\star},\dots,c_{k}^{\star})}\mathbb{E}\left(\prod_{j=1}^{k}w_{c_{j}^{\star}}^{n_{j}}\right)$

where the vector $(c_{1}^{\star},\dots,c_{k}^{\star})$ ranges over all
permutations of $k$ positive integers in $\\{1,\dots,m\\}$. Recall that
$w_{m}=\frac{S_{m}}{T}$ as defined in Eq.(7). Under the assumptions in Section
3 and Section 4 the joint law of the unnormalized jumps of $P$ is
$\mathcal{L}(dS_{0},\dots,dS_{m}|M=m)=\prod_{l=1}^{m}h(S_{l})dS_{l}$. Then,
using the identity $1/T^{n}=\int_{0}^{+\infty}1/\Gamma(n)u^{n-1}e^{-uT}du$, we
have:

$\displaystyle\pi(n_{1},..,n_{k}|M=m)$ $\displaystyle=\mathbb{I}_{\\{k\leq
m\\}}\sum_{(c_{1}^{\star},\dots,c_{k}^{\star})}\int\prod_{j=1}^{k}\frac{S_{c_{j}^{\star}}^{n_{j}}}{T^{n_{j}}}\mathcal{L}(dS_{1},\dots,dS_{m}|m)$
$\displaystyle=\mathbb{I}_{\\{k\leq
m\\}}\sum_{(c_{1}^{\star},\dots,c_{k}^{\star})}\int_{0}^{+\infty}\left(\frac{1}{\Gamma(n)}u^{n-1}\prod_{j=1}^{k}\int_{0}^{+\infty}S_{c_{j}^{\star}}^{n_{j}}\mathrm{e}^{-S_{c_{j}^{\star}}u}h(S_{c_{j}^{\star}})dS_{c_{j}^{\star}}\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\left.\prod_{c\notin(c_{1}^{*},\dots,c_{k}^{\star})}\int_{0}^{+\infty}\mathrm{e}^{-S_{c}u}h(S_{c})dS_{c}\right)du$
$\displaystyle=\mathbb{I}_{\\{k\leq
m\\}}\sum_{(c_{1}^{\star},\dots,c_{k}^{\star})}\int_{0}^{+\infty}\frac{1}{\Gamma(n)}u^{n-1}\prod_{j=1}^{k}\mathbb{E}\left(S_{c_{i}^{\star}}^{n_{j}}\mathrm{e}^{-S_{c_{j}^{\star}}u}\right)\prod_{c\notin(c_{1}^{*},\dots,c_{k}^{\star})}\mathbb{E}(\mathrm{e}^{-S_{c}u})du$
$\displaystyle=\mathbb{I}_{\\{k\leq
m\\}}\sum_{(c_{1}^{\star},\dots,c_{k}^{\star})}\int_{0}^{+\infty}\frac{1}{\Gamma(n)}u^{n-1}\prod_{j=1}^{k}\kappa(n_{j},u)\psi(u)^{m-k}du$
$\displaystyle=\mathbb{I}_{\\{k\leq
m\\}}\frac{m!}{(m-k)!}\int_{0}^{+\infty}\frac{1}{\Gamma(n)}u^{n-1}\prod_{j=1}^{k}\kappa(n_{j},u)\psi(u)^{m-k}du$

where $\psi(u)=\int_{0}^{\infty}e^{-us}h(s)ds$ is the Laplace transform of the
density $h(s)$, while for each integer $n$,
$\kappa(n_{,}u):=\int_{0}^{\infty}s^{n}\mathrm{e}^{-us}h(s)ds=(-1)^{n}\frac{d}{du^{n}}\psi(u)$.
We exploit the fact that in the second last equation, the summation term does
not depend on the indices $(c_{1}^{\star},\dots,c_{k}^{\star})$, with
$\sum_{c_{1}^{\star},\dots,c_{k}^{\star}}1=\frac{m!}{(m-k)!}$ since
$(c_{1}^{\star},\dots,c_{k}^{\star})$ ranges over all permutations of $k$
positive integers between $1$ and $m$.

Now, combining Eq. (A.27) with the equality derived above, we obtain

$\displaystyle\pi(n_{1},\dots,n_{k})$
$\displaystyle=\sum_{m=1}^{+\infty}\pi(n_{1},\dots,n_{k}|M=m)q_{m}$
$\displaystyle=\sum_{m=1}^{+\infty}\mathbb{I}_{\\{k\leq
m\\}}\frac{m!}{(m-k)!}\int_{0}^{+\infty}\frac{1}{\Gamma(n)}u^{n-1}\prod_{j=1}^{k}\kappa(n_{j},u)\psi(u)^{m-k}du$
$\displaystyle=\sum_{m=k}^{+\infty}\frac{m!}{(m-k)!}\int_{0}^{+\infty}\frac{1}{\Gamma(n)}u^{n-1}\prod_{j=1}^{k}\kappa(n_{j},u)\psi(u)^{m-k}du$
$\displaystyle=\sum_{m=0}^{+\infty}\frac{(m-k)!}{(m)!}\int_{0}^{+\infty}\frac{1}{\Gamma(n)}u^{n-1}\prod_{j=1}^{k}\kappa(n_{j},u)\psi(u)^{m}du$

which gives Eq. (8) and concludes the proof.

∎

###  A.2 Proof of Corollary 1

###### Proof.

Let $(n_{1},\dots,n_{k})$ be a composition of $n$, i.e. a vector of $k$
positive integers $(n_{j}>0\ i=1,\dots,k)$ such that $n=\sum_{j=1}^{k}n_{i}$.
There are

${{n}\choose{n_{1},\dots,n_{k}}}\frac{1}{k!}$

partitions of $\\{1,\dots,n\\}$ such that their cluster size is given by with
$(n_{1},\dots,n_{k})$. Thus it is clear that the joint prior probability of
sampling a partition (clustering configuration) with $k$ clusters and cluster
sizes $(n_{1},\dots,n_{k})$ is

${{n}\choose{n_{1},\dots,n_{k}}}\frac{1}{k!}\pi(n_{1},\dots,n_{k})$

The marginal prior probability of sampling a partition with $k$ cluster is
computed summing over all the possible cluster sizes $(n_{1},\dots,n_{k})$:

$\displaystyle p^{\star}_{k}$
$\displaystyle=\sum_{n_{1},\dots,n_{k}}{{n}\choose{n_{1},\dots,n_{k}}}\frac{1}{k!}\pi(n_{1},\dots,n_{k})$
$\displaystyle=\int_{0}^{+\infty}\dfrac{u^{n-1}}{\Gamma(n)}\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}\frac{n!}{k!}\sum_{n_{1},\dots,n_{k}}\prod_{j=1}^{k}\frac{\kappa(n_{j},u)du}{n_{j}!}$
$\displaystyle=\int_{0}^{+\infty}\dfrac{u^{n-1}}{\Gamma(n)}\left\\{\sum_{m=0}^{\infty}\frac{(m+k)!}{m!}\psi(u)^{m}q_{m+k}\right\\}B_{n,k}(\kappa(\cdot,u))$

where $B_{n,k}(\kappa(\cdot,u))$ is the _partial Bell polynomial_ (Pitman,
2006) for the sequence of coefficient $\\{\kappa(n,u),\ n=1,2,\dots\\}$.

∎

###  A.3 Proof of Theorem 2

###### Proof.

Let $\bm{\theta}=(\theta_{1},\dots,\theta_{n})$ be a sample from
$P\sim\text{Norm-IPPF}(h,p_{0},q_{M})$:

$\displaystyle\theta_{1},\dots,\theta_{n}\mid P$
$\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}$ $\displaystyle P$
$\displaystyle P$ $\displaystyle\sim$ $\displaystyle\text{Norm-
IPPF}(h,p_{0},q_{M})$

With a slight abuse of notation, we use $\mathcal{L}(X)$ to denote the pdf or
pmf of a random variable $X$ and $\mathcal{L}(\tilde{P})$ to indicate the
Janossy density of a point process. We need to show that the posterior density
$\mathcal{L}(\tilde{P}|\theta_{1},\dots,\theta_{n})$ is still the Janossy
density of a finite point process. Indeed,

$\begin{split}\mathcal{L}(\tilde{P}|\bm{\theta})&\propto\mathcal{L}(\bm{\theta},\tilde{P})=\mathcal{L}(\theta_{1},\dots,\theta_{n}\mid\tilde{P})\mathcal{L}(\tilde{P})=\prod_{i=1}^{n}P(\theta_{i})\mathcal{L}(\tilde{P})\\\
&=\left\\{\frac{1}{T^{n}}\prod_{i=1}^{n}\left(\sum_{m\in\cal
M}S_{m}\delta_{\tau_{m}}(\theta_{i})\right)\right\\}M!p_{M}\prod_{m\in\cal
M}h(S_{m})p_{0}(\tau_{m})\\\ \end{split}$ (A.28)

We introduce the variable $U_{n}\in{\mathbb{R}}^{+}$ such that the
$\mathcal{L}(U_{n}\mid T)$ is a Gamma$(n,T)$ and whose marginal distribution
is given by Eq. (11). Recall that

$\displaystyle T$ $\displaystyle=$ $\displaystyle\sum_{m\in\cal M}S_{m}$
$\displaystyle\frac{1}{T^{n}}$ $\displaystyle=$
$\displaystyle\int_{0}^{n}\frac{1}{\Gamma(n)}u^{n-1}e^{-uT}\text{d}u$

Then

$\begin{split}\mathcal{L}(\tilde{P},u|\bm{\theta})&\propto\frac{1}{\Gamma(n)}u^{n-1}\mathrm{e}^{-uT}\left\\{\prod_{i=1}^{n}\left(\sum_{m\in\cal
M}S_{m}\delta_{\tau_{m}}(\theta_{i})\right)\right\\}M!q_{M}\prod_{m\in\cal
M}h(S_{m})p_{0}(\tau_{m})\\\
&=\frac{1}{\Gamma(n)}u^{n-1}\mathrm{e}^{-u\sum_{m\in\cal
M}S_{m}}\left\\{\prod_{i=1}^{n}\left(\sum_{m\in\cal
M}S_{m}\delta_{\tau_{m}}(\theta_{i})\right)\right\\}M!q_{M}\prod_{m\in\cal
M}h(S_{m})p_{0}(\tau_{m})\\\
&=\frac{1}{\Gamma(n)}u^{n-1}\left\\{\prod_{i=1}^{n}\left(\sum_{m\in\cal
M}S_{m}\delta_{\tau_{m}}(\theta_{i})\right)\right\\}M!p_{M}\prod_{m\in\cal
M}\mathrm{e}^{-uS_{m}}h(S_{m})p_{0}(\tau_{m})\\\ \end{split}$ (A.29)

Since $(\theta_{1},\dots,\theta_{n})$ is a sample from a discrete
distribution, there is a positive probability of ties among the $\theta_{i}$s.
We denote with
$\bm{\theta}^{\star}=(\theta^{\star}_{1},\dots,\theta^{\star}_{k})$ the vector
of distinct values, with $k\leq\min\\{M,n\\}$. Moreover we denote with
$(n_{1},\dots,n_{k})$ of induced clusters in the sample
$(\theta_{1},\dots,\theta_{n})$, i.e.
$n_{j}=\\#\\{i:\theta_{i}=\theta^{\star}_{j}\\}$, $j=1,\dots,k$.
Let ${\cal M}^{(a)}$ be the set of indexes of the allocated points,
$k=M^{(a)}\leq M$. Moreover, conditionally on $M$, the process
$\widetilde{P}=\\{(S_{1},\tau_{1}),\ldots,(S_{M},\tau_{M})\\}$ is defined on
an $M$-dimensional space. Note that $M=M^{(a)}+M^{(na)}$. Without loss of
generality and for ease of notation, we assume that the set ${\cal M}^{(a)}$
is ordered, i.e. ${\cal M}^{(a)}=\\{1,\dots,k=M^{(a)}\\}$. This assumption is
inconsequential as ${\cal M}^{(a)}$ is a set of indices and the names of the
labels are irrelevant. If $m\in{\cal M}^{(a)}$, then
$\tau_{m}=\theta^{\star}_{m}$ and there is a one to one correspondence between
the set $\\{\theta_{m}^{\star},m=1,\ldots,k\\}$ and
$\\{\tau_{m},m=1,\ldots,M^{(a)}\\}$. Let ${\cal M}^{(na)}$ be the set of
indexes corresponding to unallocated jumps. A posteriori, conditionally on
$(\theta_{1},\ldots,\theta_{n})$, ${\cal M}^{(na)}=\\{k+1,\ldots,M\\}$ and
$\widetilde{P}=\widetilde{P}^{(a)}\cup\widetilde{P}^{(na)}$. If $M=k$, then
${\cal M}^{(na)}$ coincides with the empty set. Conditionally on $u$, we can
now obtain the Janossy measure of the posterior distribution:

$\displaystyle\mathcal{L}(\widetilde{P}|u,\bm{\theta})$ $\displaystyle\propto$
$\displaystyle\mathcal{L}(\widetilde{P},u\mid\bm{\theta})$ (A.30)
$\displaystyle\propto$
$\displaystyle\left\\{\prod_{m\in\mathcal{M}^{(a)}}S_{m}^{n_{m}}\delta_{\tau_{m}}(\theta^{\star}_{m})\right\\}(M^{(na)}+k)!q_{M^{(na)}+k}\prod_{m\in\mathcal{M}}\mathrm{e}^{-uS_{m}}h(S_{m})p_{0}(\tau_{m})$
$\displaystyle=$
$\displaystyle\left\\{\prod_{m\in\mathcal{M}^{(a)}}\mathrm{e}^{-uS_{m}}S_{m}^{n_{m}}\delta_{\tau_{m}}(\theta^{\star}_{m})h(S_{m})p_{0}(\tau_{m})\right\\}(M^{(na)}+k)!q_{M^{(na)}+k}$
$\displaystyle\cdot\prod_{m\in\mathcal{M}^{(na)}}\mathrm{e}^{-uS_{m}}h(S_{m})p_{0}(\tau_{m})$
$\displaystyle=$
$\displaystyle{{M^{(na)}+k}\choose{M^{(na)},k}}\left\\{k!\prod_{m\in\mathcal{M}^{(a)}}\mathrm{e}^{-uS_{m}}S_{j}^{n_{m}}\delta_{\tau_{m}}(\theta^{\star}_{m})h(S_{m})p_{0}(\tau_{m})\right\\}$
$\displaystyle\times\left\\{M^{(na)}!q_{M^{(na)}+k}\prod_{m\in\mathcal{M}^{(na)}}\mathrm{e}^{-uS_{m}}h(S_{m})p_{0}(\tau_{m})\right\\}$

This implies that, conditionally on the auxiliary variable $U_{n}=u$ and on
the sample $(\theta_{1},\dots,\theta_{n})$, $\tilde{P}$ is the superpostition
of two spatial point processes. In fact, the Jannossy density of $\tilde{P}$
can be factorized as

$\mathcal{L}(\tilde{P}|u,\bm{\theta})={{M^{(na)}+k}\choose{M^{(na)},k}}\mathcal{L}(\tilde{P}_{a}|u,\bm{\theta})\times\mathcal{L}(\tilde{P}^{(na)}\mid
u,\bm{\theta})$ (A.31)

We now determine the density of the two processes, the one corresponding to
the allocated jumps and the process defining the unallocated jumps. In the
first case, we have that:

$\mathcal{L}(\widetilde{P}_{a}\mid\mathcal{M}^{(a)},u,\bm{\theta})\propto
k!\prod_{m\in\mathcal{M}^{(a)}}\mathrm{e}^{-uS_{m}}S_{j}^{n_{m}}\delta_{\theta^{\star}_{m}}(\tau_{m})h(S_{m}).$

This implies that $\widetilde{P}^{(a)}$ is a set of $k$ independent points. In
particular, we have
$\widetilde{P}^{(a)}=\\{(S_{1},\theta_{1}^{\star}),\dots,(S_{k},\theta_{k}^{\star})\\}$,
for fixed $\theta^{\star}_{m}$s. The $S_{m}$s have density given by

$h_{n_{m}}(s;u)=\left((-1)^{n_{m}}\frac{d}{du^{n_{m}}}\psi(u)\right)^{-1}\mathrm{e}^{-us}s^{n_{m}}h(s),\quad
s>0,$

where

$\psi(u):=\int_{0}^{\infty}e^{-us}h(s)ds,\quad u>0$

is the Laplace transform of the density $h(s)$.

We now focus on the process of unallocated jumps. It is straightforward to
derive the Janossy density of the process $\widetilde{P}^{(na)}$:

$\begin{split}\mathcal{L}(\widetilde{P}^{(na)}\mid\mathcal{M}^{(na)},u,\bm{\theta})&\propto(M^{(na)}+k)!q_{M^{(na)}+k}\prod_{m\in\mathcal{M}_{na}}\mathrm{e}^{-uS_{m}}h(S_{m})p_{0}(\tau_{m})\\\
&=\left\\{\frac{(M^{(na)}+k)!}{M^{(na)}!}\psi(u)^{M^{(na)}}q_{M^{(na)}+k}\right\\}M^{(na)}!\prod_{m\in\mathcal{M}^{(na)}}\psi(u)^{-1}\mathrm{e}^{-uS_{m}}h(S_{m})p_{0}(\tau_{m})\\\
&\propto\widetilde{q}_{M^{(na)}}M^{(na)}!\prod_{m\in\mathcal{M}^{(na)}}\psi(u)^{-1}\mathrm{e}^{-uS_{m}}h(S_{m})p_{0}(\tau_{m})\end{split}$
(A.32)

where

$\widetilde{q}_{M^{(na)}}=\frac{1}{C_{k}}\frac{(M^{(na)}+k)!}{M^{(na)}!}\psi(u)^{M^{(na)}}q_{M^{(na)}+k}$

Since

$C_{k}:=\sum_{M^{(na)}=0}^{\infty}\frac{(M^{(na)}+k)!}{M^{(na)}!}\psi(u)^{M^{(na)}}q_{M^{(na)}+k}<\sum_{M^{(na)}=0}^{\infty}\frac{(M^{(na)}+k)!}{M^{(na)}!}\psi(u)^{M^{(na)}}=\frac{k!}{(1-\psi(u))^{k+1}}$

we can conclude that $\\{\widetilde{q}_{M^{(na)}}\\}$ is a discrete
probability on $\\{0,1,\dots\\}$ and Eq. (A.32) defines a proper Janossy
density. A posteriori, conditionally on $U_{n}$, $M^{(na)}$ is distributed as
$\widetilde{q}_{M^{(na)}}$ and, given $M^{(na)}$, the jumps $S_{m}$s are
i.i.d. from the exponentially tilted distribution defined by $h(s)e^{-uS}$.

This concludes the proof of item (a) and (b). Proof of item (c) follows
directly from Eq. (A.31). Derivation of the posterior distribution of $U_{n}$
follows from the expression of the eppf obtained in Theorem 1. ∎

###  A.4 Proof of Eq. (11)

$\displaystyle
f_{U_{n}}(u,n)=\frac{u^{n-1}}{\Gamma(n)}\mathbb{E}(T^{n}\mathrm{e}^{-Tu})$
$\displaystyle=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\frac{d}{du^{n}}\mathbb{E}\left(\mathrm{e}^{-Tu}\right)$
$\displaystyle=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\frac{d}{du^{n}}\mathbb{E}\left(\mathbb{E}(\mathrm{e}^{-Tu}\mid
M)\right)$
$\displaystyle=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\frac{d}{du^{n}}\mathbb{E}\left(\prod_{m=1}^{M}\mathbb{E}(\mathrm{e}^{-S_{m}u})\right)$
$\displaystyle=\frac{u^{n-1}}{\Gamma(n)}(-1)^{n}\frac{d}{du^{n}}\mathbb{E}\left(\psi(u)^{M}\right)$

###  A.5 Proof of Proposition 1

###### Proof.

The Laplace transform of $S_{m}$ can be rewritten via the Lévy-Khintchine
formula:

$\displaystyle\psi(u)$ $\displaystyle:=\mathbb{E}(\mathrm{e}^{-\lambda
T})=\exp\left\\{-\alpha\int_{0}^{+\infty}(1-\mathrm{e}^{-\lambda
s})\rho(s;\beta)ds\right\\}$
$\displaystyle=\exp\left\\{-\alpha\left(\log\left(\frac{\beta+\lambda}{\beta}\right)+\sum_{m=1}^{+\infty}\frac{\Gamma(2m)}{2^{2m}(m!)^{2}\beta^{m}}-\sum_{m=1}^{+\infty}\frac{\Gamma(2m)}{2^{2m}(m!)^{2}(\beta+\lambda)^{m}}\right)\right\\}$
$\displaystyle=\exp\left\\{-\alpha\log\left(\frac{\beta+\lambda+\sqrt{(\beta+\lambda)^{2}-1}}{\beta+\sqrt{\beta^{2}-1}}\right)\right\\}.$
$\displaystyle=\left(\frac{\beta+\sqrt{\beta^{2}-1}}{\beta+u+\sqrt{(\beta+u)^{2}-1}}\right)^{\alpha}$

The same expression is obtained when $S_{m}\sim
h(s)=\alpha(\beta+\sqrt{\beta^{2}-1})^{\alpha}\dfrac{\mathrm{e}^{-\beta
s}}{s}I_{\alpha}(s)$, $s>0$ (see Gradshteyn and Ryzhik, 2007, Equation
(17.13.112)) and this proves both item (a) and (b). Note that, when $\beta=1$,
$f_{T}$ is called Bessel function density (Feller, 1971).

To prove (c), we use the definition of the function $\kappa$ and obtain that
for each $n_{j}\geq 1$

$\displaystyle\kappa(n_{j},u)$
$\displaystyle=\mathbb{E}(S_{j}^{n_{j}}\mathrm{e}^{-uS_{j}})=\int_{0}^{\infty}s^{n_{j}}\mathrm{e}^{-us}h(s)ds$
$\displaystyle=\alpha(\beta+\sqrt{\beta^{2}-1})^{\alpha}\int_{0}^{\infty}s^{n_{j}-1}\mathrm{e}^{-(u+\beta)s}I_{\alpha}(s)ds$
$\displaystyle=\alpha(\beta+\sqrt{\beta^{2}-1})^{\alpha}\sum_{m=0}^{\infty}\frac{1}{m!\Gamma(\alpha+m+1)2^{2m+\alpha}}\int_{0}^{\infty}s^{2m+\alpha+n_{j}-1}\mathrm{e}^{-(u+\beta)s}ds$
$\displaystyle=\alpha(\beta+\sqrt{\beta^{2}-1})^{\alpha}\sum_{m=0}^{\infty}\frac{1}{m!\Gamma(\alpha+m+1)2^{2m+\alpha}}\frac{\Gamma(2m+\alpha+n_{j})}{\left(u+\beta\right)^{2m+\alpha+n_{j}}}$
$\displaystyle\frac{\alpha\left(\beta+\sqrt{\beta^{2}-1}\right)^{\alpha}}{2^{\alpha}(u+\beta)^{n_{j}+\alpha}}\frac{\Gamma(\alpha+n_{j})}{\Gamma(\alpha+1)}\
_{2}F_{1}\left(\frac{n_{j}+\alpha}{2},\frac{n_{j}+\alpha+1}{2};\alpha+1;\frac{1}{(u+\beta)^{2}}\right)$

∎

###  A.6 Proof of Convergence of the Finite DP

To prove convergence of the finite Dirichlet process to the infinite DP when
$\gamma=\alpha/\Lambda$ goes to infinity, we need the following results.

  * •

Let $\\{P_{t}\\}_{t=1,2\dots}$ be a sequence of species sampling processes
defined on the same space $\Theta$ as $P$. We denote by
$\pi_{t}(n_{1},\dots,n_{k})$ the eppf of $P_{t}$ for each $t=1,2,\dots$ and
$\pi_{0}(n_{1},\dots,n_{k})$ the eppf of $P_{0}$. Then, if

$\displaystyle a)\
\pi_{t}(n_{1},\dots,n_{k})=\int_{0}^{\infty}\tilde{\pi}_{t}(u;n_{1},\dots,n_{k})du,\quad\text{for
each }t$ $\displaystyle b)\
\pi_{0}(n_{1},\dots,n_{k})=\int_{0}^{\infty}\tilde{\pi}_{0}(u;n_{1},\dots,n_{k})du$
$\displaystyle c)\
\lim_{t\rightarrow\infty}\pi_{t}(u,n_{1},\dots,n_{k})=\pi_{0}(u,n_{1},\dots,n_{k})\text{
for each }u>0\ \text{and }n_{1},\dots,n_{k}$

then
$\lim_{t\rightarrow\infty}\pi_{t}(n_{1},\dots,n_{k})=\pi_{0}(n_{1},\dots,n_{k})\text{
for each }n_{1},\dots,n_{k}$, moreover this latter implies that $\\{P_{t}\\}$
converges in law to $P$ (see Argiento et al., 2016).

Now let $n_{1},\dots,n_{k}$ be a composition of $n$. the eppf of the Dirichlet
process is

$\pi_{0}(n_{1},\dots,n_{k})=\frac{\Gamma(\alpha)}{\Gamma(\alpha+n)}\prod_{j=1}^{k}\alpha\Gamma(n_{j})$
(A.33)

It is trivial to observe that

$\pi_{0}(n_{1},\dots,n_{k})=\int_{0}^{\infty}\tilde{\pi}_{0}(u;n_{1},\dots,n_{k})du=\int_{0}^{\infty}\frac{u^{n-1}(u+1)^{-(\alpha+n)}}{\Gamma(n)}\prod_{j=1}^{k}\alpha\Gamma(n_{j})du$
(A.34)

From Eq. (A.34) we see that condition $b)$ is satisfied. We now check that
also condition $a)$ is satisfied. Let $\pi_{\Lambda}(n_{1},\dots,n_{k})$ be
the eppf of a Finite DP with parameter $\Lambda$ and $\gamma=\alpha/\Lambda$.
From Eq. (A.33) we can easily derive

$\begin{split}\pi_{\Lambda}(n_{1},\dots,n_{k})&=\int_{0}^{\infty}\tilde{\pi}_{\Lambda}(u;n_{1},\dots,n_{k})du\\\
&=\int_{0}^{\infty}\frac{u^{n-1}}{\Gamma(n)}\frac{\Lambda+k(u+1)^{\alpha/\Lambda}}{\Lambda(u+1)^{\alpha(k+1)/\Lambda}(u+1)^{n}}\exp\left\\{-\Lambda\frac{(u+1)^{\alpha/\Lambda}-1}{(u+1)^{\alpha/\Lambda}}\right\\}\\\
&\hskip
227.62204pt\times\prod_{j=1}^{k}\Lambda\frac{\Gamma(\alpha/\Lambda+n_{j})}{\Gamma(\alpha/\Lambda)}du\end{split}$

We observe now that

  1. i.

$\lim_{\Lambda\rightarrow\infty}\frac{\Lambda+k(u+1)^{\alpha/\Lambda}}{\Lambda(u+1)^{\alpha(k+1)/\Lambda}(u+1)^{n}}=\frac{1}{(u+1)^{n}}$

  2. ii.

$\lim_{\Lambda\rightarrow\infty}\exp\left\\{-\Lambda\frac{(u+1)^{\alpha/\Lambda}-1}{(u+1)^{\alpha/\Lambda}}\right\\}=\frac{1}{(1+u)^{\alpha}}$

  3. iii.

For each $j=1,\dots,k$

$\lim_{\Lambda\rightarrow\infty}\Lambda\frac{\Gamma(\alpha/\Lambda+n_{j})}{\Gamma(\alpha/\Lambda)}=\alpha\Gamma(n_{j})$

From the previous three items, we conclude that

$\lim_{\Lambda\rightarrow\infty}\tilde{\pi}_{\Lambda}(u;n_{1},\dots,n_{k})=\tilde{\pi}_{0}(u;n_{1},\dots,n_{k})$

which proves that condition $c)$ is also verified.

##  Appendix B Appendix: The conditional Gibbs sampler for Finite Dirichlet
Mixture model

Here we describe the conditional algorithm for Finite Mixture models in the
particular case of Subsection 9.1.1. In what follows we also place a prior on
$\gamma$ and $\Lambda$ as implemented in the Example of Section 10. The full
model is:

$\displaystyle Y_{i}\mid c_{i},\tau_{1},\ldots,\tau_{M},M$ $\displaystyle\sim$
$\displaystyle f(y\mid\tau_{c_{i}})$ $\displaystyle\tau_{i}\mid M$
$\displaystyle\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}$ $\displaystyle
p_{0}(\tau)$ $\displaystyle\mathbf{w}\mid M$ $\displaystyle\sim$
$\displaystyle\text{Dirichlet}_{M}(\gamma,\cdots,\gamma)$ $\displaystyle
c_{i}\mid M,\mathbf{w}$ $\displaystyle\sim$
$\displaystyle\text{Multinomial}_{M}(1,w_{1},\ldots,w_{M})$
$\displaystyle\gamma$ $\displaystyle\sim$
$\displaystyle\text{Gamma}(a_{1},b_{1})$ $\displaystyle M$ $\displaystyle\sim$
$\displaystyle\text{Poisson}(\Lambda)$ $\displaystyle\Lambda$
$\displaystyle\sim$ $\displaystyle\text{Gamma}(a_{2},b_{2})$

We build a blocked Gibbs sampler to update blocks of parameters, which are
drawn from multivariate distributions. In particular, the parameters of
interest are $(\tilde{P},\mathbf{c},U)$, where $\tilde{P}$ is the unnormalized
finite point process and $U$ is an auxiliary variable introduced in Theorem 2.
Full conditionals can be derived for most of the parameters. The main steps of
the algorithms are:

  1. 1.

Sampling from $\bm{\mathcal{L}}(U|{\bm{Y}},\mathbf{c},\widetilde{P})$: by
construction, conditionally on $\tilde{P}$, the random variable $U$ is
distributed as Gamma with parameters $(n,T)$.

  2. 2.

Sampling from $\bm{\mathcal{L}}(c_{i}\mid u,{\bm{Y}},\tilde{P})$: each
$c_{i}$, for $i=1,\dots,n$, has a discrete law with support $\\{1,\dots,M\\}$,
and probabilities $\mathbb{P}(c_{i}=m)\propto S_{m}f(Y_{i};\tau_{m})$. Le
$k\leq M$ be the number of allocated components after resampling the entire
allocation vector $\mathbf{c}$. Rename the allocated components from 1 to $k$,
so that the allocated clusters correspond to the first $k$ components in the
mixture and the remaining $(M-k)$ are empty for $k<M$ .

  3. 3.

Sampling from $\bm{\mathcal{L}}(M,\mathbf{w}|u,{\bm{c}},{\bm{Y}})$: first we
observe that conditionally on ${\bm{c}}$, the weights $\mathbf{w}$ do not
depend on the observations $\bm{Y}$. Therefore, we have to sample from
$\mathcal{L}(M,\mathbf{w}\mid u,{\bm{c}})$. As stated in Theorem 2, we can
split this step into three sub-steps. Recall that we partition the vector
$\mathbf{w}=(\mathbf{w}_{a},\mathbf{w}_{na})$, where
$\mathbf{w}_{a}=(w_{1},\ldots,w_{k})$ and
$\mathbf{w}_{na})=(w_{k+1},\ldots,w_{M})$ correspond to the allocated and
unallocated components, respectively. Moreover, in our construction $k$ is
fixed and determined in Step 2 and $w_{j}\propto S_{j}$.

    1. 3.a

Sampling from $\bm{\mathcal{L}}(M\mid u,\mathbf{c},\bm{Y})$: Note that
$M=k+M^{(na)}$, where $M^{(na)}$ is the number of unallocated components. We
sample $M^{(na)}$ from the discrete probability measure defined
on$\\{0,1,2,\ldots\\}$:

$\displaystyle\Pr(M^{(na)}=m)$ $\displaystyle=$ $\displaystyle q^{\star}_{m}$
$\displaystyle=$
$\displaystyle\frac{(u+1)^{\gamma}k}{(u+1)^{\gamma}k+\Lambda}\mathcal{P}_{0}(\Lambda/(u+1)^{\gamma})+\frac{\Lambda}{(u+1)^{\gamma}k+\Lambda}\mathcal{P}_{1}(\Lambda/(u+1)^{\gamma})$

which corresponds to a two-components mixture. Here $\mathcal{P}_{i}(\lambda)$
denotes the shifted Poisson pmf on $\\{i,i+1,i+2,\ldots\\}$ with mean
$i+\lambda$, $i=0,1$.

    2. 3.a

Sampling from $\bm{\mathcal{L}}(\mathbf{w}^{(a)}|u,\mathbf{c},\bm{Y})$: the
allocated process is a set of independent pairs of variables
$\left\\{(S_{1},\tau_{1}),\dots(S_{k},\tau_{k})\right\\}$ such that, for
$m=1,\dots,k$:

$\displaystyle S_{m}$ $\displaystyle\sim$
$\displaystyle\text{Gamma}(n_{m}+\gamma,u+1)$
$\displaystyle\mathbb{P}(\tau_{m}=d\tau_{m}\mid\mathbf{c},\mathbf{y})$
$\displaystyle\propto$
$\displaystyle\prod_{\\{i:c_{i}=m\\}}f(y_{i}\mid\tau_{m})p_{0}(\tau_{m})d\tau_{m}$

where $n_{m}$ is the cardinality of the set $\\{i:c_{i}=m\\}$, i.e is the
number of allocation variables $c_{i}$ such that $c_{i}=m$. When the density
$p_{0}$ and the family of kernels $f(y\mid\tau)$ are conjugate, then all the
full conditionals $\mathbb{P}(\tau_{m}=d\tau_{m}\mid\mathbf{c},\mathbf{y})$
are available in closed form, simplifying the implementation of the algorithm.

    3. 3.c

Sampling from $\bm{\mathcal{L}}(\mathbf{w}^{(na)}\mid
u,\mathbf{c},\bm{Y},M_{na})$: the weights corresponding to unallocated jumps
can be sampled as follows

$\displaystyle S_{m}$ $\displaystyle\sim$
$\displaystyle\text{Gamma}(\gamma,u+1)$
$\displaystyle\tau_{m}\mid\mathbf{c},\mathbf{y}$ $\displaystyle\sim$
$\displaystyle p_{0}(\tau_{m})$

for $m=k+1,\ldots,M$.

  4. 4.

Sampling from $\bm{\mathcal{L}}(\Lambda\mid u,k,\gamma)$: this steps involves
a conjugate update:

$\displaystyle\Lambda\mid\text{rest}$ $\displaystyle\sim$
$\displaystyle\frac{\psi(u)}{1+b_{2}}\text{Gamma}(a_{2}^{\star}+1,1-\psi(u)+b_{2})+\frac{1-\psi(u)+b_{2}}{1+b_{2}}\text{Gamma}(a_{2}^{\star},1-\psi(u)+b_{2})$
$\displaystyle a_{2}^{\star}$ $\displaystyle=$ $\displaystyle k+a_{2}$
$\displaystyle\psi(u)$ $\displaystyle=$
$\displaystyle\frac{1}{(u+1)^{\gamma}}$

where $\psi(u)$ is the Laplace transform of a Gamma$(\gamma,1)$ density.

  5. 5.

Sampling from $\bm{\mathcal{L}}(\gamma\mid u,k,\mathbf{c},\Lambda)$: we have
to implement a Metropolis-Hasting step to sample from:

$\pi(\gamma\mid\text{rest})\propto\left(\Lambda\psi(u)+k\right)\mathrm{e}^{\Lambda\psi(u)}\frac{1}{\psi(u)^{k}}\prod_{j=1}^{k}\frac{\Gamma(\gamma+n_{j})}{\Gamma(\gamma)}$

as the full conditional is not available in closed form. For the example of
Section 10, we have used a random walk as proposal distribution. Adaptive
strategies can be easily implemented to improve mixing.

We note that we could have opted for a Negative-Binomial distribution as prior
on $M$, with a Beta hyper-prior on the probability of success. The Negative-
Binomial can be used to induce more sparsity. In this case step (3.b),(4) and
(5) need to modified accordingly.

##  Appendix C Appendix: Additional Figures

Figure C.8: Posterior distribution of $k$, when $\lambda$ and $\gamma$ are
random. Figure C.9: Posterior distribution of $M^{(na)}$, when $\lambda$ and
$\gamma$ are random.

## References

  * Aldous (1985) Aldous, D. J. (1985) Exchangeability and related topics.  In École d’Été de Probabilités de Saint-Flour XIII-1983, 1–198. Springer.
  * Argiento et al. (2016) Argiento, R., Bianchini, I. and Guglielmi, A. (2016) Posterior sampling from $\varepsilon$-approximation of normalized completely random measure mixtures.  Electron. J. Statist., 10, 3516–3547.  URL: https://doi.org/10.1214/16-EJS1168.
  * Argiento et al. (2019) Argiento, R., Cremaschi, A. and Vannucci, M. (2019) Hierarchical normalized completely random measures to cluster grouped data.  Journal of the American Statistical Association, to appear.  URL: https://www.carloalberto.org/research/working-papers/.
  * Biernacki et al. (2000) Biernacki, C., Celeux, G. and Govaert, G. (2000) Assessing a mixture model for clustering with the integrated completed likelihood.  IEEE transactions on pattern analysis and machine intelligence, 22, 719–725.
  * Callens et al. (2011) Callens, T., Galbusera, P., Matthysen, E., Durand, E. Y., Githiru, M., Huyghe, J. R. and Lens, L. (2011) Genetic signature of population fragmentation varies with mobility in seven bird species of a fragmented kenyan cloud forest.  Molecular Ecology, 20, 1829–1844.
  * Charalambides (2005) Charalambides, C. A. (2005) Combinatorial methods in discrete distributions, vol. 600.  John Wiley & Sons.
  * Daley and Vere-Jones (2007) Daley, D. J. and Vere-Jones, D. (2007) An introduction to the theory of point processes: volume II: general theory and structure.  Springer Science & Business Media.
  * Dellaportas and Papageorgiou (2006) Dellaportas, P. and Papageorgiou, I. (2006) Multivariate mixtures of normals with unknown number of components.  Statistics and Computing, 16, 57–68.
  * DeVore and Lorentz (1993) DeVore, R. A. and Lorentz, G. G. (1993) Constructive approximation, vol. 303.  Springer Science & Business Media.
  * Devroye (2009) Devroye, L. (2009) Random variate generation for exponentially and polynomially tilted stable distributions.  ACM Transactions on Modeling and Computer Simulation (TOMACS), 19, 18.
  * Dey et al. (2012) Dey, D. D., MüIler, P. and Sinha, D. (2012) Practical nonparametric and semiparametric Bayesian statistics, vol. 133.  Springer Science & Business Media.
  * Eddelbuettel and François (2011) Eddelbuettel, D. and François, R. (2011) Rcpp: Seamless R and C++ integration.  Journal of Statistical Software, 40, 1–18.  URL: http://www.jstatsoft.org/v40/i08/.
  * Erdélyi et al. (1953) Erdélyi, A., Magnus, W., Oberhettinger, F., Tricomi, F. G. and Bateman, H. (1953) Higher transcendental functions, vol. 2.  McGraw-Hill New York.
  * Falush et al. (2003) Falush, D., Stephens, M. and Pritchard, J. K. (2003) Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies.  Genetics, 164, 1567–1587.
  * Favaro et al. (2015) Favaro, S., Nipoti, B., Teh, Y. W. et al. (2015) Random variate generation for Laguerre-type exponentially tilted $\alpha-$stable distributions.  Electronic Journal of Statistics, 9, 1230–1242.
  * Favaro and Teh (2013) Favaro, S. and Teh, Y. W. (2013) MCMC for normalized random measure mixture models.  Statistical Science, 335–359.
  * Feller (1971) Feller, W. (1971) An Introduction to Probability Theory and Its Applications, vol. II.  John Wiley, New York, second edition edn.
  * Ferguson (1973) Ferguson, T. S. (1973) A bayesian analysis of some nonparametric problems.  The annals of statistics, 209–230.
  * Frühwirth-Schnatter (2006) Frühwirth-Schnatter, S. (2006) Finite mixture and Markov switching models.  Springer Science & Business Media.
  * Fruhwirth-Schnatter et al. (2019) Fruhwirth-Schnatter, S., Celeux, G. and Robert, C. P. (2019) Handbook of mixture analysis.  Chapman and Hall/CRC.
  * Frühwirth-Schnatter and Malsiner-Walli (2018) Frühwirth-Schnatter, S. and Malsiner-Walli, G. (2018) From here to infinity: sparse finite versus dirichlet process mixtures in model-based clustering.  Advances in Data Analysis and Classification.  URL: https://doi.org/10.1007/s11634-018-0329-y.
  * Galbusera et al. (2000) Galbusera, P., Lens, L., Schenck, T., Waiyaki, E. and Matthysen, E. (2000) Genetic variability and gene flow in the globally, critically-endangered taita thrush.  Conservation Genetics, 1, 45–55.
  * Geisser and Eddy (1979) Geisser, S. and Eddy, W. F. (1979) A predictive approach to model selection.  Journal of the American Statistical Association, 74, 153–160.
  * Gelfand and Kottas (2002) Gelfand, A. E. and Kottas, A. (2002) A computational approach for full nonparametric bayesian inference under dirichlet process mixture models.  Journal of Computational and Graphical Statistics, 11, 289–305.
  * Goutis and Robert (1998) Goutis, C. and Robert, C. P. (1998) Model choice in generalised linear models: A bayesian approach via kullback-leibler projections.  Biometrika, 85, 29–37.
  * Gradshteyn and Ryzhik (2007) Gradshteyn, I. and Ryzhik, L. (2007) Table of integrals, series, and products - Seventh Edition.  San Diego (USA): Academic Press, sixth edn.
  * Green (1995) Green, P. J. (1995) Reversible jump markov chain monte carlo computation and bayesian model determination.  Biometrika, 82, 711–732.
  * Hofert (2011) Hofert, M. (2011) Sampling exponentially tilted stable distributions.  ACM Transactions on Modeling and Computer Simulation (TOMACS), 22, 3.
  * Huelsenbeck and Andolfatto (2007) Huelsenbeck, J. P. and Andolfatto, P. (2007) Inference of population structure under a dirichlet process model.  Genetics, 175, 1787–1802.
  * Huelsenbeck et al. (2011) Huelsenbeck, J. P., Andolfatto, P. and Huelsenbeck, E. T. (2011) Structurama: Bayesian inference of population structure.  Evolutionary Bioinformatics, 7, EBO–S6761.
  * Ishwaran and James (2001) Ishwaran, H. and James, L. F. (2001) Gibbs sampling methods for stick-breaking priors.  Journal of the American Statistical Association, 96, 161–173.
  * Ishwaran and James (2003) — (2003) Some further developments for stick-breaking priors: finite and infinite clustering and classification.  Sankhyā: The Indian Journal of Statistics, 577–592.
  * Ishwaran and Zarepour (2002) Ishwaran, H. and Zarepour, M. (2002) Exact and approximate sum representations for the dirichlet process.  Canadian Journal of Statistics, 30, 269–283.
  * Jacod and Shiryaev (2013) Jacod, J. and Shiryaev, A. (2013) Limit theorems for stochastic processes, vol. 288.  Springer Science & Business Media.
  * James et al. (2009) James, L. F., Lijoi, A. and Prünster, I. (2009) Posterior analysis for normalized random measures with independent increments.  Scandinavian Journal of Statistics, 36, 76–97.
  * Jordan (2010) Jordan, M. I. (2010) Hierarchical models, nested models and completely random measures.  Frontiers of statistical decision making and Bayesian analysis: In honor of James O. Berger. New York: Springer, 207–218.
  * Kalli et al. (2011) Kalli, M., Griffin, J. E. and Walker, S. G. (2011) Slice sampling mixture models.  Statistics and computing, 21, 93–105.
  * Kim et al. (2006) Kim, S., Tadesse, M. G. and Vannucci, M. (2006) Variable selection in clustering via dirichlet process mixture models.  Biometrika, 93, 877–893.
  * Komárek (2009) Komárek, A. (2009) A new r package for bayesian estimation of multivariate normal mixtures allowing for selection of the number of components and interval-censored data.  Computational Statistics & Data Analysis, 53, 3932–3947.
  * Kong (1992) Kong, A. (1992) A note on importance sampling using standardized weights.  University of Chicago, Dept. of Statistics, Tech. Rep, 348.
  * Lau and Green (2007) Lau, J. W. and Green, P. J. (2007) Bayesian model-based clustering procedures.  Journal of Computational and Graphical Statistics, 16, 526–558.
  * Lijoi et al. (2007) Lijoi, A., Mena, R. H. and Prünster, I. (2007) Controlling the reinforcement in bayesian non-parametric mixture models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69, 715–740.
  * Lijoi et al. (2010) Lijoi, A., Prünster, I. et al. (2010) Models beyond the dirichlet process.  Bayesian nonparametrics, 28, 3.
  * Lo (1984) Lo, A. Y. (1984) On a class of bayesian nonparametric estimates: I. density estimates.  The Annals of Statistics, 351–357.
  * Malsiner-Walli et al. (2016) Malsiner-Walli, G., Frühwirth-Schnatter, S. and Grün, B. (2016) Model-based clustering based on sparse finite gaussian mixtures.  Statistics and computing, 26, 303–324.
  * Malsiner-Walli et al. (2017) — (2017) Identifying mixtures of mixtures using bayesian estimation.  Journal of Computational and Graphical Statistics, 26, 285–295.
  * McLachlan et al. (2000) McLachlan, G. J., Lee, S. X. and Rathnayake, S. I. (2000) Finite mixture models.  Annual Review of Statistics and Its Application.
  * Miller and Harrison (2013) Miller, J. W. and Harrison, M. T. (2013) A simple example of dirichlet process mixture inconsistency for the number of components.  In Advances in neural information processing systems, 199–206.
  * Miller and Harrison (2018) — (2018) Mixture models with a prior on the number of components.  Journal of the American Statistical Association, 113, 340–356.
  * Møller and Waagepetersen (2003) Møller, J. and Waagepetersen, R. P. (2003) Statistical inference and simulation for spatial point processes.  Chapman and Hall/CRC.
  * Muliere and Tardella (1998) Muliere, P. and Tardella, L. (1998) Approximating distributions of random functionals of ferguson-dirichlet priors.  Canadian Journal of Statistics, 26, 283–297.
  * Neal (2000) Neal, R. M. (2000) Markov chain sampling methods for dirichlet process mixture models.  Journal of Computational and Graphical Statistics, 9, 249–265.
  * Nobile (1994) Nobile, A. (1994) Bayesian Analysis of Finite Mixture Distributions.  Ph.D. thesis, Department of Statistics, Carnegie Mellon University.
  * Nobile et al. (2004) Nobile, A. et al. (2004) On the posterior distribution of the number of components in a finite mixture.  The Annals of Statistics, 32, 2044–2073.
  * Pitman (1996) Pitman, J. (1996) Blackwell-macqueen urn scheme.  Statistics, Probability, and Game Theory: Papers in Honor of David Blackwell, 30, 245.
  * Pitman (2003) — (2003) Poisson-Kingman Partitions.  In Science and Statistics: a Festschrift for Terry Speed, vol. 40 of IMS Lecture Notes-Monograph Series, 1–34. Hayward (USA): Institute of Mathematical Statistics.
  * Pitman (2006) — (2006) Combinatorial stochastic processes.  In École d’Été de Probabilités de Saint-Flour XXXII-2002, 1–255. Springer.
  * Pollard (1946) Pollard, H. (1946) The representation of $\mathrm{e}^{-x^{\lambda}}$ as Laplace integral.  Bulletin of the American Mathematical Society, 52, 908–910.
  * Pritchard et al. (2000) Pritchard, J. K., Stephens, M. and Donnelly, P. (2000) Inference of population structure using multilocus genotype data.  Genetics, 155, 945–959.
  * Pritchard and Wen (2003) Pritchard, J. K. and Wen, W. (2003) Documentation for STRUCTURE software: Version 2.3.X.  Available at https://web.stanford.edu/group/pritchardlab/structure.html.
  * Regazzini et al. (2003) Regazzini, E., Lijoi, A. and Prünster, I. (2003) Distributional results for means of normalized random measures with independent increments.  Annals of Statistics, 560–585.
  * Richardson and Green (1997) Richardson, S. and Green, P. J. (1997) On bayesian analysis of mixtures with an unknown number of components.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59, 731–792.
  * Roeder (1990) Roeder, K. (1990) Density estimation with confidence sets exemplified by superclusters and voids in the galaxies.  Journal of the American Statistical Association, 85, 617–624.
  * Rousseau and Mengersen (2011) Rousseau, J. and Mengersen, K. (2011) Asymptotic behaviour of the posterior distribution in overfitted mixture models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73, 689–710.
  * Sokal (1997) Sokal, A. (1997) Monte carlo methods in statistical mechanics: foundations and new algorithms.  In Functional integration, 131–192. Springer.
  * Stephens (2000) Stephens, M. (2000) Bayesian analysis of mixture models with an unknown number of components-an alternative to reversible jump methods.  Annals of statistics, 40–74.
  * Tadesse et al. (2005) Tadesse, M. G., Sha, N. and Vannucci, M. (2005) Bayesian variable selection in clustering high-dimensional data.  Journal of the American Statistical Association, 100, 602–617.
  * Wiper et al. (2001) Wiper, M., Insua, D. R. and Ruggeri, F. (2001) Mixtures of gamma distributions with applications.  Journal of Computational and Graphical Statistics, 10, 440–454.