capbtabboxtable[][]

# Nonparametric priors with
full-range borrowing of information

Filippo Ascolani  Bocconi Institute for Data Science and Analytics (BIDSA),
Bocconi University Milano, Italy  Beatrice Franzolini  Bocconi Institute for
Data Science and Analytics (BIDSA),
Bocconi University Milano, Italy  Antonio Lijoi  Bocconi Institute for Data
Science and Analytics (BIDSA),
Bocconi University Milano, Italy  Igor Prünster  Bocconi Institute for Data
Science and Analytics (BIDSA),
Bocconi University Milano, Italy

###### Abstract

Modeling of the dependence structure across heterogeneous data is crucial for
Bayesian inference since it directly impacts the borrowing of information.
Despite the extensive advances over the last two decades, most available
proposals allow only for non–negative correlations. We derive a new class of
dependent nonparametric priors that can induce correlations of any sign, thus
introducing a new and more flexible idea of borrowing of information. This is
achieved thanks to a novel concept, which we term hyper-tie, and represents a
direct and simple measure of dependence. We investigate prior and posterior
distributional properties of the model and develop algorithms to perform
posterior inference. Illustrative examples on simulated and real data show
that our proposal outperforms alternatives in terms of prediction and
clustering.

Keywords: Bayesian nonparametrics; Borrowing of information; Completely random
measure; Dependent nonparametric prior; Negative correlation; Partial
exchangeability.

##  1 Introduction

Bayesian nonparametric methods are increasingly popular, mainly thanks to
their flexibility and strong foundations. The most common assumption
underlying Bayesian models is exchangeability, which corresponds to invariance
of the joint distribution of the observations with respect to finite
permutations. However, real phenomena often present a level of heterogeneity
that makes exchangeability unrealistic: collected data may refer to different
features, populations, or, in general, may be collected under different
experimental conditions. Such situations entail a significant level of
heterogeneity and opportunities for borrowing information, that can be
exploited through the notion of partial exchangeability, which implies
exchangeability within each experimental condition, but not across. Two
sequences of observations ${X}=(X_{i})_{i\geq 1}$ and $Y=(Y_{j})_{j\geq 1}$,
taking values in a space $\mathbb{X}$, are partially exchangeable if and only
if, for all sample sizes $(n,m)$ and all permutations $(\pi_{1},\pi_{2})$,

$\bigl{(}(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m}\bigr{)}\overset{d}{=}\bigl{(}(X_{\pi_{1}(i)})_{i=1}^{n},(Y_{\pi_{2}(j)})_{j=1}^{m}\bigr{)}.$

with $\overset{d}{=}$ denoting equality in distribution. From an inferential
point of view, partial exchangeability entails that the order of the
observations within each sample is non-informative, while the belonging to a
specific sample is relevant and has to be taken into account. Moreover, by de
Finetti’s representation theorem (de Finetti, 1938) $X$ and $Y$ are partially
exchangeable if and only if there exist random probabilities
$(\tilde{p}_{1},\tilde{p}_{2})$ such that for any $i,j=1,\ldots,n$

$\displaystyle(X_{i},Y_{j})\mid(\tilde{p}_{1},\tilde{p}_{2})\overset{iid}{\sim}\tilde{p}_{1}\times\tilde{p}_{2}\qquad(\tilde{p}_{1},\,\tilde{p}_{2})\sim
Q$ (1)

with $Q$ playing the role of the prior. The dependence induced by $Q$ at the
level of the observables defines the Bayesian learning mechanism and it
connects to the notion of borrowing of information. This term was first coined
by John Tukey (Brillinger, 2002) and popularized with reference to Stein’s
paradox and empirical Bayes techniques in Efron and Morris (1977). More
generally, statisticians refer to borrowing of information when many samples
contribute to inference related to just one sample. Imagine collecting the
samples $(X_{i})_{i=1}^{n}$ and $(Y_{j})_{j=1}^{m}$, while being interested
only in the parameter $\tilde{p}_{1}$ associated to $X$. The simplest approach
could be to disregard the second sample $(Y_{j})_{j=1}^{m}$, with the drawback
of losing potentially useful information. The typical borrowing instead
consists in shrinking the estimates for different samples towards each other:
shrinkage is justified by the fact that distributions of different, but
related, populations are expected to be similar in terms of shape and/or
location. However, many contexts may still require borrowing of information
between $(X_{i})_{i=1}^{n}$ and $(Y_{j})_{j=1}^{m}$, but without necessarily
resulting in shrinkage. Indeed, one’s available prior information may imply
that the responses in different groups have a negative association and, thus,
tend to be dissimilar in location, which makes shrinkage undesirable.
Similarly, when there is no pre-experimental knowledge on the dependence
between $X_{i}$ and $Y_{j}$, a flexible prior specification allowing also for
negative association would be more appropriate. A toy parametric example to
further clarify that borrowing does not necessarily imply classic shrinkage is
provided in Section S2. of the Appendix. Some applied scenarios of borrowing
of information not resulting in shrinkage are, for instance, the study of
survival times and abundances of competitive species (Lee et al., 2020), the
incorporation of retrospective data to study associations between biomarkers
(Gong et al., 2021), the association between dental caries and dental
fluorosis (Lorenz et al., 2018), the analysis of stocks and bonds returns (see
Bhardwaj and Dunsby, 2013, and Section 6.3), and the clustering of
multivariate responses with missing entries (see Section 6.4). In this paper
we introduce a class of nonparametric priors that allows for a more general
version of borrowing, which includes shrinkage as a special case. These can be
used as core building blocks for models tailored to specific applications.

Starting from the pioneering works of Cifarelli and Regazzini (1978) and
MacEachern (1999, 2000), Bayesian nonparametric contributions for
non–exchangeable data have grown substantially, see Foti and Williamson
(2013), Müller et al. (2015) and Quintana et al. (2022) for insightful
reviews. The vast majority of nonparametric models for partially exchangeable
data entails that the random probabilities in (13) are such that

$\begin{cases}\tilde{p}_{1}\overset{a.s.}{=}\sum_{k\geq
1}\bar{J}_{k}\delta_{\theta_{k}}\\\\[6.0pt]
\tilde{p}_{2}\overset{a.s.}{=}\sum_{k\geq
1}\bar{W}_{k}\delta_{\phi_{k}}\end{cases}\quad\quad\theta_{k}\overset{\text{i.i.d.}}{\sim}P_{0},\quad\phi_{k}\overset{\text{i.i.d.}}{\sim}P_{0}$
(2)

where the random weights $\left((\bar{J}_{k}),(\bar{W}_{k})\right)$ and the
atoms $\left((\theta_{k}),(\phi_{k})\right)$ are independent and
$\theta_{k}\perp\phi_{h}$ for $k\neq h$. In this paper we focus on this class
of models and, for ease of exposition, take $\tilde{p}_{1}$ and
$\tilde{p}_{2}$ with the same marginal distribution.

A first prominent strategy for defining $Q$ is to explicitly assign the
distribution of the weights and the atoms in (2) so to create dependence
between $\tilde{p}_{1}$ and $\tilde{p}_{2}$: this approach has led to
dependent Dirichlet processes (MacEachern, 1999, 2000; Quintana et al., 2022),
dependent stick-breaking processes, kernel stick-breaking processes (Dunson
and Park, 2008), probit stick-breaking processes (Rodriguez and Dunson, 2011)
and others. Despite their flexibility and the availability of posterior
sampling schemes, the derivation of analytical results is very difficult for
these models; it is often not clear how the dependence of the series reflects
at the level of the observables and therefore such methods may lack
transparency.

A second popular strategy, analytically more tractable, relies on completely
random measures (CRMs) either working directly on the law of multi-dimensional
vectors of CRMs (Epifani and Lijoi, 2010; Griffin and Leisen, 2017; Riva-
Palacio and Leisen, 2021) or combining conditionally independent CRMs, using
additive structures (Müller et al., 2004; Griffin et al., 2013; Lijoi and
Nipoti, 2014; Lijoi et al., 2014a, b), nested structures (Rodriguez et al.,
2008; Camerlenghi et al., 2019), or hierarchical structures (Teh et al., 2006;
Camerlenghi et al., 2019). CRMs are then suitably transformed to obtain the
random probabilities in (2).

Dependent random probabilities clearly induce dependence across groups of
observations. The simplest and most intuitive way to quantify the dependence
structure is through correlations. Therefore, when considering correlations
among observables, we will implicitly assume real-valued $X_{i}$’s and
$Y_{j}$’s, namely $\mathbb{X}=\mathbb{R}$. All other results and concepts are
valid for general spaces $\mathbb{X}$. A first result in this direction shows
that, regardless of the specific dependent model, observations in different
groups cannot be more correlated (in absolute sense) than the ones in the same
group.

######  Proposition 1.

Suppose $X$ and $Y$ are partially exchangeable sequences, such that
$\tilde{p}_{1}$ and $\tilde{p}_{2}$ in (13) have the same marginal
distribution. Then

$-\mbox{\rm corr}(X_{i},X_{i^{\prime}})\leq\mbox{\rm
corr}(X_{i},Y_{j})\leq\mbox{\rm corr}(X_{i},X_{i^{\prime}}),$

for any $i,i^{\prime}$ and $j$.

Due to exchangeability within each group, the upper bound in Proposition 1 is
always non–negative and it can be shown that, for all the models as in (2),
the correlation between observations in the same sample, $\mbox{\rm
corr}(X_{i},X_{i^{\prime}})$, is determined by the probability of a tie. As
for the correlation across samples $\mbox{\rm corr}(X_{i},Y_{j})$, we show
that a similar result holds true, with _hyper-ties_ , the new notion we
introduce, replacing ties.

Moreover, note that for most models based on CRMs, which allow for the
computation of the correlation, $\mbox{\rm corr}(X_{i},Y_{j})$ turns out to be
positive: this happens in particular when the interaction between two or more
groups is of interest. Therefore, the literature available to date within the
partially exchangeable setting is focused on models that attain a limited
range of possible values of the correlation, when it can be evaluated. Here we
aim to overcome this limitation and introduce a novel class of priors which
yield a wider range of correlation values among the observables, including
those with negative sign. The next result shows that the sign of the
correlation is only determined by the dependence structure between the atoms.

######  Proposition 2.

Suppose $X$ and $Y$ are partially exchangeable sequences, such that the
underlying $\tilde{p}_{1}$ and $\tilde{p}_{2}$ are as in (2). Moreover, for
any $k$ and $k^{\prime}$, let $\mbox{\rm
corr}(\theta_{k},\phi_{k^{\prime}})\geq 0$. Then $\mbox{\rm
corr}\,(X_{i},Y_{j})\geq 0$, for any $i$ and $j$.

For instance, hierarchical processes (Teh et al., 2006; Camerlenghi et al.,
2019), which represent one of the most popular dependent models, induce
dependence by the sharing of atoms across groups. However, by Proposition 2,
this means that achieving negative correlation is impossible. Hence, a
flexible joint distribution for the sequence of atoms must be specified. This
task is accomplished by our proposal, termed normalized CRMs with Full-Range
Borrowing of Information (n-FuRBI), that allows to attain any possible value
for the correlation specified in Proposition 1. Moreover, it encompasses many
previous constructions as special cases. We will show that it nicely combines
the flexibility of the random series construction with the analytical
tractability featured by CRMs. Our proposal allows to consider any interesting
choice of borrowing of information: independence, classical shrinkage, but
also repulsion of estimates for different samples, generating what we term
full–range borrowing of information. Note that the repulsive behaviour of
n-FuRBI is different from the one featured by the priors introduced in
Petralia et al. (2012) and Quinlan et al. (2017), that induce repulsion among
the atoms of a single random probability measure.

The appendix includes all the analytical derivations and proofs, the
simulation algorithms for the implementation of the proposed class of models,
additional examples and numerical studies. In the following we use the prefix
S to indicate sections of the Appendix. The code to allow full replication of
the numerical results is available at
https://github.com/beatricefranzolini/FuRBI.

##  2 General results on dependent processes

The vast majority of dependent processes introduced in the literature are
almost surely discrete and therefore admit a series representation as in (2).
A key preliminary step leading to the definition of hyper-tie and n-FuRBI
priors is the observation that the random probabilities in (2) can be embedded
into

$\begin{cases}p_{1}\overset{a.s.}{=}\sum_{k\geq
1}\bar{J}_{k}\delta_{(\theta_{k},\phi_{k})}\\\
p_{2}\overset{a.s.}{=}\sum_{k\geq
1}\bar{W}_{k}\delta_{(\theta_{k},\phi_{k})}\end{cases}\quad(\theta_{k},\phi_{k})\overset{i.i.d.}{\sim}G_{0},$
(3)

with $G_{0}$ a probability distribution on $\mathbb{X}\times\mathbb{X}$, whose
marginals equal $P_{0}$. While $p_{1}$ and $p_{2}$ share the same atoms, the
weights and the atoms are independent and the pair of random probability
measures $\tilde{p}_{1}$ and $\tilde{p}_{2}$ in (2) are obtained as the
projections over different coordinates of $p_{1}$ and $p_{2}$, namely
$\tilde{p}_{1}(\cdot)=p_{1}(\cdot\times\mathbb{X})$ and
$\tilde{p}_{2}(\cdot)=p_{2}(\mathbb{X}\times\cdot)$. The structure of popular
models is recovered by letting either $G_{0}=P_{0}^{2}$, which corresponds to
independence, or
$G_{0}(\mathrm{d}\theta,\mathrm{d}\phi)=P_{0}(\mathrm{d}\theta)\delta_{\\{\theta\\}}(\mathrm{d}\phi)$,
that is $\theta_{k}=\phi_{k}$ for any $k$ as happens for, e.g., hierarchical
processes (see Camerlenghi et al., 2019). Almost sure discreteness implies
that a sample from the random probability measure $\tilde{p}_{1}$ (or
$\tilde{p}_{2}$) will display ties with positive probability. The probability
of a tie, i.e. a coincidence of any two observations $i$ and $j$ in the same
sample, is

$\beta:=\mathbb{P}(X_{i}=X_{j})=\sum_{k\geq
1}\text{E}(\bar{J}^{2}_{k})=\sum_{k\geq
1}\text{E}(\bar{W}^{2}_{k})=\mathbb{P}(Y_{i}=Y_{j})$ (4)

with $(\bar{J}_{k})_{k\geq 1}$ and $(\bar{W}_{k})_{k\geq 1}$ equal in
distribution since we are assuming, for simplicity, that $\tilde{p}_{1}$ and
$\tilde{p}_{2}$ are equal in distribution. When considering jointly the two
samples, the concept of tie can be replaced by the one of _hyper-tie_ , that
is two observations in different samples coinciding with components having the
same label. According to (13), its probability is

$\gamma:=\sum_{k\geq
1}\mathbb{P}(X_{i}=\theta_{k},\,Y_{j}=\phi_{k})=\sum_{k\geq
1}\text{E}(\bar{J}_{k}\bar{W}_{k}).$ (5)

Sampling from components with the same label is equivalent to sampling the
same atom at the level of the underlying $(p_{1},p_{2})$ in (3). Clearly, when
the atoms are shared between $\tilde{p}_{1}$ and $\tilde{p}_{2}$, i.e.
$G_{0}(\mathrm{d}\theta,\mathrm{d}\phi)=P_{0}(\mathrm{d}\theta)\delta_{\\{\theta\\}}(\mathrm{d}\phi)$,
a hyper-tie corresponds to an actual tie between observations in different
samples.

The next result shows the relationship between $\beta$ and $\gamma$, the
probabilities of a tie and hyper-tie, respectively: in particular, the
probability of a tie is always larger and equality is attained if and only if
the probability masses of $p_{1}$ and $p_{2}$ are almost surely equal.

######  Proposition 3.

Let $(\tilde{p}_{1},\tilde{p}_{2})$ be as in (2) and $\beta,\gamma$ as in (4)
and (5), respectively. Then $0\leq\gamma\leq\beta$ and $\beta=\gamma$ if and
only if $\bar{W}_{k}\overset{a.s.}{=}\bar{J}_{k}$ for any $k$.

Hyper-ties play a crucial role in determining the dependence between
observables across groups, as the ties do for the dependence between
observables within groups, as shown by the next proposition.

######  Proposition 4.

Consider model (13) with $(\tilde{p}_{1},\tilde{p}_{2})$ as in (2). Then, for
any $i\neq i^{\prime}$ and any $j\neq j^{\prime}$

$\mbox{\rm corr}(X_{i},X_{i^{\prime}})=\mbox{\rm
corr}(Y_{j},Y_{j^{\prime}})=\beta\qquad\qquad\mbox{\rm
corr}(X_{i},Y_{j})=\gamma\,\rho_{0}$

with $\rho_{0}$ the correlation between two random variables jointly sampled
from $G_{0}$.

Thus, while the correlation between observations in the same sample equals the
probability of a tie, the correlation between observations from different
samples is determined by the probability of a hyper-tie, corrected by the
correlation between atoms. Clearly a suitable choice of the joint distribution
of the atoms makes the latter negative. Thus, by choosing $G_{0}$
appropriately, for instance as a bivariate normal, it is easy to tune the
correlation according to the available prior knowledge. The following
Corollary shows the values that can be attained, once the marginal law is
specified.

######  Corollary 1.

Consider model (13) with $(\tilde{p}_{1},\tilde{p}_{2})$ as in (2). If the
marginal distribution of $\tilde{p}_{1}$ and $\tilde{p}_{2}$ is fixed, then
$\mbox{\rm corr}(X_{i},Y_{j})\in[-\beta,\beta]$ and the extreme values are
attained if and only if the jumps are equal and $\rho_{0}=\pm 1$.

Unsurprisingly, with equal weights and jumps, which corresponds to full
exchangeability, one achieves the extreme case of $\mbox{\rm
corr}(X_{i},Y_{j})=\beta$. Null correlation, instead, is attained when atoms
are uncorrelated or when the probability of hyper-ties is zero. Lastly,
maximum negative correlation $\mbox{\rm corr}(X_{i},Y_{j})=-\beta$ is attained
with equal weights and negatively correlated atoms and can be thought of as
the opposite case with respect to exchangeability, at least in terms of
correlation. Ties and hyper-ties play a similar role also in the predictive
structure, as the next result shows.

######  Proposition 5.

Consider model (13) with $(\tilde{p}_{1},\tilde{p}_{2})$ as in (2). Then

$\mathbb{P}\left(X_{1}\in A,X_{2}\in B\right)=\beta P_{0}(A\cap
B)+\left(1-\beta\right)P_{0}(A)P_{0}(B).$

and

$\mathbb{P}\left(X_{1}\in A,Y_{1}\in B\right)=\gamma G_{0}(A\times
B)+\left(1-\gamma\right)P_{0}(A)P_{0}(B).$

The result is indeed quite intuitive. If $X_{1}$ and $Y_{1}$ form a hyper-tie
(with probability $\gamma$) they come from the same pair of atoms and need to
be sampled jointly; otherwise they refer to different atoms and are sampled
independently. The same happens inside each group, where $X_{1}$ and $X_{2}$
are equal with probability $\beta$.

######  Example 1.

The hierarchical Dirichlet process (Teh et al., 2006) is characterized by the
hierarchical representation
$\tilde{p}_{i}\mid\tilde{p}_{0}\overset{\text{i.i.d.}}{\sim}\text{DP}(\theta,\tilde{p}_{0})$,
with $\tilde{p}_{0}\sim\text{DP}(\theta_{0},P_{0})$, where $P_{0}$ is a
diffuse measure and DP$(\alpha,H)$ denotes the law of a Dirichlet process with
concentration parameter $\alpha>0$ and baseline distribution $H$. Since the
$\tilde{p}_{i}$’s share the atoms, an hyper-tie corresponds to an actual tie
between observations in different samples, so that with simple computations we
get

$\beta=\mbox{\rm
corr}(X_{i},X_{j})=1-\frac{\theta\theta_{0}}{(1+\theta)(1+\theta_{0})},\qquad\gamma=\mbox{\rm
corr}(X_{i},Y_{j})=\frac{1}{1+\theta_{0}}.$

Thus, the correlation among the observables is forced to be positive, with
$\theta_{0}$ tuning the dependence; see Example $1$ in Camerlenghi et al.
(2019) for more details.

Given the above results and considerations, it should be clear that $\gamma$
defined in (5) is crucial for tuning the level of dependence. However, closed
form expressions of $\gamma$ are available only for a few cases and, in fact,
we are facing a trade–off: on the one hand we have dependent processes based
on the stick-breaking representation, that allow for high flexibility while
sacrificing the availability of analytical results; on the other hand we have
constructions based on CRMs, for which an extensive theory has been developed,
though they are not as effective for tuning the dependence, since all the
existing instances produce non-negative correlation across samples. In the
following we combine the best of both approaches through n-FuRBI: they are
flexible processes that can attain any value for the correlation between the
observables, while at the same time a posterior representation can be derived.
Their construction is based on CRMs and completely random vectors, reviewed in
the next section.

##  3 Some basics on completely random measures

As shown in Lijoi and Prünster (2010), many Bayesian nonparametric models can
be obtained as suitable transformations of CRMs; among others, these include
the Dirichlet process, the Pitman-Yor process and the neutral-to-the-right
priors. The extension of CRMs to the bivariate setting is provided by
completely random vectors $\mu=(\mu_{1},\mu_{2})$, whose components take
values in the space of boundedly finite measures on $\mathbb{X}$ and are such
that, for every collection of pairwise disjoint sets $(A_{i})_{i\geq 1}^{n}$,
the random vectors
$(\mu_{1}(A_{1}),\mu_{2}(A_{1})),\ldots,(\mu_{1}(A_{n}),\mu_{2}(A_{n}))$ are
mutually independent. We focus on the case of no fixed atoms and no
deterministic component, so that the marginal CRMs $\mu_{1}$ and $\mu_{2}$ are
almost surely discrete and can be written as sum of $\mathbb{X}$–valued random
atoms with random weights, i.e.

$\mu_{1}\overset{a.s.}{=}\sum_{i\geq
1}J_{i}\delta_{\kappa_{i}},\quad\mu_{2}\overset{a.s.}{=}\sum_{i\geq
1}W_{i}\delta_{\kappa_{i}}.$

In the following section it will be convenient to use the reparametrization
$\kappa_{i}=(\theta_{i},\phi_{i})\in\mathbb{X}=\mathbb{X}_{1}\times\mathbb{X}_{2}$.
Such completely random vectors are characterized by the Lévy-Khintchine
representation

$\mathbb{E}\left\\{e^{-\mu_{1}(f_{1})-\mu_{2}(f_{2})}\right\\}=\exp\left[-\int\limits_{\mathbb{R}^{2}_{+}\times\mathbb{X}}\\{1-e^{-s_{1}\,f_{1}(x)-s_{2}\,f_{2}(x)}\\}\,v(\mathrm{d}s_{1},\mathrm{d}s_{2},\mathrm{d}x)\right]$
(6)

where $\mu_{i}(f_{i})=\int_{\mathbb{X}}f_{i}(x)\mu_{i}(\mathrm{d}x)$ for
$\mathbb{R}^{+}$-valued $f_{i}$ and
$v(\mathrm{d}s_{1},\mathrm{d}s_{2},\mathrm{d}x)$ is the joint Lévy intensity.
We shall focus on the homogeneous case, in which jumps $(J_{j})_{j\geq 1}$ and
locations $(X_{j})_{j\geq 1}$ are independent. In terms of Lévy intensity it
reads
$v(\mathrm{d}s_{1},\mathrm{d}s_{2},\mathrm{d}x)=\rho(\mathrm{d}s_{1},\mathrm{d}s_{2})\alpha(\mathrm{d}x)$
for some finite measure $\alpha$ on $\mathbb{X}$ and measure $\rho$. Moreover,
in the sequel we will also need the joint and marginal Laplace exponents given
by

$\displaystyle\psi_{b}(\lambda_{1},\lambda_{2}):=$
$\displaystyle\int\limits_{\mathbb{R}^{2}_{+}\times\mathbb{X}}(1-e^{-\lambda_{1}s_{1}-\lambda_{2}s_{2}})\rho(\mathrm{d}s_{1},\mathrm{d}s_{2})\alpha(\mathrm{d}x),\quad\lambda_{1}>0,\lambda_{2}>0.$
$\displaystyle\psi(\lambda):=$
$\displaystyle\int\limits_{\mathbb{R}_{+}\times\mathbb{X}}(1-e^{-\lambda
s})\rho(\mathrm{d}s)\alpha(\mathrm{d}x)\quad\lambda>0,$

For an exhaustive account on CRMs, we refer to Kingman (1967, 1993).
Completely random vectors and CRMs are often normalized to obtain random
probability measures, as introduced in Regazzini et al. (2003), i.e.
$p(\cdot)={\mu(\cdot)}/{\mu(\mathbb{X})}$. Notice that in principle any random
measure $\mu$ such that $\mathbb{P}(0<\mu(\mathbb{X})<\infty)=1$ can be
normalized in order to define a random probability measure. However, the
strength of completely random vectors and measures lies in their
Lévy–Khintchine representations and unique correspondence with the associated
Lévy intensity, which allow a high degree of analytical tractability. CRMs and
the corresponding normalized probabilities have been extensively studied to
model exchangeable data (see, for instance, James et al., 2006, 2009, 2010;
Lijoi and Prünster, 2010; Favaro et al., 2016; Camerlenghi et al., 2018).
Similarly, a completely random vector can be used to model dependence between
two groups. For more details on completely random vectors and an interesting
account of their dependence structure, we refer to Catalano et al. (2021);
Catalano et al. (2023). Since the two measures in the vector share all the
atoms, by virtue of Proposition 2 the induced model yields non–negative
correlation between samples. The issue is addressed in the next section, by
means of a novel class of random probability measures that leverage the
dependence structure specifed for the atoms.

##  4 Full-range borrowing of information nonparametric prior

###  4.1 Definition and first properties

In this section we introduce n-FuRBI and for simplicity we still consider only
the case of two samples with the same a priori marginal distribution.

######  Definition 1.

Consider a completely random vector $(\mu_{1},\mu_{2})$ on $\mathbb{X}^{2}$
with Lévy intensity

$v(\mathrm{d}s_{1},\mathrm{d}s_{2},\mathrm{d}x_{1},\mathrm{d}x_{2})=\rho(\mathrm{d}s_{1},\mathrm{d}s_{2})\>\alpha(\mathrm{d}x_{1},\mathrm{d}x_{2}),$

where $\alpha(\mathrm{d}x_{1},\mathrm{d}x_{2})=\theta
G_{0}(\mathrm{d}x_{1},\mathrm{d}x_{2})$, where
$\theta=\alpha(\mathbb{X}^{2})\in(0,+\infty)$, and $G_{0}$ is a non-atomic
probability measure on $\mathbb{X}^{2}$ such that
$G_{0}(\cdot\times\mathbb{X})=G_{0}(\mathbb{X}\times\cdot)=P_{0}(\cdot)$. Then
$\tilde{\mu}_{1}$ and $\tilde{\mu}_{2}$ defined as

$\tilde{\mu}_{1}(\cdot)=\mu_{1}(\mathbb{X}\times\cdot)\qquad\tilde{\mu}_{2}(\cdot)=\mu_{2}(\cdot\times\mathbb{X})$

are CRMs with Full-Range Borrowing of Information (FuRBI CRMs) and underlying
Lévy intensity $v$. The normalized versions
$\tilde{p}_{j}(\cdot)=\tilde{\mu}_{j}(\cdot)/\tilde{\mu}_{j}(\mathbb{X})$ for
$j=1,2$ are said normalized CRMs with Full-Range Borrowing of Information
(n-FuRBI).

Essentially, first a pair of random measures endowed with the same locations
is constructed on the product space $\mathbb{X}^{2}$; as a second step, the
coordinates of each pair of atoms are split. Thus, the n-FuRBI admit a
representation as in (2) and (3). In general FuRBI CRMs are not completely
random vectors, because the joint sampling of the atoms forbids the
independence of the vector evaluated on pairwise disjoint sets. However, the
representation in terms of a completely random vector in the product space is
useful to characterize the joint law of the FuRBI CRMs, as shown in the
following proposition.

######  Proposition 6.

Let $(\tilde{\mu}_{1},\tilde{\mu}_{2})$ be a vector of FuRBI CRMs. Then

  * (i)

$\tilde{\mu}_{1}$ and $\tilde{\mu}_{2}$ are CRMs with intensity
$\rho(\mathrm{d}s)\theta P_{0}(\mathrm{d}x)$, where
$\rho(\mathrm{d}s)=\int_{\mathbb{R}_{+}}\rho(\mathrm{d}s_{1},\mathrm{d}s)$.

  * (ii)

For any $A$ and $B$, the following equality holds

$\begin{aligned}
\text{E}\bigl{[}\mathrm{e}^{-\lambda_{1}\tilde{\mu}_{1}(A)-\lambda_{2}\tilde{\mu}_{2}(B)}\bigr{]}=\exp\\{-G_{0}(A\times
B^{c})\psi(\lambda_{1})-G_{0}(A^{c}\times B)\psi(\lambda_{2})\\}\\\
\times\>\exp\\{-G_{0}(A\times
B){\psi_{b}}(\lambda_{1},\lambda_{2})\\},\end{aligned}$

where $\psi$ denotes the common marginal Laplace exponent and $\psi_{b}$ the
joint Laplace exponent of $(\mu_{1},\mu_{2})$.

  * (iii)

The joint law of $(\tilde{\mu}_{1},\tilde{\mu}_{2})$ is characterized by the
joint Lévy intensity of $(\mu_{1},\mu_{2})$.

The next proposition shows that the $\beta$ and $\gamma$ associated to any
couple of n-FuRBI can be computed through their Laplace exponents.

######  Proposition 7.

Consider $(\tilde{p}_{1},\tilde{p}_{2})$ n-FuRBI. Then the probability of a
tie and of a hyper-tie are respectively

$\beta=-\int_{\mathbb{R}_{+}}u\left\\{\frac{\mathrm{d}^{2}}{\mathrm{d}u^{2}}\psi(u)\right\\}e^{-\psi(u)}\,\mathrm{d}u,\quad\gamma=-\int_{\mathbb{R}^{2}_{+}}\left\\{\frac{\partial^{2}}{\partial
u_{1}\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}.$

Thus, the crucial value of $\gamma$ can be obtained by computing, analytically
or numerically, a bivariate integral. The two results above show a recurrent
trait of our approach: interesting quantities will be usually rewritten in
terms of the original completely random vector, in order to exploit its
analytical tractability. We conclude this section with two examples of FuRBI
CRMs, that also show how some existing constructions can be obtained as
special cases.

######  Example 2 (FuRBI CRMs with equal jumps).

Let
$\rho(\mathrm{d}s_{1})\delta_{s_{1}}(\mathrm{d}s_{2})\,\theta\,G_{0}(\mathrm{d}x_{1},\mathrm{d}x_{2})$
be the underlying Lévy intensity. The series representation of the
corresponding FuRBI CRMs is

$\tilde{\mu}_{1}(\cdot)\overset{a.s.}{=}\sum_{k\geq
1}W_{k}\delta_{\theta_{k}}\qquad\qquad\tilde{\mu}_{2}(\cdot)\overset{a.s.}{=}\sum_{k\geq
1}W_{k}\delta_{\phi_{k}}\qquad\qquad\text{with}\enskip(\theta_{k},\phi_{k})\overset{i.i.d}{\sim}G_{0}.$

Therefore, $\gamma=\beta$, so that a tie and a hyper–tie are observed with the
same probability.

######  Example 3 (Extended Compound FuRBI CRMs).

Consider the Lévy intensity

$v(\mathrm{d}s_{1},\mathrm{d}s_{2},\mathrm{d}x_{1},\mathrm{d}x_{2})=\int
z^{-2}h(s_{1}/z,\,s_{2}/z)\,\mathrm{d}s_{1}\mathrm{d}s_{2}v^{*}(\mathrm{d}z)\,\theta\,G_{0}(\mathrm{d}x_{1},\mathrm{d}x_{2}),$

where $h$ is some density and $v^{*}$ is a Lévy intensity that satisfies

$\int
z^{-2}\int\min\\{1,||{s}||\\}h(s_{1}/z,\,s_{2}/z)\,\mathrm{d}s_{1}\mathrm{d}s_{2}v^{*}(\mathrm{d}z)<\infty,\quad||{s}||=\sqrt{s_{1}^{2}+s_{2}^{2}}.$

The series representation of the corresponding FuRBI CRMs is

$\tilde{\mu}_{1}(\cdot)\overset{a.s.}{=}\sum_{k\geq
1}m_{1,k}W_{k}\delta_{\theta_{k}}\qquad\qquad\tilde{\mu}_{2}(\cdot)\overset{a.s.}{=}\sum_{k\geq
1}m_{2,k}W_{k}\delta_{\phi_{k}}$

where $(\theta_{k},\phi_{k})\overset{i.i.d}{\sim}G_{0}$ and
$(m_{1,k},m_{2,k})\overset{iid}{\sim}h$. When $G_{0}$ is degenerate on the
main diagonal, one retrieves the class of compound random measures introduced
by Griffin and Leisen (2017).

###  4.2 Correlation structure between n-FuRBI

In order to analyze the dependence between the marginal n-FuRBI priors
$\tilde{p}_{1}$ and $\tilde{p}_{2}$, it is useful to compute the correlation
of the random probability measures evaluated on the same set $A$. In all the
existing CRM-based models such a correlation does not depend on the specific
set considered and, hence, it is often used as a global measure of dependence.
The next proposition provides the covariance structure between two n-FuRBI.

######  Proposition 8.

Let $\tilde{p}_{1}$ and $\tilde{p}_{2}$ be n-FuRBI. Then for any $A,B$, such
that $0\leq P_{0}(A)\leq 1$ and $0\leq P_{0}(B)\leq 1$, we have $\mbox{\rm
cov}(\tilde{p}_{1}(A),\tilde{p}_{2}(B))=\gamma\left[G_{0}(A\times
B)-P_{0}(A)P_{0}(B)\right]$ and

$\mbox{\rm
corr}(\tilde{p}_{1}(A),\tilde{p}_{2}(B))=\frac{\gamma}{\beta}\frac{G_{0}(A\times
B)-P_{0}(A)P_{0}(B)}{\sqrt{P_{0}(A)(1-P_{0}(A))P_{0}(B)(1-P_{0}(B))}}.$

By setting $A=B$, from the previous results one immediately deduces that
$\mbox{\rm cov}(\tilde{p}_{1}(A),\tilde{p}_{2}(A))=\gamma\left[G_{0}(A\times
A)-P_{0}(A)^{2}\right]$ and

$\mbox{\rm
corr}(\tilde{p}_{1}(A),\tilde{p}_{2}(A))=\frac{\gamma}{\beta}\frac{G_{0}(A\times
A)-P_{0}(A)^{2}}{P_{0}(A)(1-P_{0}(A))}.$

Unlike what usually happens with existing models, here the correlation can be
negative, when $A$ is such that $G_{0}(A\times A)<P_{0}(A)^{2}$, that is when
$G_{0}$ exhibits a repulsive behaviour between the coordinates in
$\mathbb{X}^{2}$. Moreover, the correlation depends on the specific set on
which the two measures are evaluated and, therefore, it has to be interpreted
as a local measure of dependence. See Section S3. for an illustration of this
phenomenon on sets of the form $(-\infty,x)$.

######  Example 4 (n-FuRBI with equal jumps).

In this case, Proposition 3 entails $\beta=\gamma$. Therefore

$\mbox{\rm
corr}\left(\tilde{p}_{1}(A),\tilde{p}_{2}(A)\right)=\frac{G_{0}(A\times
A)-P_{0}(A)^{2}}{P_{0}(A)(1-P_{0}(A))}.$

Moreover, still by virtue of Proposition 3, for a given $G_{0}$ this is the
highest possible correlation in absolute value.

Proposition 4 then provides the correlation between the observables, which is
even more important from a modeling perspective.

######  Example 5 (Gamma n-FuRBI with equal jumps).

If the common marginal is the law of a Dirichlet process, then $\mbox{\rm
corr}(X_{i},Y_{j})=\rho_{0}/(1+\theta)$. Choosing appropriately $\rho_{0}$ and
$\theta$ the entire range $(-1,1)$ becomes available.

Note that hyper-ties allow to perform a more general type of borrowing,
compared to ties, even when the correlation is positive. While ties are a
useful construction to model multiple samples that share certain values/latent
parameters, hyper-ties can borrow information even when the two samples have
no common values/latent parameter. This aspect will play a crucial role in the
data-analyses of Sections 6.3 and 6.4; for these the assumption of common
values would be highly unrealistic.

##  5 Inference

###  5.1 Posterior Characterization

Having provided an exhaustive description of the a priori properties of
n-FuRBI, the following key step is to provide a tractable posterior
characterization. Conjugacy is out of question here: even in the exchangeable
context it is a property characterizing the Dirichlet process (see James et
al., 2006). Nevertheless, conditional on a set of suitable latent variables,
the posterior distribution of the original completely random vector
$(\mu_{1},\mu_{2})$ turns out to be again a completely random vector leading
to a neat posterior characterization and viable methods for sampling.

Consider a sample of $n$ observations $(X_{i})_{i=1}^{n}$ from $\tilde{p}_{1}$
with unique values
$\underline{X}_{n}^{*}=\left(X^{*}_{1},\dots,X^{*}_{k}\right)$ and associated
multiplicities $(n_{1},\dots,n_{k})$; analogously, consider $m$ observations
$(Y_{j})_{j=1}^{m}$ from $\tilde{p}_{2}$ with unique values
$\underline{Y}_{m}^{*}=\left(Y^{*}_{1},\dots,Y^{*}_{c}\right)$ and
multiplicities $(m_{1},\dots,m_{c})$. While it is immediate to check for ties,
hyper-ties cannot be identified from the data. To this end, we define a latent
random element ${p}$ encoding the hyper-ties, such that
${p}=\left\\{(i_{l},j_{l})\right\\}_{l}$, where $(i,j)$, with $1\leq i\leq k$
and $1\leq j\leq c$, denotes a hyper-tie between $X_{i}^{*}$ and $Y_{j}^{*}$.
Moreover $(i,0)$, with $1\leq i\leq k$, denotes that $X_{i}^{*}$ does not form
a hyper-tie with any value in $\underline{Y}_{m}^{*}$ and $(0,j)$, with $1\leq
j\leq c$, denotes that $Y_{j}^{*}$ does not form an hyper-tie with any value
in $\underline{X}_{n}^{*}$.

Therefore, if $(i,j)\in{p}$ with $i\neq 0$ and $j\neq 0$, it means that
$X_{i}^{*}$ and $Y_{j}^{*}$ come from the same pair of atoms in representation
(3). Instead, $(i,0)\in{p}$ implies that $X_{i}^{*}$ is the only value
associated to a specific pair, and similarly for $Y_{j}^{*}$ if $(0,j)\in{p}$.
Since we are working with unique values, it is clear that each $X_{i}^{*}$ and
$Y_{j}^{*}$ can form at most one hyper-tie, i.e. it is associated to a unique
member of ${p}$. This justifies the following formal definition.

######  Definition 2.

We say that ${p}=\left\\{(i_{l},j_{l})\right\\}_{l}$ is a compatible hyper-
ties structure for $(X_{i})_{i=1}^{n}$ and $(Y_{j})_{j=1}^{m}$ if, firstly,
for any $1\leq i\leq k$, there exists exactly one $i_{l}$ such that $i_{l}=i$,
thus each element of $\underline{X}_{n}^{*}$ forms at most one hyper-tie;
secondly, for any $1\leq j\leq c$, there exists exactly one $j_{l}$ such that
$j_{l}=j$, thus each element of $\underline{Y}_{m}^{*}$ forms at most one
hyper-tie; lastly, for any $l$, if $i_{l}=0$ then $j_{l}\neq 0$, thus at least
one coordinate refers to an element of $\underline{X}_{n}^{*}$ or
$\underline{Y}_{m}^{*}$.

As a simple example, suppose that $\underline{X}_{n}$ and $\underline{Y}_{m}$
contain respectively $2$ and $1$ unique values. Then $k=2$, $c=1$ and the
support of ${p}$ is

$\mathcal{P}=\biggl{\\{}\\{(1,1),(2,0)\\},\\{(1,0),(2,1)\\},\\{(1,0),(2,0),(0,1)\\}\biggr{\\}}.$

Once the latent structure ${p}$ is identified, its elements can be
conveniently partitioned into the set $\Delta_{{p}}=\left\\{(i,j)\in{p}\mid
i\neq 0\text{ and }j\neq 0\right\\}$, which includes all the hyper-ties, and
the sets $\Delta^{1}_{{p}}=\left\\{(i,j)\in{p}\mid j=0\right\\}$ and
$\Delta^{2}_{{p}}=\left\\{(i,j)\in{p}\mid i=0\right\\}$. If $X_{i}^{*}$ and
$Y_{j}^{*}$ form a hyper-tie, it means that $(X_{i}^{*},Y_{j}^{*})$ is an
actual atom in representation (3). Instead, if $X_{i}^{*}$ does not form a
hyper-tie, we have a partial knowledge of the original pair: the unknown
second coordinate can be sampled from $P_{X_{i}^{*}}(\cdot)$, that is the
conditional distribution given $X_{i}^{*}$, induced by the joint measure
$G_{0}$, which will henceforth be assumed to be non–atomic. A similar argument
applies if $Y_{j}^{*}$ does not form a hyper-tie.

In order to simplify notation, we set $g_{i,j}=g_{0}(X_{i}^{*},Y_{j}^{*})$,
$g_{i,0}=p_{0}(X_{i}^{*})$, and $g_{0,j}=p_{0}(Y_{j}^{*})$, where $g_{0}$ and
$p_{0}$ are the density functions of $G_{0}$ and $P_{0}$ respectively, that we
assume exist with respect to suitable dominating measures. Finally, we
consider the following integrals

$\tau_{n,m}(\underline{u})=\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}^{n}s_{2}^{m}\,\rho(\mathrm{d}s_{1},\mathrm{d}s_{2}),\quad\underline{u}=(u_{1},u_{2}),$

where often $n$ and $m$ will be equal to $n_{i}$ and $m_{j}$, with $1\leq
i\leq k$ and $1\leq j\leq c$. For consistency, we set $n_{0}=m_{0}=0$.

The key result of the section relies on a latent structure that is identified
by random variables whose conditional distributions, given $(X_{i})_{i=1}^{n}$
and $(Y_{j})_{j=1}^{m}$, are available. Indeed, these random variables are
given by ${p}$, whose probability mass function is proportional to

$\left(\prod_{(i,j)\in
p}g_{i,j}\right)\int_{\mathbb{R}_{+}^{2}}u_{1}^{n-1}u_{2}^{m-1}\prod_{(i,j)\in{p}}\tau_{n_{i},m_{j}}(\underline{u})\,\mathrm{e}^{-\psi_{b}(\underline{u})}\,\mathrm{d}\underline{u},$

the vector $(U_{1},U_{2})$, whose density on $\mathbb{R}_{+}^{2}$ is
proportional to
$u_{1}^{n-1}u_{2}^{m-1}\prod_{(i,j)\in{p}}\tau_{n_{i},m_{j}}(\underline{u})e^{-\psi_{b}({u})}$,
the variables $\\{Z_{i}^{x}\\}_{i}$, whose distribution is
$P_{X_{i}^{*}}(\cdot)$, for any $i=1,\dots k$, and $\\{Z_{j}^{y}\\}_{j}$,
whose distribution is $P_{Y_{j}^{*}}(\cdot)$, for any $j=1,\dots,c$. We are
now ready to state the key posterior characterization.

######  Theorem 1.

Let$(X_{i})_{i=1}^{n}$ and $(Y_{j})_{j=1}^{m}$ be from model (13), with $Q$
being the law of a n-FuRBI. Then, the distribution of $(\mu_{1},\mu_{2})$
conditional on $(X_{i})_{i=1}^{n}$, $(Y_{j})_{j=1}^{m}$ and the set of latent
variables $(p,U_{1},U_{2},$ $\\{Z_{i}^{x}\\}_{i},\\{Z_{j}^{y}\\}_{j})$ is

$(\hat{\mu}_{1},\hat{\mu}_{2})+\sum_{(i,j)\in\Delta_{{p}}}{J}_{i,j}\delta_{\left(X_{i}^{*},Y_{j}^{*}\right)}+\sum_{(i,j)\in\Delta^{1}_{{p}}}{J}_{i,0}\delta_{\left(X_{i}^{*},Z_{i}^{x}\right)}+\sum_{(i,j)\in\Delta^{2}_{{p}}}{J}_{0,j}\delta_{\left(Z_{j}^{y},Y_{j}^{*}\right)},$

where $(\hat{\mu}_{1},\hat{\mu}_{2})$ is a completely random vector with
intensity
$e^{-U_{1}s_{1}-U_{2}s_{2}}\rho(\mathrm{d}s_{1},\mathrm{d}s_{2})G_{0}(\mathrm{d}x)$
and ${J}_{i,j}=(J^{1}_{i,j},J^{2}_{i,j})$, with $i=0,\ldots,k$ e
$j=0,\ldots,c$, are jumps with density proportional to

$s_{1}^{n_{i}}s_{2}^{m_{j}}e^{-U_{1}s_{1}-U_{2}s_{2}}\rho(\mathrm{d}s_{1},\mathrm{d}s_{2}).$

Moreover $(\hat{\mu}_{1},\hat{\mu}_{2})$ and ${J}_{i,j}$ are independent.

Conditional on the latent variables, the structure is quite intuitive: the
posterior is the law of a completely random vector with modified intensity and
fixed locations, given by the pairs formed by the hyper-ties. This is somehow
reminiscent of the posterior structures of exchangeable models (James et al.,
2009; Lijoi and Prünster, 2010), with the key novelty played by the new notion
of hyper-ties, in addition to the identification of a suitable latent
structure.

The distribution of the latent variables admits a nice interpretation. For
instance, the mass function of the latent structure ${p}$ is the product of
two terms: the probability of observing the number of hyper-ties identified by
${p}$ times the likelihood that exactly those pairs are formed, through the
density function $g_{0}$. Thus, thanks to the homogeneity of the original
completely random vector, we observe a separate effect for jumps and locations
on this hidden clustering structure. The next corollary shows how the
posterior distribution of the normalized measures can be deduced from Theorem
1. The statement focuses on $p_{1}$, though an analogous representation holds
also for $p_{2}$.

######  Corollary 2.

Under the same assumptions of Theorem 1, conditional on $(X_{i})_{i=1}^{n}$,
$(Y_{j})_{j=1}^{m}$ and the latent variables
$(p,U_{1},U_{2},\\{Z_{i}^{x}\\}_{i},\\{Z_{j}^{y}\\}_{j})$, the random
probability measure $p_{1}$ in (3) equals in distribution

$\displaystyle w_{1}\frac{\hat{\mu}_{1}}{T_{1}}$
$\displaystyle+w_{2}\frac{\sum_{(i,j)\in\Delta_{{p}}}J^{1}_{i,j}\delta_{\left(X_{i}^{*},Y_{j}^{*}\right)}}{\sum_{(i,j)\in\Delta_{{p}}}J^{1}_{i,j}}+w_{3}\frac{\sum_{(i,j)\in\Delta^{1}_{{p}}}J^{1}_{i,0}\delta_{\left(X_{i}^{*},Z_{i}^{x}\right)}}{\sum_{(i,j)\in\Delta^{1}_{{p}}}J^{1}_{i,0}}+w_{4}\frac{\sum_{(i,j)\in\Delta^{2}_{{p}}}J^{1}_{0,j}\delta_{\left(Z_{j}^{y},Y_{j}^{*}\right)}}{\sum_{(i,j)\in\Delta^{2}_{{p}}}J^{1}_{0,j}},$

where $T_{1}=\hat{\mu}_{1}(\mathbb{X}\times\mathbb{X})$, while

$w_{1}\propto T_{1},\quad
w_{2}\propto\sum_{(i,j)\in\Delta_{{p}}}J^{1}_{i,j},\quad
w_{3}\propto\sum_{(i,j)\in\Delta^{1}_{{p}}}J^{1}_{i,0},\quad
w_{4}\propto\sum_{(i,j)\in\Delta^{2}_{{p}}}J^{1}_{0,j},$

with the constraint $\sum_{i=1}^{4}w_{i}=1$.

###  5.2 Predictive structure

Prediction of new observations arises naturally within the Bayesian framework,
since it coincides with the estimate of the distribution under a square loss
function. Moreover, it has the merit of providing intuition on how the model
behaves and learns and it can be used to develop marginal algorithms that
avoid the direct sampling of $\tilde{p}_{1}$ and $\tilde{p}_{2}$, which are
infinite-dimensional objects. In Proposition 5 we saw how to sample the first
pair of observations. The next result tackles the general case.

######  Theorem 2.

Consider samples $(X_{i})_{i=1}^{n}$ and $(Y_{j})_{j=1}^{m}$ from model (13),
with the same setting of Theorem 1. Then there exist probability weights
$\xi_{0}$, $\\{\xi_{i}^{x}\\}$ and $\\{\xi_{j}^{y}\\}$ such that

$\mathbb{P}\bigl{(}X_{n+1}\in
C\mid(X_{i})_{i=1}^{n},\,(Y_{j})_{j=1}^{m}\bigr{)}=\xi_{0}P_{0}(C)+\sum_{i=1}^{k}\xi_{i}^{x}\delta_{X_{i}^{*}}(C)+\sum_{j=1}^{c}\xi_{j}^{y}P_{Y_{j}^{*}}\left(C\right).$

Analogously, there exist probability weights $\eta_{0}$, $\\{\eta_{i}^{x}\\}$
and $\\{\eta_{j}^{y}\\}$ such that for any $C\in\mathcal{X}$

$\mathbb{P}\bigl{(}Y_{m+1}\in
C\mid(X_{i})_{i=1}^{n},\,(Y_{j})_{j=1}^{m}\bigr{)}=\eta_{0}P_{0}(C)+\sum_{j=1}^{c}\eta_{j}^{y}\delta_{Y_{j}^{*}}(C)+\sum_{i=1}^{k}\eta_{i}^{x}P_{X_{i}^{*}}\left(C\right).$

Explicit formulae for the weights are available in the proof of Theorem 2, in
Section $S1$. In specific cases they can be computed in closed form,
conditional to the latent variables: see e.g. example $S1$ in Section $S4$ for
the Inverse Gaussian case with equal jumps.

Hence, the marginal predictive distributions have a quite intuitive form: they
are linear combinations of the centering distribution $P_{0}$, a weighted
version of the empirical distribution and a last term that depends on the
other sample. The crucial differences with respect to prediction rules arising
in the exchangeable case (Lijoi and Prünster, 2010; De Blasi et al., 2015) is
the addition of the last term, which clearly shows how posterior inference
changes when incorporating heterogeneous information and performing borrowing
of information.

######  Example 6 (n-FuRBI with equal atoms).

If the joint distribution $G_{0}$ is degenerate such that the atoms are
completely shared between $\tilde{p}_{1}$ and $\tilde{p}_{2}$, then
$P_{Z}(\cdot)=\delta_{Z}(\cdot)$. Therefore, the last term in Theorem 2
becomes a weighted version of the empirical distribution relative to the other
sample.

Algorithms for posterior inference and prediction are derived in Section S4.

##  6 Numerical Illustrations and Real Data Analyses

###  6.1 Bayesian mixture models

Discrete Bayesian models, as the one specified in (13), are usually not
employed directly on the data, but as a building block in hierarchical mixture
models: in this setting ${X}$ and ${Y}$ are hidden values that describes the
clustering structure within the data. Such models have been introduced by Lo
(1984) for the Dirichlet processes and gained popularity thanks also to the
availability of sampling methods for posterior inference (Escobar and West,
1995; Ishwaran and James, 2001; Neal, 2000). Suppose $\\{f(\cdot\mid
x)\,:\,x\in\mathbb{X}\\}$ is a family of probability density kernels on a
space $\mathbb{W}$. Then the model can be formulated in the context of (13) as

$\begin{aligned} W_{i}\mid&X_{i}\overset{\text{ind}}{\sim}f(\cdot\mid
X_{i})\\\
&X_{i}\mid\tilde{p}_{1}\overset{\text{i.i.d.}}{\sim}\tilde{p}_{1}\end{aligned},\quad\begin{aligned}
V_{j}\mid&Y_{j}\overset{\text{ind}}{\sim}f(\cdot\mid Y_{j})\\\
&Y_{j}\mid\tilde{p}_{2}\overset{\text{i.i.d.}}{\sim}\tilde{p}_{2}\end{aligned},\quad(\tilde{p}_{1},\tilde{p}_{2})\sim\text{n-FuRBI}.$

where $(W_{i})_{i=1}^{n}$ and $(V_{j})_{j=1}^{m}$ are the observable samples
and are assumed to be conditionally independent, given $(X_{i})_{i=1}^{n}$ and
$(Y_{j})_{j=1}^{m}$. Integrating out the latent variables $(X_{i})_{i=1}^{n}$
and $(Y_{j})_{j=1}^{m}$, the data are random draws from suitable countable
mixtures, i.e.

$W_{i}\mid\tilde{p}_{1}\overset{iid}{\sim}\int f(\cdot\mid
x)\,\tilde{p}_{1}(\mathrm{d}x),\quad
V_{j}\mid\tilde{p}_{2}\overset{iid}{\sim}\int f(\cdot\mid
y)\,\tilde{p}_{2}(\mathrm{d}y).$

######  Example 7 (Gaussian mixtures).

We assume $f(\cdot\mid x):=N(\cdot\mid x,\sigma^{2})$, with $\sigma^{2}$
positive known constant, to be the normal density. Thus, the latent parameter
is the mean, i.e. $\mathbb{X}=\mathbb{R}$. In this case
$\text{cov}(X_{i},Y_{j})=\text{cov}(W_{i},V_{j})$, so that the joint behavior
of the latent means is reflected on the observations: this shows the
importance of the correlation structure given by Proposition 4 also for
hierarchical models. Alternatively, the latent parameters could specify both
the mean and the variance, with $\mathbb{X}=\mathbb{R}\times\mathbb{R}_{+}$.

The goal is then to draw samples from the posterior distribution given
$(W_{i})_{i=1}^{n}$ and $(V_{j})_{j=1}^{m}$: however this requires to
integrate out all the possible partitions of the $n+m$ latent variables. As
detailed in Section S4., it is possible to devise a Gibbs sampler for drawing
from the posterior distribution of $(X_{i})_{i=1}^{n}$ and
$(Y_{j})_{j=1}^{m}$.

Once a posterior sample $(X_{i})_{i=1}^{n}$ and $(Y_{j})_{j=1}^{m}$ is
generated, relevant quantities of interest can be approximated by exploiting
the conditional independence of $(W_{i})_{i=1}^{n}$ and $(V_{j})_{j=1}^{m}$,
given the latent variables.

###  6.2 Simulation study for density estimation

We consider a simple application with simulated data, in order to understand
how inference changes when taking into account heterogeneous sources of
information. Assume the following generating mechanism:
$W_{i}\overset{\text{i.i.d.}}{\sim}N(\cdot\mid 10,1)$, for $i=1,\dots,20$, and
$V_{j}\overset{\text{i.i.d.}}{\sim}N(\cdot\mid-10,1)$, for $j=1,\dots,100$.
Supposing only the phenomenon associated to the first sample is of interest,
hierarchical mixtures are considered to make prediction on the unknown density
of $W_{i}$. The kernel considered is the one specified in Example 7, with
known $\sigma^{2}=1$ and latent mean $\mu$. Four different approaches for
modelling dependence between $(W_{i})_{i\geq 1}$ and $(V_{i})_{i\geq 1}$ are
devised: the exchangeable approach, according to which sequences ${W}$ and
${V}$ are supposed to form one exchangeable sequence, inducing the highest
positive correlation between $W_{i}$ and $V_{j}$; the independent approach,
according to which the sample $(V_{i})_{i\geq 1}$ is disregarded entirely,
that is $(W_{i})_{i\geq 1}$ and $(V_{i})_{i\geq 1}$ are treated independently;
the hierarchical approach, where we use a hierarchical Dirichlet process (see
Example 1) that corresponds to a classical borrowing of information; the FuRBI
approach, where the underlying random probability measures $\tilde{p}_{1}$ and
$\tilde{p}_{2}$ are n-FuRBI with equal weights and the distribution on the
atoms is
$G_{0}(\cdot\mid\rho_{0})=N_{2}\left(\cdot\mid\underline{0},1,\rho_{0}\right)$
with $\rho_{0}\sim\mbox{Unif}([-1,1])$, where
$N_{2}(\cdot\mid\underline{m},\sigma_{0}^{2},\rho_{0})$ denotes the bivariate
normal distribution with mean vector $\underline{m}$, common variance
$\sigma_{0}^{2}$ and correlation $\rho_{0}$. It can be proven that under this
specification $\mbox{\rm corr}(W_{i},V_{j})=0$, so that a priori ${W}$ and
${V}$ are marginally uncorrelated. The prior specification is purposely
simple, especially regarding the base measure and the concentration parameter,
in order to single out the effect of the borrowing between the two groups as
much as possible.

For the first two cases and the n-FuRBI, the marginal distribution is given by
a Dirichlet process with $\theta=1$ and $P_{0}(\cdot)=N(\cdot\mid 0,1)$;
instead for the hierarchical process the concentration parameters are fixed in
order to match the expected number of different clusters with the other
methods, for a fair comparison. As highlighted in Example 5, n-FuRBI with
equal jumps lead to the most general setting in terms of achievable
correlation between samples; moreover, choosing the marginal processes to
derive from a Gamma process, we can achieve any value in the interval
$(-1,1)$, tuning appropriately the concentration parameter $\theta$.

Figure 1:  Left: mean posterior densities for the case with opposite true
means. Right: mean integrated error (computed on a grid and as the median over
$50$ different samples) for the four estimates, varying the true mean of $V$.

The left panel of Figure 1 shows the performances of the four methods, after
the application of the blocked Gibbs sampler provided in the supporting
material: the mean posterior density (computed pointwise) is depicted. The
exchangeable approach behaves very badly, as expected, because the two samples
clearly have a different distribution. The independent choice leads to a
reasonable estimate, even if it still overestimates the probability mass
around the prior mean (because of the small sample size of the first sample).
The hierarchical estimate is quite good, but our proposal, instead, fits
almost perfectly the target density and seems to exploit the opposite
behaviour of the two phenomena: this is clearly highlighted by the posterior
distribution of $\rho_{0}$, whose approximated mean is close to $-0.9$.

One may wonder whether these superior performances follow from the precise
specification above, with opposite true means. Therefore, we have repeated the
experiment by keeping the same generating mechanism for $W$, but with the true
mean of ${V}$ ranging in the set $\\{-16,-14,\dots,14,16\\}$: the mean
integrated absolute error (computed on a grid and as the median over $50$
different samples) is depicted in the right panel of Figure 1. It is apparent
that the FuRBI approach almost always yields the smallest error, regardless of
the true value. Its performance is close to the exchangeable case only when
the two true means are equal, that is when exchangeability actually holds;
analogously, the n-FuRBI priors yield the highest error when the mean of $V$
corresponds to the prior mean, i.e., when the other group provides less
additional information. The hierarchical process captures the right dependence
when the two means coincide, but can be misled when they are close; finally,
when the second sample is very far from the first one it performs better than
the independent model, probably thanks to the different inner clustering
structure. The results are also summarized in Table 1. Thus, n-FuRBI seem to
be always capable of combining heterogeneous information in the right way; in
particular, at least in this example, they recognize the most useful type of
borrowing of information. In Section S5.1 similar experiments are conducted,
using different data generating distributions: they show that the conclusions
hold even when the data display significantly different features, as
multimodality or heavy tails.

[capbesideposition=right,center,capbesidewidth=6cm]figure[]

Figure 2: Posterior median of the correlation (obtained through $100$
simulation studies) between the three unknown means. Black with triangular
shapes: correlation between the first and second components. Red with square
shapes: correlation between the first and third components. Green with
circular shapes: correlation between the second and third components. Table 1:
Mean integrated absolute error associated to the four methods for some values
of the mean of ${V}$. The values in bold are the smallest ones for each row.

Mean of  V Exch. Ind. FuRBI Hier. -16 1.769 0.995 0.163 0.604 -10 1.769 0.995
0.189 0.592 0 1.737 0.995 0.489 0.587 10 0.205 0.995 0.338 0.397 16 1.666
0.995 0.435 0.592

Finally, we consider a similar application with three groups, in order to see
whether n-FuRBI are able to discern more complex types of dependence. We
assume to observe $W_{1,i}\overset{\text{i.i.d.}}{\sim}N(\cdot\mid 10,1)$,
$W_{2,i}\overset{\text{i.i.d.}}{\sim}N(\cdot\mid-10,1)$, and
$W_{3,i}\overset{\text{i.i.d.}}{\sim}N(\cdot\mid x,1)$, where $i=1,\ldots,20$
and $x\in\\{-10,-9,\dots,10\\}$. Then, for each value of $x$ we apply the same
n-FuRBI with the same weights described above, but where the atoms are
distributed according to

$\displaystyle G_{0}\left(\cdot\right)$
$\displaystyle=N_{3}\left(\cdot\biggl{\lvert}\underline{0},1,\begin{bmatrix}1&\rho_{12}&\rho_{13}\\\
\rho_{12}&1&\rho_{23}\\\ \rho_{13}&\rho_{23}&1\end{bmatrix}\right),$

where $N_{3}(\cdot\mid\mu_{0},\sigma^{2},\Psi)$ denotes a multivariate normal
distribution with mean $\mu_{0}$, all the variances equal to $\sigma^{2}$ and
correlation matrix $\Psi$ and
$\rho_{12},\rho_{13},\rho_{23}\overset{\text{i.i.d.}}{\sim}\mbox{Unif}([-1,1])$.
The posterior medians of $\rho_{12},\rho_{13}$ and $\rho_{23}$ are depicted in
Figure 2, for any value of $x$. The results are in line with our intuition:
the correlation between the first and second component is always close to $-1$
(indeed they have opposite behaviour relative to the prior), while $\rho_{13}$
and $\rho_{23}$ vary linearly with $x$, being positive when the means have the
same sign.

###  6.3 Predicting stocks and bonds returns

Figure 3: Empirical correlation between average stock return and average
commodity return computed on a moving window of 12 months using data from
March 2011 to January 2021.

Findings from the previous section and Section S5.1 suggest that n-FuRBI may
be used to enhance density estimates and prediction in multi-sample data.
Here, the performance is showcased on a real dataset of stock and bond
returns. We collected monthly returns of January 2021 for a sample of 49
stocks portfolios from the Kenneth R. French’s Data Library (data available at
http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html) and
for a sample of 55 commodities from the Primary Commodity Prices Database of
the International Monetary Fund (data available at
https://www.imf.org/en/Research/commodity-prices).

We employ a Bayesian mixture model and assume that stock and bonds returns,
denoted by $W_{i}$ and $V_{j}$, respectively, are sampled from mixtures of
normals where the mixing distributions act on mean and variance of the kernel,
i.e.,

$\displaystyle W_{i}\mid\tilde{p}_{1}\overset{iid}{\sim}\int N(\cdot\mid
x,\sigma^{2}_{w})\,\tilde{p}_{1}(\mathrm{d}x,\mathrm{d}\sigma^{2}_{w})$
$\displaystyle\quad\quad V_{j}\mid\tilde{p}_{2}\overset{iid}{\sim}\int
N(\cdot\mid
y,\sigma^{2}_{v})\,\tilde{p}_{2}(\mathrm{d}y,\mathrm{d}\sigma^{2}_{v}).$

Stocks and commodities exhibit correlation that largely varies over time
ranging from positive to negative values (see, for instance, Bhardwaj and
Dunsby, 2013, and Figure 3). As a consequence, commodities returns contain
useful information to make inference over the distribution of stocks
portfolios, and viceversa. Thus, borrowing of information represents a natural
strategy to improve inference. However, returns may differ even largely in
value between the two sets of financial instruments, especially in periods of
negative correlation. For instance, in our dataset, 53% of the observed stocks
returns are negative, while only 16% of the bonds returns have negative sign.
As such, classical nonparametric borrowing, consisting in sharing of mixture
components, is not appropriate and, as shown in the following, possibly
harmful. We instead make use of n-FuRBI models as prior distribution, i.e.,

$\displaystyle(\tilde{p}_{1},\tilde{p}_{2})\mid\theta,z,G_{0}$
$\displaystyle\sim\text{n-FuRBI}(\theta,\rho,G_{0})$ $\displaystyle\theta$
$\displaystyle\sim\text{Gamma}(\alpha,\beta)$

The base measure $G_{0}$ is chosen so that marginal distributions are given by
normalized CRMs with conjugate Normal-InverseGamma base measure, i.e.

$\displaystyle
G_{0}(\mathrm{d}_{x},\mathrm{d}_{y},\mathrm{d}\sigma^{2}_{w},\mathrm{d}\sigma^{2}_{v}\mid\rho_{0})=$
$\displaystyle\,N_{2}(\mathrm{d}x,\mathrm{d}y\mid{m},\Sigma(\lambda_{1},\lambda_{2},\sigma^{2}_{w},\sigma^{2}_{v}\rho_{0}))$
$\displaystyle\times\text{InvGamma}(\mathrm{d}\sigma^{2}_{w}\mid\alpha_{1},\beta_{1})\times\text{InvGamma}(\mathrm{d}\sigma^{2}_{v}\mid\alpha_{2},\beta_{2})$

with

$m=(m_{1},m_{2})^{\prime}\qquad\text{and}\qquad\Sigma=\begin{bmatrix}\frac{\sigma^{2}_{w}}{\lambda_{1}}&\rho_{0}\,\frac{\sigma_{w}}{\lambda_{1}^{1/2}}\frac{\sigma_{v}}{\lambda_{2}^{1/2}}\\\
\rho_{0}\,\frac{\sigma_{w}}{\lambda_{1}^{1/2}}\frac{\sigma_{v}}{\lambda_{2}^{1/2}}&\frac{\sigma^{2}_{v}}{\lambda_{2}}\end{bmatrix}$

and we use the following joint underlying Lévy intensity
$v(\mathrm{d}s_{1},\mathrm{d}s_{2},\mathrm{d}x_{1},\mathrm{d}x_{2})=\\{z\,[\rho(\mathrm{d}s_{1})\delta_{0}(\mathrm{d}s_{2})+\rho(\mathrm{d}s_{2})\delta_{0}(\mathrm{d}s_{1})]+(1-z)\,\rho(\mathrm{d}s_{1})\delta_{s_{1}}(\mathrm{d}s_{2})\\}\,\theta\,G_{0}(\mathrm{d}x_{1},\mathrm{d}x_{2})$,
with $z\sim\text{Unif}([0,\,1])$. We term the resulting n-FuRBI _additive
n-FuRBI_ , since the series representation of the corresponding FuRBI CRMs is

$\tilde{\mu}_{1}(\cdot)\overset{a.s.}{=}\sum_{k\geq
1}W_{k}\delta_{\theta_{0,k}}+\sum_{k\geq
1}J_{k}\delta_{\theta_{1,k}}\qquad\qquad\tilde{\mu}_{2}(\cdot)\overset{a.s.}{=}\sum_{k\geq
1}W_{k}\delta_{\phi_{0,k}}+\sum_{k\geq 1}V_{k}\delta_{\phi_{2,k}},$

where $(\theta_{0,k},\phi_{0,k})\overset{i.i.d}{\sim}G_{0}$,
$\theta_{1,k}\overset{i.i.d}{\sim}P_{0}$ and
$\phi_{2,k}\overset{i.i.d}{\sim}P_{0}$. When $G_{0}$ is degenerate on the main
diagonal (i.e. $\rho_{0}=1$), one retrieves GM-dependent completely random
measures (Lijoi et al., 2014a, b; Lijoi and Nipoti, 2014). In order to obtain
two Dirichlet processes marginally we set $\rho(s)=s^{-1}e^{-s}$, so that
$\beta=1/(1+\theta)$ and
$\gamma=(1-z)\,{}_{3}F_{2}(\theta-\theta\,z+2,\,1,\,1;\,\theta+2,\,\theta+2;\,1){\theta}/{(1+\theta)^{2}}$,
where ${}_{3}F_{2}$ is the generalized hypergeometric function.

(a) FuRBI with $\rho_{0}\in[-1,1]$

(b) FuRBI with $\rho_{0}=-0.95$

(c) FuRBI with $\rho_{0}=0.95$

(d) Exchangeable model

(e) GM-dependent model

(f) Independent model

Figure 4:  Posterior density estimates for stocks returns.

As for the hyperparameters of the model, we set the a priori expectations
$m_{1}$ and $m_{2}$ in the two groups equal to the empirical averages of the
two groups in December 2020, i.e., the month preceding the data collection,
leading to $m_{1}=5.8591$ and $m_{2}=3.9731$. In the following, we say that a
financial instrument is _outperforming_ if its observed return is higher than
its a priori expected value. In order to assign $\rho_{0}$, we use the results
of Propositions 4 and 5. The elicited $\rho_{0}$ should reflect our prior
opinion about the correlation, which means that it should induce a learning
mechanism agreeing with the following principle: under positive/negative
correlation, conditioning on the event of outperforming commodities, the prior
probability of outperforming/underperforming stocks should increase. Prior
opinion about the correlation can be formulated working with financial experts
and, thanks to n-FuRBI, incorporated through an informative prior on the
parameter $\rho_{0}$. Here, we consider three scenarios: in the first and
second, we derive inferential results under a prior opinion of negative and
positive correlation, respectively, while in the third scenario we assume that
no information on the correlation is available. The three scenarios are
obtained with, respectively, $\rho_{0}=0.95$, $\rho_{0}=-0.95$, and using a
uniform prior on $\rho_{0}$. After standardizing the data, we set the
remaining hyperparameters in a weakly informative way, i.e.
$\lambda_{1}=\lambda_{2}=1$, $\alpha_{1}=\alpha_{2}=2$, and
$\beta_{1}=\beta_{2}=4$. Sensitivity analysis, carried out in Section S6.2,
shows that results are robust with respect to different choices for
$\lambda_{j}$, $\alpha_{j}$ and $\beta_{j}$ for $j=1,2$. We perform $50,000$
iterations of the marginal algorithm (Section S4.1) and discard the first
$10,000$ as burn–in. Section S8. contains results about convergence
diagnostic, mixing performance, and computational times of the algorithm.

Finally, we compare our approach with three alternative models: the
independent model and the exchangeable model, described in the previous
section, and the GM-dependent model from Lijoi et al. (2014b), which performs
classical borrowing based on ties and shares the same additive structure of
additive n-FuRBI.

| ALCPO | MLCPO
---|---|---
FuRBI $\rho_{0}\in[-1,1]$ | -1.2347 | -0.9627
FuRBI $\rho_{0}=-0.95$ | -1.2925 | -1.0115
FuRBI $\rho_{0}=0.95$ | -1.2896 | -1.0149
Exch | -1.5024 | -1.1521
GM-dep | -1.4864 | -1.1557
Ind | -1.3495 | -1.1017
Table 2: ALCPO and MLCPO under the three models. Best performance is
highlighted in bold.

Figure S5 displays the posterior density estimates for stocks returns. The
analogous figure for bonds returns can be found in Section S6.1. Models
employing additive n-FURBI produce density estimates that better resemble the
empirical distribution. The best performance is attained by placing a (non-
informative) prior over the correlation $\rho_{0}$, which leads to a posterior
skewed towards negative values but still quite dispersed (see Figure S7)
reflecting the direction and intensity of the borrowing of information. The
FuRBI models with fixed $\rho_{0}$ perform worse compared to full-borrowing;
nonetheless, thanks to their flexibility, they still produce better results
than other competitors. The GM-dependent and the exchangeable models yield the
worst density estimates in terms of resemblance of the histogram, as expected.
Indeed, the type of borrowing they perform differ from the one allowed by
FuRBIs (even when $\rho_{0}=0.95$), as it is based on ties, which are not
appropriate for the specific problem at hand. Lastly, we note that the
independent model appears to provide a reasonable density estimation, but
presents significantly higher uncertainty.

While Figure S5 provides insight on the model performance, an important caveat
is in order: a too close resemblance of the empirical distribution may
indicate overfitting.

To evaluate the predictive performance, we resort to the conditional
predictive ordinates (CPOs) statistics (see, e.g. Gelfand et al., 1992;
Barrios et al., 2013). Essentially, for each value $i$, we train the model
without the $i$-th observation and compute the predictive density at the
observed point. For the first sample it reads
$\text{CPO}^{w}_{i}=\tilde{f}(w_{i}\mid{w}^{-i},{v})$, for $i=1,\ldots,n$ and
analogously for the second sample we have
$\text{CPO}^{v}_{j}=\tilde{f}(v_{j}\mid{w},{v}^{-j})$, for $j=1,\ldots,m$,
where ${w}$ and ${v}$ denote the vectors of observed returns for,
respectively, stocks and commodities.

Table 2 displays the average logarithmic CPO (ALCPO) and the median
logarithmic CPO (MLCPO) in the overall sample. Higher values correspond to a
better performance, and the n-FuRBI exhibits the best performance.

###  6.4 Clustering of multivariate data with missing entries

Figure 5: Simulated data: left panel shows true clusters locations, right
panel shows complete simulated data for $n=1000$ before applying the
missingness mechanisms.

We now show how to leverage on our methodology to perform borrowing of
information and clustering with multivariate data affected by missing entries.
The n-FuRBI priors are very well suited for this problem: indeed, incomplete
observations can be interpreted as projections of latent complete observations
and, in particular, hyper-ties between incomplete observations can be thought
of as actual ties between complete observations.

We consider a $P$-variate ($P>1$) dataset with missing entries and divide the
dataset into distinct samples based on the missing entries: denote by
$(\underline{W}^{(j_{1},\ldots,j_{l})}_{i},i=1,\ldots,n_{(j_{1},\ldots,j_{l})})$
the sample where $l$ outcomes with labels $(j_{1},\ldots,j_{l})$ are missing.
The dimension of the vector $\underline{W}^{(j_{1},\ldots,j_{l})}_{i}$ is
therefore $P_{j_{1},\ldots,j_{l}}=P-l$. Denote by
$\tilde{q}_{j_{1},\ldots,j_{l}}$ the corresponding unknown distribution, i.e.,

$\underline{W}^{(x)}_{i}\mid\tilde{q}_{x}\overset{iid}{\sim}\tilde{q}_{x}\qquad\text{for
}i=1,\ldots,n_{x}\text{ and }x\in I,$

where $I$ is the index set of all the possible combinations of missing
variables identifying different samples, which are at most $2^{P}-1$.
Independent analyses for each sample should clearly be avoided and classical
nonparametric borrowing cannot even be specified because the support spaces of
different samples differ one from the other.

To perform clustering, we assume that each $\tilde{q}_{x}$ is a mixture of
multivariate normal kernels with diagonal covariance matrix and mixing measure
$\tilde{p}_{x}$ on locations, i.e.

$\underline{W}^{(x)}_{i}\mid\tilde{p}_{x},\underline{\sigma}^{2}\overset{iid}{\sim}\int
N_{P_{x}}(\cdot\mid\underline{\mu}_{x},\underline{\sigma}^{2}_{x})\,\tilde{p}_{x}(\mathrm{d}\underline{\mu}_{x}),$

where
$\underline{\sigma}^{2}=\left(\sigma^{2}_{1},\dots,\sigma^{2}_{P}\right)$,
$\underline{\sigma}^{2}_{x}$ is the restriction of $\underline{\sigma}^{2}$ to
all the elements besides $x$ and
$N_{K}(\cdot\mid\underline{\mu},\underline{\tau}^{2})$ denotes the $K$-variate
normal distribution with mean vector $\underline{\mu}$ and diagonal covariance
matrix given by $\underline{\tau}^{2}$. Independence of the kernel (implied by
the diagonal covariance matrix) is a common assumption in clustering models
for multivariate responses (see, for instance, Gao et al., 2020; Franzolini et
al., 2023): in this way we are forcing the clustering structure to encode all
the dependence across responses. The $\tilde{p}_{x}$ are distributed as

$(\tilde{p}_{x},x\in I)\sim\text{additive n-FuRBI},$

described in Section $6.3$. The atoms of $(\tilde{p}_{x},x\in I)$ are
costrained so that an hyper-tie can be interpreted as an actual tie between
complete observations: moreover the choice of dependent weights allows to
recover group-specific features, if the missingness mechanism is informative.
Section S7.1 provides a discussion of this and contains the details about the
choice of the hyperparameters.

simul | missing | % of missing | n-FuRBI | n-FuRBI | n-FuRBI | mice + | mice +
---|---|---|---|---|---|---|---
number | mechanism | entries | $z=0.2$ | $z=0.5$ | $z=0.8$ | k-means | DPM
n.1 | MCAR | 16.1% | 0.7883 | 0.7882 | 0.7881 | 0.7408 | 0.7734
n.2 | MNAR | 16.7% | 0.7703 | 0.7704 | 0.7706 | 0.6323 | 0.7617
n.3 | MCAR | 35.9% | 0.7292 | 0.7285 | 0.7283 | 0.6786 | 0.7165
n.4 | MNAR | 34% | 0.7304 | 0.7301 | 0.7432 | 0.6391 | 0.7328
Table 3: Rand indexes for $5$ competing methods: 3 n-FuRBI models with varying parameter $z$, mice+k-means and mice+DPM. For n-FuRBI and mice+DPM the posterior expected value is computed averaging over the Rand indexes of all clustering configurations visited by the MCMC chain after burn-in. simul | missing | % of missing | n-FuRBI | n-FuRBI | n-FuRBI | mice + | mice +
---|---|---|---|---|---|---|---
number | mechanism | entries | $z=0.2$ | $z=0.5$ | $z=0.8$ | k-means | DPM
n.1 | MCAR | 16.1% | 4.24 | 4.19 | 4.22 | 3 | 5.48
n.2 | MNAR | 16.7% | 4.59 | 3.29 | 3.37 | 2 | 5.36
n.3 | MCAR | 35.9% | 4.38 | 4.18 | 4.20 | 3 | 7.01
n.4 | MNAR | 34.0% | 4.28 | 4.17 | 4.59 | 2 | 5.85
Table 4:  Estimated number of clusters for $5$ competing methods. The
posterior mean is used for n-FuRBI and mice+DPM, while the number of clusters
is selected by maximizing the average silhouette for mice+k-means. The true
number of clusters is equal to $4$.

First, we conduct a simulation study where data for $n=1,000$ items, $P=3$
responses, and $K=4$ clusters are simulated from a mixture of Gaussian
distributions. Figure 5 shows the locations of the true clusters and the
complete simulated data before deleting entries. Then, different missingness
mechanisms are applied to determine the entries to be treated as missing.
Missing completely at random (MCAR) scenarios are obtained by sampling missing
entries uniformly, while, in missing non at random (MNAR) scenarios the
probability of being missing depends on the true cluster allocation. Different
combinations of missing variables define different samples: the number of
samples ranges from 3 to 6 among simulation scenarios. The detailed
distributions of missing values are provided in Section S7.2. Different values
of the hyperparameter $z$ of the Lévy intensity are considered. Our results
are compared with those obtained with two alternative approaches, called “mice
+ k-means” and “mice + DPM”, which follow a two-steps procedure: first one
imputes missing data by chained equations as implemented in the R package mice
(van Buuren and Groothuis-Oudshoorn, 2011), then, the clustering structure is
estimated with, respectively, k-means and a Dirichlet process mixture. Note
that the number of clusters for k-means is chosen to maximize the average
silhouette. For each run of the n-FuRBI model, we perform $25,000$ iterations
of the MCMC chain and discard the first half as burn-in. Section S8. contains
results about convergence diagnostics, mixing performance, and computational
times of the algorithm. Tables 3 and 5 summarize the performance of the
models. The n-FuRBI priors outperform the alternatives in all scenarios
considered, in term of estimating both the number of clusters and the
clustering configuration, measured by Rand indexes between the estimated
configuration and the true clustering structure. Moreover, the posterior
distribution of n-FuRBI models reflects uncertainty both about the estimated
clustering configuration and about the imputation mechanism, which is instead
ignored by two-step procedures.

Figure 6: Scatter plots of the four scores (after standardization) for the
brandsma dataset. Coordinates of missing data are set equal to their
respective posterior median. Different colors and symbols denote the three
estimated clusters obtained minimizing the variation of information loss with
respect to the posterior distribution.

Finally, we apply the model also on the brandsma dataset (Snijders and Bosker
(2012)), which refers to grade 8 students (age about 11 years) in elementary
schools in the Netherlands (see, Brandsma and Knuver, 1989). The goal is to
cluster $n=4,106$ pupils, based on their IQ verbal score (IQV), IQ performance
score (IQP), language score (LRP), and arithmetic score (APR). The number of
subjects presenting missing entries is $339$ out of $4,106$ (i.e., $8.26\%$).
As before, different combinations of missing variables define different
samples: the number of samples is 7 in the brandsma dataset. In this real data
analysis, the final clustering configuration provides a lower-dimensional
description of the data rather than an estimate of ideal true clusters. Data
are standardized before running the model, so that the sample means and
variances are equal to $0$ and $1$. Figure 6 shows the estimated clustering
configuration obtained minimizing the variation of information loss with
respect to the posterior distribution. The model identifies three clusters,
which show as major tendency that groups of students performing above/below
average for one of the four scores tends to perform above/below average also
for the other scores. In particular, a first cluster includes 53% of the
subjects, which have lower performances: indeed cluster averages of the
standardized scores are IQV$=-0.371$, IQP$=-0.398$, LRP$=-0.387$, and
APR$=-0.435$. Instead, the second cluster, including 44% of the subjects,
retains the best students: the cluster averages of the standardized scores are
IQV$=0.609$, IQP$=0.595$, LRP$=0.629$, and APR$=0.642$. Finally, the students
with the worst scores are allocated to a third cluster whose averages are
IQV$=-2.01$, IQP$=-1.43$, LRP$=-1.90$, and APR$=-1.34$.

##  7 Conclusion

Hyper-ties play a crucial role in driving the Bayesian learning mechanism and
the borrowing of information across samples. However existing nonparametric
priors either do not allow an explicit evaluation of the probability of a
hyper-tie or, when they do, often only non-negative correlation is induced. On
the contrary, n-FuRBIs allow for analytical tractability and may induce either
positive or negative correlation between the random probabilities as well as
across samples resulting in a novel and flexible idea of borrowing of
strength. They are immediately applicable to model multi-sample data through
mixture models, as shown in Section 6.3. Morever, n-FuRBIs also allow for a
variety of interesting extensions, since they can be seen as an effective
building block to model non-trivial dependencies in more complex data
analyses. Future work will further explore these applications.

## Acknowledgements

F. Ascolani, A. Lijoi and I. Prünster are partially supported by MIUR, PRIN
Project P2022H5WZ9. Part of this work was carried out while B. Franzolini was
a Research Fellow at the Agency for Science, Technology and Research, in
Singapore, Republic of Singapore. B. Franzolini is supported by PNRR - PE1
FAIR - CUP B43C22000800006.

## References

  * Arbel and Prünster (2017) Arbel, J. and I. Prünster (2017).  A moment-matching Ferguson & Klass algorithm.  Statistics and Computing 27(1), 3–17.
  * Barrios et al. (2013) Barrios, E., A. Lijoi, L. E. Nieto-Barajas, and I. Prünster (2013).  Modeling with normalized random measure mixture models.  Statistical Science 28(3), 313–334.
  * Bhardwaj and Dunsby (2013) Bhardwaj, G. and A. Dunsby (2013).  The business cycle and the correlation between stocks and commodities.  Journal of Investment Consulting 14(2), 14–25.
  * Brandsma and Knuver (1989) Brandsma, H. and J. Knuver (1989).  Effects of school and classroom characteristics on pupil progress in language and arithmetic.  International Journal of Educational Research 13(7), 777–788.
  * Brillinger (2002) Brillinger, D. R. (2002).  John W. Tukey: his life and professional contributions.  The Annals of Statistics 30(6), 1535–1575.
  * Camerlenghi et al. (2019) Camerlenghi, F., D. B. Dunson, A. Lijoi, I. Prünster, and A. Rodriguez (2019).  Latent nested nonparametric priors.  Bayesian Analysis 14(4), 1303–1356.
  * Camerlenghi et al. (2019) Camerlenghi, F., A. Lijoi, P. Orbanz, and I. Prünster (2019).  Distribution theory for hierarchical processes.  The Annals of Statistics 47(1), 67–92.
  * Camerlenghi et al. (2018) Camerlenghi, F., A. Lijoi, and I. Prünster (2018).  Bayesian nonparametric inference beyond the Gibbs-type framework.  Scandinavian Journal of Statistics 45(4), 1062–1091.
  * Catalano et al. (2023) Catalano, M., H. Lavenant, A. Lijoi, and I. Prünster (2023).  A Wasserstein index of dependence for random measures.  Journal of the American Statistical Association, forthcoming.
  * Catalano et al. (2021) Catalano, M., A. Lijoi, and I. Prünster (2021).  Measuring dependence in the Wasserstein distance for Bayesian nonparametric models.  The Annals of Statistics 49(5), 2916–2947.
  * Cifarelli and Regazzini (1978) Cifarelli, D. M. and E. Regazzini (1978).  Nonparametric statistical problems under partial exchangeability: The role of associative means.  Quaderni Istituto Matematica Finanziaria dell’Università di Torino Serie III 12, 1–36.
  * De Blasi et al. (2015) De Blasi, P., S. Favaro, A. Lijoi, R. H. Mena, I. Prünster, and M. Ruggiero (2015).  Are Gibbs-type priors the most natural generalization of the Dirichlet process?  IEEE Transactions on Pattern Analysis & Machine Intelligence 37(2), 212–229.
  * de Finetti (1938) de Finetti, B. (1938).  Sur la condition d’equivalence partielle.  Actualités Scientifiques et Industrielles 739, 5–18, Translated In: Studies in Inductive and Probability, II. Jeffrey, R. (ed.) University of California Press: Berkeley 1980.
  * Dunson and Park (2008) Dunson, D. and J. Park (2008).  Kernel stick-breaking processes.  Biometrika 95(2), 307–323.
  * Efron and Morris (1977) Efron, B. and C. Morris (1977).  Stein’s paradox in statistics.  Scientific American 236(5), 119–127.
  * Epifani and Lijoi (2010) Epifani, I. and A. Lijoi (2010).  Nonparametric priors for vectors of survival functions.  Statistica Sinica 20(4), 1455–1484.
  * Escobar and West (1995) Escobar, M. D. and M. West (1995).  Bayesian density estimation and inference using mixtures.  Journal of the American Statistical Association 90(430), 577–588.
  * Favaro et al. (2016) Favaro, S., A. Lijoi, C. Nava, B. Nipoti, I. Pruenster, and Y. W. Teh (2016).  On the stick-breaking representation for homogeneous NRMIs.  Bayesian Analysis 11(3), 697–724.
  * Ferguson and Klass (1972) Ferguson, T. S. and M. J. Klass (1972).  A representation of independent increment processes without gaussian components.  The Annals of Mathematical Statistics 43(5), 1634–1643.
  * Foti and Williamson (2013) Foti, N. J. and S. A. Williamson (2013).  A survey of non-exchangeable priors for Bayesian nonparametric models.  IEEE Transactions on Pattern Analysis & Machine Intelligence 37(2), 359–371.
  * Franzolini et al. (2023) Franzolini, B., A. Cremaschi, W. v. d. Boom, and M. De Iorio (2023).  Bayesian clustering of multiple zero-inflated outcomes.  Philosophical Transactions of the Royal Society A 381, 1–16.
  * Gao et al. (2020) Gao, L. L., J. Bien, and D. Witten (2020).  Are clusterings of multiple data views independent?  Biostatistics 21(4), 692–708.
  * Gelfand et al. (1992) Gelfand, A. E., D. K. Dey, and H. Chang (1992).  Model determination using predictive distributions with implementation via sampling-based methods.  Technical Report, Department of Statistics, Stanford University.
  * Gong et al. (2021) Gong, M., P. Liu, F. C. Sciurba, P. Stojanov, D. Tao, G. C. Tseng, K. Zhang, and K. Batmanghelich (2021).  Unpaired data empowers association tests.  Bioinformatics 37(6), 785–792.
  * Griffin et al. (2013) Griffin, J. E., M. Kolossiatis, and M. F. Steel (2013).  Comparing distributions by using dependent normalized random-measure mixtures.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75(3), 499–529.
  * Griffin and Leisen (2017) Griffin, J. E. and F. Leisen (2017).  Compound random measures and their use in Bayesian non-parametrics.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79, 525–545.
  * Ishwaran and James (2001) Ishwaran, H. and L. F. James (2001).  Gibbs sampling methods for stick-breaking priors.  Journal of the American Statistical Association 96(453), 161–173.
  * James et al. (2006) James, L. F., A. Lijoi, and I. Prünster (2006).  Conjugacy as a distinctive feature of the Dirichlet process.  Scandinavian Journal of Statistics 33(1), 105–120.
  * James et al. (2009) James, L. F., A. Lijoi, and I. Prünster (2009).  Posterior analysis for normalized random measures with independent increments.  Scandinavian Journal of Statistics 36(1), 76–97.
  * James et al. (2010) James, L. F., A. Lijoi, and I. Prünster (2010).  On the posterior distribution of classes of random means.  Bernoulli 16(1), 155–180.
  * Kingman (1967) Kingman, J. (1967).  Completely random measures.  Pacific Journal of Mathematics 21(1), 59–78.
  * Kingman (1993) Kingman, J. (1993).  Poisson Processes.  Clarendon Press, Oxford.
  * Lee et al. (2020) Lee, A. M., B.-E. Sæther, and S. Engen (2020).  Spatial covariation of competing species in a fluctuating environment.  Ecology 101(1), e02901.
  * Lijoi et al. (2005) Lijoi, A., R. H. Mena, and I. Prünster (2005).  Hierarchical mixture modeling with normalized inverse-Gaussian priors.  Journal of the American Statistical Association 100(472), 1278–1291.
  * Lijoi and Nipoti (2014) Lijoi, A. and B. Nipoti (2014).  A class of hazard rate mixtures for combining survival data from different experiments.  Journal of the American Statistical Association 109(506), 802–814.
  * Lijoi et al. (2014a) Lijoi, A., B. Nipoti, and I. Prünster (2014a).  Bayesian inference with dependent normalized completely random measures.  Bernoulli 20(3), 1260–1291.
  * Lijoi et al. (2014b) Lijoi, A., B. Nipoti, and I. Prünster (2014b).  Dependent mixture models: clustering and borrowing information.  Computational Statistics & Data Analysis 71, 417–433.
  * Lijoi and Prünster (2010) Lijoi, A. and I. Prünster (2010).  Models beyond the Dirichlet process.  In Bayesian nonparametrics (Hjort, N.L., Holmes, C.C., Müller, P., Walker, S.G. Eds.), pp. 80–136. Cambridge University Press, Cambridge.
  * Lo (1984) Lo, A. Y. (1984).  On a class of Bayesian nonparametric estimates: I. Density estimates.  The Annals of Statistics 12(1), 351–357.
  * Lorenz et al. (2018) Lorenz, D. J., S. Levy, and S. Datta (2018).  Inferring marginal association with paired and unpaired clustered data.  Statistical methods in medical research 27(6), 1806–1817.
  * MacEachern (1999) MacEachern, S. N. (1999).  Dependent nonparametric processes.  In ASA Proceedings of the Section on Bayesian Statistical Science., Alexandria, VA: American Statistical Association.
  * MacEachern (2000) MacEachern, S. N. (2000).  Dependent Dirichlet processes.  Technical Report, The Ohio State University.
  * Müller et al. (2004) Müller, P., F. Quintana, and G. Rosner (2004).  A method for combining inference across related nonparametric Bayesian models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 66(3), 735–749.
  * Müller et al. (2015) Müller, P., F. A. Quintana, A. Jara, and T. Hanson (2015).  Bayesian nonparametric data analysis.  Springer.
  * Neal (2000) Neal, R. M. (2000).  Markov chain sampling methods for Dirichlet process mixture models.  Journal of Computational and Graphical Statistics 9(2), 249–265.
  * Papaspiliopoulos and Roberts (2008) Papaspiliopoulos, O. and G. O. Roberts (2008).  Retrospective Markov chain Monte Carlo methods for Dirichlet process hierarchical models.  Biometrika 95(1), 169–186.
  * Petralia et al. (2012) Petralia, F., V. Rao, and D. B. Dunson (2012).  Repulsive mixtures.  In Advances in Neural Information Processing Systems - NIPS.
  * Quinlan et al. (2017) Quinlan, J. J., F. A. Quintana, and G. L. Page (2017).  Parsimonious hierarchical modeling using repulsive distributions.  arXiv preprint arXiv:1701.04457.
  * Quintana et al. (2022) Quintana, F., P. Müller, A. Jara, and S. MacEachern (2022).  The dependent Dirichlet process and related models.  Statistical Science 37(1), 24–41.
  * Regazzini et al. (2003) Regazzini, E., A. Lijoi, and I. Prünster (2003).  Distributional results for means of normalized random measures with independent increments.  The Annals of Statistics 31(2), 560–585.
  * Rigon and Durante (2021) Rigon, T. and D. Durante (2021).  Logit stick-breaking priors for bayesian density regression.  Journal of Statistical Planning and inference 211, 131–142.
  * Riva-Palacio and Leisen (2021) Riva-Palacio, A. and F. Leisen (2021).  Compound vectors of subordinators and their associated positive Lévy copulas.  Journal of Multivariate Analysis 183, 104728.
  * Rodriguez and Dunson (2011) Rodriguez, A. and D. Dunson (2011).  Nonparametric Bayesian models through probit stick-breaking processes.  Bayesian Analysis 6(1), 145–178.
  * Rodriguez et al. (2008) Rodriguez, A., D. B. Dunson, and A. E. Gelfand (2008).  The nested Dirichlet process.  Journal of the American Statistical Association 103(483), 1131–1154.
  * Sethuraman (1994) Sethuraman, J. (1994).  A constructive definition of Dirichlet priors.  Statistica sinica 4(2), 639–650.
  * Snijders and Bosker (2012) Snijders, T. and R. Bosker (2012).  Multilevel analysis.  Netherlands: SAGE Publications.
  * Teh et al. (2006) Teh, Y. W., M. I. Jordan, M. J. Beal, and D. M. Blei (2006).  Hierarchical Dirichlet processes.  Journal of the American Statistical Association 101(476), 1566–1581.
  * van Buuren and Groothuis-Oudshoorn (2011) van Buuren, S. and K. Groothuis-Oudshoorn (2011).  mice: Multivariate imputation by chained equations in R.  Journal of Statistical Software 45(3), 1–67.
  * Walker (2007) Walker, S. G. (2007).  Sampling the Dirichlet mixture model with slices.  Comm. Statist. Simulation Comput. 36(1-3), 45–54.

Appendix

##  S1 Proofs

###  S1.1 Proofs of Section 1

###### Proof of Proposition 1.

Consider two partially exchangeable sequences $X$ and $Y$ whose elements take
value in $\mathbb{R}$. By de Finetti’s representation theorem, there exist two
random probability measures $\tilde{p}_{1}$ and $\tilde{p}_{2}$ such that

$\left(X_{i},Y_{j}\right)\mid\tilde{p}_{1},\tilde{p}_{2}\overset{\text{iid}}{\sim}\tilde{p}_{1}\times\tilde{p}_{2}.$

Note that
$\mbox{cov}(X_{i},Y_{j})=\mathbb{E}\\{\mbox{cov}(X_{i},Y_{j}\mid\tilde{p}_{1},\tilde{p}_{2})\\}+\mbox{cov}\\{\mathbb{E}(X_{i}\mid\tilde{p}_{1}),\mathbb{E}(Y_{j}\mid\tilde{p}_{2})\\}$,
where the first term equals 0, so that

$\mbox{cov}(X_{i},Y_{j})=\mbox{cov}\left(\int x\,\tilde{p}_{1}(dx),\int
x\,\tilde{p}_{2}(dx)\right),$

and analogously

$\mbox{cov}(X_{i},X_{i^{\prime}})=\mbox{cov}\left(\int
x\,\tilde{p}_{1}(dx),\int x\,\tilde{p}_{1}(dx)\right)=\mbox{var}\left(\int
x\,\tilde{p}_{1}(dx)\right).$

Lastly assume that $\tilde{p}_{1}\overset{d}{=}\tilde{p}_{2}$, where
$\overset{d}{=}$ indicates equality in distribution. By the Cauchy-Schwartz
inequality

$-\mbox{var}\left(\int x\,\tilde{p}_{1}(dx)\right)\leq\mbox{cov}\left(\int
x\,\tilde{p}_{1}(dx),\int x\,\tilde{p}_{2}(dx)\right)\leq\mbox{var}\left(\int
x\,\tilde{p}_{1}(dx)\right),$

which, in terms of the observables, can be equivalently rewritten as

$-\mbox{cov}(X_{i},X_{i^{\prime}})\leq\mbox{cov}(X_{i},Y_{j})\leq\mbox{cov}(X_{i},X_{i^{\prime}}).$

∎

###### Proof of Proposition 2.

By definition of covariance we have

$\text{cov}(X_{i},Y_{j})=\text{cov}\left(\sum_{j\geq
1}J_{j}\theta_{j},\sum_{k\geq 1}W_{k}\phi_{k}\right)=\sum_{j\geq 1}\sum_{k\geq
1}\text{cov}\left(J_{j}\theta_{j},W_{k}\phi_{k}\right).$

For arbitrary $j$ and $k$ we have

$E\left(J_{j}W_{k}\theta_{j}\phi_{k}\right)=E(J_{j}W_{k})E(\theta_{j}\phi_{k})\geq
E(J_{j}W_{k})E(\theta_{j})E(\phi_{k}),$

since cov$(\theta_{j},\phi_{k})\geq 0$. Denoting
$c=E(\theta_{j})=E(\phi_{k})$, we get

$\text{cov}\left(J_{j}\theta_{j},W_{k}\phi_{k}\right)\geq
c^{2}\text{cov}(J_{j},W_{k}).$

Finally, since $\tilde{p}_{1}$ and $\tilde{p}_{2}$ are random probability
measures it holds

$\text{cov}(X_{i},Y_{j})\geq c^{2}\text{cov}\left(\sum_{j\geq
1}J_{j},\sum_{k\geq 1}W_{k}\right)=0,$

which completes the proof. ∎

###  S1.2 Proofs of Section 2

###### Proof of Proposition 3.

Recall that

$\beta:=\sum_{k\geq 1}\text{E}(\bar{J}^{2}_{k})=\sum_{k\geq
1}\text{E}(\bar{W}^{2}_{k})\qquad\gamma:=\sum_{k\geq
1}\text{E}(\bar{J}_{k}\bar{W}_{k}).$

Since

$\text{E}(\bar{J}_{k}\bar{W}_{k})\leq\sqrt{\text{E}(\bar{J}^{2}_{k})\text{E}(\bar{W}^{2}_{k})}=\text{E}(\bar{J}^{2}_{k})$

it follows that $\gamma\leq\beta$. Moreover, the equality holds if and only if
$\bar{J}_{k}\overset{a.s}{=}a_{k}+\bar{W}_{k}$, for any $k$, with
$a_{k}\in\mathbb{R}$. However the equality of marginal distributions implies
$a_{k}=0$. ∎

###### of Proposition 4.

Recall that

$\text{cov}(X_{i},Y_{j})=\text{cov}\left(\sum_{k\geq
1}\bar{J}_{k}\theta_{k},\sum_{h\geq 1}\bar{W}_{h}\phi_{h}\right)=\sum_{k\geq
1}\sum_{h\geq
1}\text{cov}\left(\bar{J}_{k}\theta_{k},\bar{W}_{h}\phi_{h}\right).$

and for arbitrary $k$ and $h$, we have

$\displaystyle\mathbb{E}(\bar{J}_{k}\bar{W}_{h}\theta_{k}\phi_{h})=$
$\displaystyle\text{E}(\bar{J}_{k}\bar{W}_{h})\text{E}(\theta_{k}\phi_{h})$
$\displaystyle=$
$\displaystyle\text{E}(\bar{J}_{k}\bar{W}_{h})\left[\text{E}(\theta_{k}\phi_{k})\mathbbm{1}_{\\{k=h\\}}+\text{E}(\theta_{k})\text{E}(\phi_{h})\mathbbm{1}_{\\{k\neq
h\\}}\right],$

while

$\mathbb{E}(\bar{J}_{k}\theta_{k})=\text{E}(\bar{J}_{k})\text{E}(\theta_{k}).$

Thus, setting $c=E(\theta_{k})=E(\phi_{h})$, we have

$\text{cov}(X_{i},Y_{j})=\sum_{k\geq
1}\text{E}(\bar{J}_{k}\bar{W}_{h})\text{E}(\theta_{k}\phi_{k})-c^{2}\sum_{k\geq
1}\text{E}(\bar{J}_{k})\text{E}(\bar{W}_{k})++c^{2}\sum_{k\geq 1}\sum_{h\neq
k}\text{cov}\left(\bar{J}_{k},\bar{W}_{h}\right)$

where

$\displaystyle\sum_{k\geq 1}\sum_{h\neq
k}\text{cov}\left(\bar{J}_{k},\bar{W}_{h}\right)=$
$\displaystyle\text{cov}\left(\sum_{k\geq 1}\bar{J}_{k},\sum_{h\geq
1}\bar{W}_{h}\right)-\sum_{k\geq
1}\text{cov}\left(\bar{J}_{k}\bar{W}_{k}\right)$ $\displaystyle=$
$\displaystyle-\sum_{k\geq 1}\text{E}(\bar{J}_{k}\bar{W}_{h})+\sum_{k\geq
1}\text{E}(\bar{J}_{k})\text{E}(\bar{W}_{k})$

Putting everything together we obtain

$\displaystyle\text{cov}(X_{i},Y_{j})=\sum_{k\geq
1}\text{E}[\bar{J}_{k}\bar{W}_{k}]\text{cov}(\theta_{k},\phi_{k}).$

Moreover

$\mbox{var}(X_{i})=\mbox{var}(Y_{j})=\int\int
x\,G_{0}(\mathrm{d}x,\mathrm{d}y)=\text{\mbox{var}}(\theta_{k})$

Thus, $\mbox{corr}(X_{i},Y_{j})=\gamma\,\rho_{0}$ proving the second statement
in Proposition 4. Finally, applying the same procedure marginally, we get

$\text{cov}(X_{i},X_{i}^{\prime})=\sum_{k\geq
1}\text{E}(\bar{J}_{k}^{2})\,\text{\mbox{var}}(\theta_{k})$

which proves the first statement in Proposition 4. ∎

###### Proof of Corollary 1.

The result immediately follows from Propositions 3 and 4. ∎

###### Proof of Proposition 5.

Let $\beta$ be the probability of a tie. By definition we get

$\displaystyle\mathbb{P}\left(X_{1}\in A,X_{2}\in B\right)=$
$\displaystyle\,\mathbb{P}(X_{1}\in A,X_{2}\in B\mid X_{1}=X_{2})\beta+$
$\displaystyle+\mathbb{P}(X_{1}\in A,X_{2}\in B\mid X_{1}\neq
X_{2})(1-\beta),$

which, by independence of the atoms, equals

$\displaystyle\mathbb{P}\left(X_{1}\in A,X_{2}\in B\right)=$
$\displaystyle\,\mathbb{P}(X_{1}\in A\in B)\beta+$
$\displaystyle+\mathbb{P}(X_{1}\in A)\mathbb{P}(X_{2}\in B)(1-\beta).$

Analogously, we have

$\displaystyle\mathbb{P}\left(X_{1}\in A,Y_{1}\in B\right)=$
$\displaystyle\,\mathbb{P}(X_{1}\in A,Y_{1}\in B\mid X_{1}\text{ and
}Y_{1}\text{ form an hyper-tie })\gamma+$ $\displaystyle+\mathbb{P}(X_{1}\in
A,Y_{1}\in B\mid X_{1}\text{ and }Y_{1}\text{ do not form an hyper-tie
})(1-\gamma),$

where $\gamma$ is the probability of a hyper-tie, which equals

$\displaystyle\mathbb{P}\left(X_{1}\in A,Y_{1}\in B\right)=$
$\displaystyle\,\mathbb{P}((X_{1},Y_{1})\in A\times B\mid X_{1}\text{ and
}Y_{1}\text{ form an hyper-tie })\gamma+$ $\displaystyle+\mathbb{P}(X_{1}\in
A)\mathbb{P}(Y_{1}\in B)(1-\gamma).$

∎

###  S1.3 Proofs of Section 4

###### Proof of Proposition 6.

The first point follows from the Lévy-Khintchine representation of the Laplace
functional of a CRV. As for (ii), one has

$\displaystyle\text{E}\left(\text{exp}\\{-\lambda_{1}\tilde{\mu}_{1}(A)-\lambda_{2}\tilde{\mu}_{2}(B)\\}\right)=$
$\displaystyle\,\text{E}\left(\text{exp}\\{-\lambda_{1}\mu_{1}(A\times\mathbb{X})-\lambda_{2}\mu_{2}(\mathbb{X}\times
B)\\}\right)$ $\displaystyle=$
$\displaystyle\,\text{E}\bigl{(}\text{exp}\\{-\lambda_{1}\mu_{1}(A\times
B^{c})-\lambda_{1}\mu_{1}(A\times B)+$
$\displaystyle\,-\lambda_{2}\mu_{2}(A^{c}\times B)-\lambda_{2}\mu_{2}(A\times
B)\\}\bigr{)}.$

By independence of evaluations on disjoint sets, $\mu_{1}(C)$ and $\mu_{2}(D)$
are independent if $C\cap D=\emptyset$, so that the right hand side reads

$\displaystyle\text{E}\bigl{(}\text{exp}\\{-\lambda_{1}\tilde{\mu}_{1}(A)-\lambda_{2}\tilde{\mu}_{2}(B)\\}\bigr{)}=$
$\displaystyle\text{E}\left(\text{exp}\\{-\lambda_{1}\mu_{1}(A\times
B^{c})\\}\right)\text{E}\left(\text{exp}\\{-\lambda_{2}\mu_{2}(A^{c}\times
B)\\}\right)\times$
$\displaystyle\times\text{E}\left(\text{exp}\\{-\lambda_{1}\mu_{1}(A\times
B)-\lambda_{2}\mu_{2}(A\times B)\\}\right).$

The result follows upon upon using the expressions of the marginal and joint
Laplace exponents of $\mu_{1}$ and $\mu_{2}$. Since from the joint Lévy
intensity it is possible to recover the joint Laplace exponent, (iii) is also
proved. ∎

In order to prove Proposition 7, we show that

$\mathbb{P}\left(X\in A,Y\in
B\right)=P_{0}(A)P_{0}(B)\left(1-\delta\right)+G_{0}(A\times B)\delta,$

where

$\delta:=-\int_{\mathbb{R}_{+}^{2}}\left\\{\frac{\partial^{2}}{\partial
u_{1}\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}.$

is the probability of a pseudo-tie. We start with three Lemmas.

######  Lemma 1.

If $\psi_{b}$ is the joint Laplace exponent of a CRV, then

$\displaystyle\int_{\mathbb{R}_{+}^{2}}$
$\displaystyle\left\\{\frac{\partial}{\partial
u_{1}}\psi_{b}(u_{1},u_{2})\right\\}\left\\{\frac{\partial}{\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}=1-\delta.$

###### Proof of Lemma S2.1.

Integrating by parts

$\displaystyle\int_{0}^{\infty}$
$\displaystyle\left\\{\frac{\partial}{\partial
u_{1}}\psi_{b}(u_{1},u_{2})\right\\}\left\\{\frac{\partial}{\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}$
$\displaystyle=-\int_{0}^{\infty}\left\\{\frac{\partial}{\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}\left\\{\frac{\partial}{\partial
u_{1}}e^{-\psi_{b}(u_{1},u_{2})}\right\\}\,\mathrm{d}u_{1}$
$\displaystyle=\biggl{[}\left[-\left\\{\frac{\partial}{\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\right]_{0}^{\infty}+\int_{0}^{\infty}\left\\{\frac{\partial^{2}}{\partial
u_{1}\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\biggr{]}$
$\displaystyle=\biggl{[}\left\\{\frac{\partial}{\partial
u_{2}}\psi_{b}(0,u_{2})\right\\}e^{-\psi_{b}(0,u_{2})}+\int_{0}^{\infty}\left\\{\frac{\partial^{2}}{\partial
u_{1}\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\biggr{]}.$

Note that
$\int_{0}^{\infty}\left\\{\frac{\mathrm{d}}{\mathrm{d}u_{2}}\psi_{b}(0,u_{2})\right\\}e^{-\psi_{b}(0,u_{2})}\,\mathrm{d}u_{2}=1$,
by the fundamental theorem of calculus. Thus the result follows immediately. ∎

######  Lemma 2.

We have

$\displaystyle\int_{\mathbb{R}_{+}^{2}}$
$\displaystyle\text{E}\left(e^{-u_{1}\mu_{1}(\mathbb{X}\times\mathbb{X})-u_{2}\mu_{2}(\mathbb{X}\times\mathbb{X})}\mu_{1}(C)\mu_{2}(C)\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}=G_{0}(C)^{2}\left(1-\delta\right)+G_{0}(C)\delta.$

###### Proof of Lemma S2.2.

By independence of evaluations on disjoint sets it follows that

$\displaystyle\int_{\mathbb{R}_{+}^{2}}\text{E}\left(e^{-u_{1}\mu_{1}(\mathbb{X}\times\mathbb{X})-u_{2}\mu_{2}(\mathbb{X}\times\mathbb{X})}\mu_{1}(C)\mu_{2}(C)\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\text{E}\left(e^{-u_{1}\mu_{1}(C)-u_{2}\mu_{2}(C)-u_{1}\mu_{1}(C^{c})-u_{2}\mu_{2}(C^{c})}\mu_{1}(C)\mu_{2}(C)\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\text{E}\left(e^{-u_{1}\mu_{1}(C)-u_{2}\mu_{2}(C)}\mu_{1}(C)\mu_{2}(C)\right\\}\text{E}\left(e^{-u_{1}\mu_{1}(C^{c})-u_{2}\mu_{2}(C^{c})}\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\text{E}\left(\frac{\partial}{\partial
u_{1}}\frac{\partial}{\partial
u_{2}}e^{-u_{1}\mu_{1}(C)-u_{2}\mu_{2}(C)}\right)\text{E}\left(e^{-u_{1}\mu_{1}(C^{c})-u_{2}\mu_{2}(C^{c})}\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\frac{\partial}{\partial
u_{1}}\frac{\partial}{\partial
u_{2}}\left[\text{E}\left(e^{-u_{1}\mu_{1}(C)-u_{2}\mu_{2}(C)}\right)\right]\text{E}\left(e^{-u_{1}\mu_{1}(C^{c})-u_{2}\mu_{2}(C^{c})}\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\frac{\partial}{\partial
u_{1}}\frac{\partial}{\partial
u_{2}}\left\\{e^{-G_{0}(C)\psi_{b}(u_{1},u_{2})}\right\\}e^{-G_{0}(C^{c})\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\frac{\partial}{\partial
u_{1}}\left\\{-G_{0}(C)\frac{\partial}{\partial
u_{2}}\psi_{b}(u_{1},u_{2})e^{-G_{0}(C)\psi_{b}(u_{1},u_{2})}\right\\}e^{-G_{0}(C^{c})\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}.$

Performing the derivative with respect to $u_{1}$, the latter expression can
be written as follows

$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\left\\{G_{0}(C)^{2}\frac{\partial}{\partial
u_{1}}\psi_{b}(u_{1},u_{2})\frac{\partial}{\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-G_{0}(C)\psi_{b}(u_{1},u_{2})}e^{-G_{0}(C^{c})\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}+$
(7)
$\displaystyle+\int_{\mathbb{R}_{+}^{2}}\left\\{-G_{0}(C)\frac{\partial}{\partial
u_{1}\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-G_{0}(C)\psi_{b}(u_{1},u_{2})}e^{-G_{0}(C^{c})\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\left\\{G_{0}(C)^{2}\frac{\partial}{\partial
u_{1}}\psi_{b}(u_{1},u_{2})\frac{\partial}{\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}+$
$\displaystyle+\int_{\mathbb{R}_{+}^{2}}\left\\{-G_{0}(C)\frac{\partial}{\partial
u_{1}\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}$

By Lemma 1 we then obtain

$\int_{\mathbb{R}_{+}^{2}}\text{E}\left(e^{-u_{1}\mu_{1}(\mathbb{X}\times\mathbb{X})-u_{2}\mu_{2}(\mathbb{X}\times\mathbb{X})}\mu_{1}(C)\mu_{2}(C)\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}=G_{0}(C)^{2}\left(1-\delta\right)+G_{0}(C)\delta,$

as desired. ∎

######  Lemma 3.

Let $C,D$ be such that $C\cap D=\emptyset$. Then

$\int_{\mathbb{R}_{+}^{2}}\text{E}\left(e^{-u_{1}\mu_{1}(\mathbb{X}\times\mathbb{X})-u_{2}\mu_{2}(\mathbb{X}\times\mathbb{X})}\mu_{1}(C)\mu_{2}(D)\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}=G_{0}(C)G_{0}(D)\left(1-\delta\right)$

###### Proof of Lemma S3.

Let $Y=\left(C\cup D\right)^{c}$. Since $C$ and $D$ are disjoint, by
independence of evaluations on disjoint sets it holds

$\displaystyle\int_{\mathbb{R}_{+}^{2}}$
$\displaystyle\text{E}\left(e^{-u_{1}\mu_{1}(\mathbb{X}\times\mathbb{X})-u_{2}\mu_{2}(\mathbb{X}\times\mathbb{X})}\mu_{1}(C)\mu_{2}(D)\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}$
$\displaystyle\text{E}\left(e^{-u_{1}\mu_{1}(C\cup D)-u_{2}\mu_{2}(C\cup
D)}\mu_{1}(C)\mu_{2}(D)\right\\}\text{E}\left\\{e^{-u_{1}\mu_{1}(Y)-u_{2}\mu_{2}(Y)}\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}$
$\displaystyle\text{E}\left(e^{-u_{1}\mu_{1}(C)-u_{2}\mu_{2}(C)}\mu_{1}(C)\right)\text{E}\left(e^{-u_{1}\mu_{1}(D)-u_{2}\mu_{2}(D)}\mu_{2}(D)\right)\times$
$\displaystyle\times\text{E}\left(e^{-u_{1}\mu_{1}(Y)-u_{2}\mu_{2}(Y)}\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}$
$\displaystyle\frac{\partial}{\partial
u_{1}}\left\\{e^{-G_{0}(C)\psi_{b}(u_{1},u_{2})}\right\\}\frac{\partial}{\partial
u_{2}}\left\\{e^{-G_{0}(D)\psi_{b}(u_{1},u_{2})}\right\\}e^{-G_{0}(Y)\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}$
$\displaystyle=G_{0}$
$\displaystyle(C)G_{0}(D)\int_{\mathbb{R}_{+}^{2}}\left\\{\frac{\partial}{\partial
u_{1}}\psi_{b}(u_{1},u_{2})\frac{\partial}{\partial
u_{2}}\psi_{b}(u_{1},u_{2})\right\\}e^{-\psi_{b}(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}$

The result follows by applying Lemma 1. ∎

###### Proof of Proposition 7.

We have

$\displaystyle\mathbb{P}\left(X\in A,Y\in
B\right)=\text{E}\left(\frac{\tilde{\mu}_{1}(A)}{\tilde{\mu}_{1}(\mathbb{X})}\frac{\tilde{\mu}_{2}(B)}{\tilde{\mu}_{2}(\mathbb{X})}\right)=\text{E}\left(\frac{\mu_{1}(A\times\mathbb{X})}{\mu_{1}(\mathbb{X}\times\mathbb{X})}\frac{\mu_{2}(\mathbb{X}\times
B)}{\mu_{2}(\mathbb{X}\times\mathbb{X})}\right)=$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\text{E}\left(e^{-u_{1}\mu_{1}(\mathbb{X}\times\mathbb{X})-u_{2}\mu_{2}(\mathbb{X}\times\mathbb{X})}\mu_{1}(A\times\mathbb{X})\mu_{2}(\mathbb{X}\times
B)\right)\,\mathrm{d}u_{1}\mathrm{d}u_{2}=$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\text{E}\biggl{(}e^{-u_{1}\mu_{1}(\mathbb{X}\times\mathbb{X})-u_{2}\mu_{2}(\mathbb{X}\times\mathbb{X})}\bigl{\\{}\mu_{1}(A\times
B)\mu_{2}(A\times B)+\mu_{1}(A\times B)\mu_{2}(A^{c}\times B)+$
$\displaystyle+\mu_{1}(A\times B^{c})\mu_{2}(A\times B)+\mu_{1}(A\times
B^{c})\mu_{2}(A^{c}\times
B)\bigr{\\}}\biggr{)}\,\mathrm{d}u_{1}\mathrm{d}u_{2}$

We compute each integral separately applying Lemmas 2 and 3 and obtain

$\displaystyle\mathbb{P}\left(X\in A,Y\in B\right)$
$\displaystyle=G_{0}(A\times\mathbb{X})G_{0}(\mathbb{X}\times
B)\left(1-\delta\right)+G_{0}(A\times B)\delta$ (8)
$\displaystyle=P_{0}(A)P_{0}(B)\left(1-\delta\right)+G_{0}(A\times B)\delta,$

as desired. Then the probability of a tie in the product space is given
exactly by $\delta$, denoted $\gamma$ in the manuscript. The probability of a
tie is given by the particular case $\psi_{b}(u_{1},u_{2})=\psi(u_{1}+u_{2})$,
since

$-\int_{\mathbb{R}_{+}^{2}}\left\\{\frac{\partial^{2}}{\partial u_{1}\partial
u_{2}}\psi_{b}(u_{1}+u_{2})\right\\}e^{-\psi_{b}(u_{1}+u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}=-\int_{0}^{\infty}\int_{0}^{u}\,\mathrm{d}v\left\\{\frac{\partial^{2}}{\partial
u^{2}}\psi_{b}(u)\right\\}e^{-\psi_{b}(u)}\,\mathrm{d}u,$

with the change of variables $u=u_{1}+u_{2}$ and $v=u_{1}$. ∎

###### Proof of Proposition 8.

Since

$\text{E}\left(\tilde{p}_{1}(A)\tilde{p}_{2}(B)\right)=\mathbb{P}\left(X\in
A,Y\in B\right),$

by (8) we have

$\text{E}\left(\tilde{p}_{1}(A)\tilde{p}_{2}(B)\right)=G_{0}(A\times\mathbb{X})G_{0}(\mathbb{X}\times
B)\left(1-\gamma\right)+G_{0}(A\times B)\gamma.$

Finally,

$\displaystyle\text{cov}\left(\tilde{p}_{1}(A),\tilde{p}_{2}(B)\right)$
$\displaystyle=G_{0}(A\times\mathbb{X})G_{0}(\mathbb{X}\times
B)\left(1-\gamma\right)+G_{0}(A\times B)\gamma-
G_{0}(A\times\mathbb{X})G_{0}(\mathbb{X}\times B)$
$\displaystyle=\gamma\left\\{G_{0}(A\times
B)-G_{0}(A\times\mathbb{X})G_{0}(\mathbb{X}\times B)\right\\}.$

From this one also obtains

$\displaystyle\text{\mbox{var}}\left(\tilde{p}_{1}(A)\right)$
$\displaystyle=\text{cov}\left(\tilde{p}_{1}(A),\tilde{p}_{1}(A)\right)=\beta\left\\{P_{0}(A)-P_{0}(A)^{2}\right\\}$
$\displaystyle=\beta P_{0}(A)\left\\{1-P_{0}(A)\right\\},$

and the desired result follows. ∎

###  S1.4 Proofs of Section 5

###### Proof of Theorem 1.

We need to compute the conditional Laplace functional of $(\mu_{1},\mu_{2})$,
i.e.

$\text{E}\left(e^{-\int_{\mathbb{X}^{2}}h_{1}(x)\,\mu_{1}(\mathrm{d}x)-\int_{\mathbb{X}^{2}}h_{2}(x)\,\mu_{2}(\mathrm{d}x)}\mid(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m}\right),$

with $h_{i}\,:\mathbb{X}^{2}\,\to\mathbb{R}^{+}$ measurable functions. Define
$A_{j}=A_{j,\epsilon}=\left\\{x\in\mathbb{X}\mid
d(x,X_{i}^{*})<\epsilon\right\\}$ and
$B_{j}=B_{j,\epsilon}=\left\\{x\in\mathbb{X}\mid
d(x,Y_{j}^{*})<\epsilon\right\\}$, with $1\leq i\leq k$ and $1\leq j\leq c$,
such that $A_{i}\cap A_{j}=\emptyset$ and $B_{i}\cap B_{j}=\emptyset$ for any
$i\neq j$. Moreover, denote

$A_{k+1}=\left(\cup_{i=1}^{k}A_{i}\right)^{c},\quad
B_{c+1}=\left(\cup_{i=1}^{c}B_{i}\right)^{c}.$

Thus our goal becomes to compute

$\displaystyle\text{E}\left(e^{-\int_{\mathbb{X}^{2}}h_{1}(x)\,\mu_{1}(\mathrm{d}x)-\int_{\mathbb{X}^{2}}h_{2}(x)\,\mu_{2}(\mathrm{d}x)}\mid(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m}\right)$
(9) $\displaystyle=\lim_{\epsilon\to
0}\text{E}\left(e^{-\int_{\mathbb{X}^{2}}h_{1}(x)\,\mu_{1}(\mathrm{d}x)-\int_{\mathbb{X}^{2}}h_{2}(x)\,\mu_{2}(\mathrm{d}x)}\mid\underline{X}_{n}^{*}\in\times_{j=1}^{k}A_{j},\underline{Y}_{m}^{*}\in\times_{j=1}^{c}B_{j}\right)$
$\displaystyle=\lim_{\epsilon\to
0}\frac{\text{E}\left(e^{-\int_{\mathbb{X}^{2}}h_{1}(x)\,\mu_{1}(\mathrm{d}x)-\int_{\mathbb{X}^{2}}h_{2}(x)\,\mu_{2}(\mathrm{d}x)}\prod_{j=1}^{k}\tilde{p}_{1}(A_{j})^{n_{j}}\prod_{j=1}^{c}\tilde{p}_{2}(B_{j})^{m_{j}}\right)}{\text{E}\left(\prod_{j=1}^{k}\tilde{p}_{1}(A_{j})^{n_{j}}\prod_{j=1}^{c}\tilde{p}_{2}(B_{j})^{m_{j}}\right)}.$

We start to evaluate

$\displaystyle\text{E}\left(\tilde{p}_{1}(A_{1})^{n_{1}}\dots\tilde{p}_{1}(A_{k})^{n_{k}}\tilde{p}_{2}(B_{1})^{m_{1}}\tilde{p}_{2}(B_{c})^{m_{c}}\right)=$
$\displaystyle=\text{E}\left(\frac{\tilde{\mu}_{1}(A_{1})^{n_{1}}\dots\tilde{\mu}_{1}(A_{k})^{n_{k}}\tilde{\mu}_{2}(B_{1})^{m_{1}}\tilde{\mu}_{2}(B_{c})^{m_{c}}}{\tilde{\mu}_{1}(\mathbb{X})^{n}\tilde{\mu}_{2}(\mathbb{X})^{m}}\right)$
$\displaystyle=\text{E}\left(\frac{\mu_{1}(A_{1}\times\mathbb{X})^{n_{1}}\dots\mu_{1}(A_{k}\times\mathbb{X})^{n_{k}}\mu_{2}(\mathbb{X}\times
B_{1})^{m_{1}}\mu_{2}(\mathbb{X}\times
B_{c})^{m_{c}}}{\mu_{1}(\mathbb{X}\times\mathbb{X})^{n}\mu_{2}(\mathbb{X}\times\mathbb{X})^{m}}\right)=\mathcal{I}.$

By Netwon’s binomial

$\displaystyle\mu_{1}(A_{h}\times\mathbb{X})=\sum_{i_{1}^{h}+\dots
i_{c+1}^{h}=n_{h}}\binom{n_{h}}{i_{1}^{h},\dots,i_{c+1}^{h}}\prod_{r=1}^{c+1}\mu_{1}^{i_{r}^{h}}(A_{h}\times
B_{r}),\quad h=1,\dots,k,$ $\displaystyle\mu_{2}(\mathbb{X}\times
B_{r})=\sum_{j_{1}^{r}+\dots
j_{k+1}^{r}=m_{r}}\binom{m_{r}}{j_{1}^{r},\dots,j_{k+1}^{r}}\prod_{h=1}^{k+1}\mu_{2}^{j_{h}^{r}}(A_{h}\times
B_{r}),\quad r=1,\dots,c.$

For ease of notation denote

$\displaystyle\sum_{\bm{i},\bm{j}}\binom{\bm{n}}{\bm{i}}\binom{\bm{m}}{\bm{j}}=$
$\displaystyle\sum_{i_{1}^{1}+\dots
i_{c+1}^{1}=n_{1}}\binom{n_{1}}{i_{1}^{1},\dots,i_{c+1}^{1}}\dots\sum_{i_{1}^{c+1}+\dots
i_{c+1}^{k+1}=n_{k+1}}\binom{n_{k+1}}{i_{1}^{k+1},\dots,i_{c+1}^{k+1}}\times$
$\displaystyle\times\sum_{j_{1}^{1}+\dots
j_{k+1}^{1}=m_{1}}\binom{m_{1}}{j_{1}^{1},\dots,j_{k+1}^{1}}\dots\sum_{j_{1}^{k+1}+\dots
j_{k+1}^{k+1}=m_{k+1}}\binom{m_{k+1}}{j_{1}^{k+1},\dots,j_{k+1}^{k+1}}.$

Thus

$\mathcal{I}=\sum_{\bm{i},\bm{j}}\binom{\bm{n}}{\bm{i}}\binom{\bm{m}}{\bm{j}}\mathcal{I}_{\bm{i},\bm{j}},$

with

$\displaystyle\mathcal{I}_{\bm{i},\bm{j}}=\text{E}\Biggl{(}$
$\displaystyle\frac{\prod_{h=1}^{k}\prod_{r=1}^{c}\mu_{1}^{i_{r}^{h}}(A_{h}\times
B_{r})\mu_{2}^{j_{h}^{r}}(A_{h}\times
B_{r})}{\mu_{1}(\mathbb{X}\times\mathbb{X})^{n}}\times$
$\displaystyle\times\frac{\prod_{h=1}^{k}\mu_{1}^{i_{c+1}^{h}}(A_{h}\times
B_{c+1})\prod_{r=1}^{c}\mu_{2}^{j_{k+1}^{r}}(A_{k+1}\times
B_{r})}{\mu_{2}(\mathbb{X}\times\mathbb{X})^{m}}\Biggr{)}$

Letting $\mu_{1}:=\mu_{1}(\mathbb{X}\times\mathbb{X})$ and
$\mu_{2}:=\mu_{2}(\mathbb{X}\times\mathbb{X})$, we have

$\frac{1}{\mu_{1}(\mathbb{X}\times\mathbb{X})^{n}\mu_{2}(\mathbb{X}\times\mathbb{X})^{m}}=\frac{1}{\Gamma(n)\Gamma(m)}\int_{\mathbb{R}_{+}^{2}}u_{1}^{n-1}u_{2}^{m-1}e^{-u_{1}\mu_{1}-u_{2}\mu_{2}}\,\mathrm{d}\underline{u},$

with $\underline{u}=(u_{1},u_{2})$. Thus, by Fubini’s Theorem

$\displaystyle\mathcal{I}_{\textbf{i},\textbf{j}}=\int_{\mathbb{R}_{+}^{2}}\frac{u_{1}^{n-1}u_{2}^{m-1}}{\Gamma(n)\Gamma(m)}$
$\displaystyle\text{E}\biggl{(}e^{-u_{1}\mu_{1}-u_{2}\mu_{2}}\left\\{\prod_{h=1}^{k}\prod_{r=1}^{c}\mu_{1}^{i_{r}^{h}}(A_{h}\times
B_{r})\mu_{2}^{j_{h}^{r}}(A_{h}\times B_{r})\right\\}\times$
$\displaystyle\times\prod_{h=1}^{k}\mu_{1}^{i_{c+1}^{h}}(A_{h}\times
B_{c+1})\prod_{r=1}^{c}\mu_{2}^{j_{k+1}^{r}}(A_{k+1}\times
B_{r})\biggr{)}\,\mathrm{d}\underline{u}=$
$\displaystyle=\int_{\mathbb{R}_{+}^{2}}\frac{u_{1}^{n-1}u_{2}^{m-1}}{\Gamma(n)\Gamma(m)}\rho_{\textbf{i},\textbf{j}}(\underline{u})\,\mathrm{d}\underline{u}.$

By independence of evaluations on disjoint sets we have

$\displaystyle\rho_{\textbf{i},\textbf{j}}(\underline{u})$
$\displaystyle=\text{E}\Biggl{(}\left\\{\prod_{h=1}^{k}\prod_{r=1}^{c}e^{-u_{1}\mu_{1}(A_{h}\times
B_{r})-u_{2}\mu_{2}(A_{h}\times B_{r})}\mu_{1}^{i_{r}^{h}}(A_{h}\times
B_{r})\mu_{2}^{j_{h}^{r}}(A_{h}\times B_{r})\right\\}\times$
$\displaystyle\times\left\\{\prod_{h=1}^{k}e^{-u_{1}\mu_{1}(A_{h}\times
B_{c+1})-u_{2}\mu_{2}(A_{h}\times B_{c+1})}\mu_{1}^{i_{c+1}^{h}}(A_{h}\times
B_{c+1})\right\\}\times$
$\displaystyle\times\left\\{\prod_{r=1}^{c}e^{-u_{1}\mu_{1}(A_{k+1}\times
B_{r})-u_{2}\mu_{2}(A_{k+1}\times B_{r})}\mu_{2}^{j_{k+1}^{r}}(A_{k+1}\times
B_{r})\right\\}\Biggr{)}$

This can be equivalently written as

$\displaystyle\prod_{h=1}^{k}\prod_{r=1}^{c}$
$\displaystyle\text{E}\left(e^{-u_{1}\mu_{1}(A_{h}\times
B_{r})-u_{2}\mu_{2}(A_{h}\times B_{r})}\mu_{1}^{i_{r}^{h}}(A_{h}\times
B_{r})\mu_{2}^{j_{h}^{r}}(A_{h}\times B_{r})\right)\times$
$\displaystyle\times\prod_{h=1}^{k}\text{E}\left(e^{-u_{1}\mu_{1}(A_{h}\times
B_{c+1})-u_{2}\mu_{2}(A_{h}\times B_{c+1})}\mu_{1}^{i_{c+1}^{h}}(A_{h}\times
B_{c+1})\right)\times$
$\displaystyle\times\prod_{r=1}^{c}\text{E}\left(e^{-u_{1}\mu_{1}(A_{k+1}\times
B_{r})-u_{2}\mu_{2}(A_{k+1}\times B_{r})}\mu_{2}^{j_{k+1}^{r}}(A_{k+1}\times
B_{r})\right).$

Considering each element separately we have

$\displaystyle\text{E}\biggl{(}e^{-u_{1}\mu_{1}(A_{h}\times
B_{r})-u_{2}\mu_{2}(A_{h}\times B_{r})}$ $\displaystyle\mu_{1}^{i}(A_{h}\times
B_{r})\mu_{2}^{j}(A_{h}\times B_{r})\biggr{)}$
$\displaystyle=\text{E}\left((-1)^{i+j}\frac{\partial^{i+j}}{\partial
u_{1}^{i}\partial u_{2}^{j}}e^{-u_{1}\mu_{1}(A_{h}\times
B_{r})-u_{2}\mu_{2}(A_{h}\times B_{r})}\right)$
$\displaystyle=(-1)^{i+j}\frac{\partial^{i+j}}{\partial u_{1}^{i}\partial
u_{2}^{j}}\text{E}\left(e^{-u_{1}\mu_{1}(A_{h}\times
B_{r})-u_{2}\mu_{2}(A_{h}\times B_{r})}\right)$
$\displaystyle=(-1)^{i+j}\frac{\partial^{i+j}}{\partial u_{1}^{i}\partial
u_{2}^{j}}\left\\{e^{-\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(ds)G_{0}(x)}\right\\}.$

Recall that we are interested in the limit as $\epsilon\to 0$, so that

$\displaystyle\frac{\partial^{i+j}}{\partial u_{1}^{i}\partial
u_{2}^{j}}\biggl{\\{}$ $\displaystyle e^{-\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}\biggr{\\}}\sim
e^{-\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}\times$
(10) $\displaystyle\times\frac{\partial^{i+j}}{\partial u_{1}^{i}\partial
u_{2}^{j}}\left\\{\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\},$

where we say $f\sim g$ if $\lim_{\epsilon\to 0}f(x)/g(x)=1$. By simple algebra
we get

$\displaystyle\frac{\partial^{i+j}}{\partial u_{1}^{i}\partial u_{2}^{j}}$
$\displaystyle\left\\{e^{-\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}\right\\}=\frac{\partial^{i+j-1}}{\partial
u_{1}^{i-1}\partial u_{2}^{j}}\biggl{\\{}-\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}\times$
$\displaystyle\times
s_{1}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)e^{-\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}\biggr{\\}}$
$\displaystyle=\frac{\partial^{i+j-2}}{\partial u_{1}^{i-2}\partial
u_{2}^{j}}\biggl{\\{}\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}^{2}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)e^{-\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}$
$\displaystyle+\left(\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right)^{2}e^{-\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}\biggr{\\}},$

and

$\lim_{\epsilon\to 0}\frac{\left(\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right)^{2}}{\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}^{2}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}=0.$

By applying this argument repeatedly we obtain (10). Thus, letting
$\rho({\underline{u}})=\sum_{\bm{i},\bm{j}}\binom{\bm{n}}{\bm{i}}\binom{\bm{m}}{\bm{j}}\rho_{\bm{i},\bm{j}}({\underline{u}})$,
by aggregating the terms we have

$\displaystyle\rho({\underline{u}})\sim$
$\displaystyle\sum_{\bm{i},\bm{j}}\binom{\bm{n}}{\bm{i}}\binom{\bm{m}}{\bm{j}}(-1)^{n+m}e^{-\psi_{b}(u)}\times$
$\displaystyle\times\prod_{h=1}^{k}\prod_{r=1}^{c}\left\\{\frac{\partial^{i_{r}^{h}+j_{h}^{r}}}{\partial
u_{1}^{i_{r}^{h}}\partial u_{2}^{j_{h}^{r}}}\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\times$
$\displaystyle\times\prod_{h=1}^{k}\left\\{\frac{\partial^{i_{c+1}^{h}}}{\partial
u_{1}^{i^{h}_{c+1}}}\int_{A_{h}\times
B_{c+1}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\times$
$\displaystyle\times\prod_{r=1}^{c}\left\\{\frac{\partial^{j_{k+1}^{r}}}{\partial
u_{2}^{i^{r}_{k+1}}}\int_{A_{k+1}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}$
$\displaystyle=\sum_{\bm{i},\bm{j}}\binom{\bm{n}}{\bm{i}}\binom{\bm{m}}{\bm{j}}(-1)^{n+m}V(\bm{i},\bm{j}).$

The following three Lemmas characterize the set of indices $(\bm{i},\bm{j})$
that are relevant once the limit is taken.

######  Lemma 4.

Consider $(\bm{i},\bm{j})$ such that $0<i_{r}^{h},i_{l}^{h}<n_{h}$, with $r>l$
and $1\leq h\leq k$. Then $\exists(\tilde{\bm{i}},\tilde{\bm{j}})$ such that
$\lim_{\epsilon\to 0}V(\bm{i},\bm{j})/V(\tilde{\bm{i}},\tilde{\bm{j}})\to 0$.

###### Proof of Lemma S2.4.

For ease of notation set $\bm{i}^{h}=(i^{h}_{1},\dots,i^{h}_{c+1})$. Then

  * •

If $r=c+1$, set
$\tilde{\bm{i}}^{h}=(i_{1}^{h},\dots,i_{l}^{h}+i^{h}_{c+1},\dots,0)$.

  * •

If $j_{h}^{r}=0$, set
$\tilde{i}^{h}=(i_{1}^{h},\dots,i_{l}^{h}+i_{r}^{h},\dots,0,\dots)$.

  * •

If $j_{h}^{l}=0$, set
$\tilde{\bm{i}}^{h}=(i_{1}^{h},\dots,0,\dots,i_{r}^{h}+i_{l}^{h},\dots)$.

  * •

If $j_{h}^{l}>0$ and $j_{h}^{r}>0$, set
$\tilde{\bm{j}}^{r}=(j_{1}^{r},\dots,0,\dots,j^{r}_{k+1}+j_{h}^{r})$ and
$\tilde{\bm{i}}^{h}=(i_{1}^{h},\dots,i_{l}^{h}+i_{r}^{h},\dots,0,\dots)$.

For example in the last case we have

$\lim_{\epsilon\to
0}\frac{\mbox{var}(\bm{i},\bm{j})}{\mbox{var}(\tilde{\bm{i}},\tilde{\bm{j}})}=\lim_{\epsilon\to
0}\frac{\int_{A_{h}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}^{i^{h}_{r}}s_{2}^{j^{r}_{h}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}{\int_{A_{c+1}\times
B_{r}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{2}^{j^{r}_{h}+j^{r}_{c+1}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}=0,$

as desired. ∎

Thus, Lemma 4 guarantees that $\bm{i}^{h}$ has exactly one element different
from $0$, that is equal to $n_{h}$.

######  Lemma 5.

Consider $(\textbf{i},\textbf{j})$ such that $i_{r}^{h}=n_{h}$ and
$j_{h}^{r}=0$. Then there exists $(\tilde{\textbf{i}},\tilde{\textbf{j}})$
such that

$\lim_{\epsilon\to
0}V(\textbf{i},\textbf{j})/V(\tilde{\textbf{i}},\tilde{\textbf{j}})\to 0.$

###### Proof of Lemma S5.

Set $(\tilde{\bm{i}},\tilde{\bm{j}})$ equal to $(\textbf{i},\textbf{j})$,
apart from $\tilde{i}^{h}_{r}=0$ and $\tilde{i}^{h}_{c+1}=n_{h}$. ∎

######  Lemma 6.

Consider $(\textbf{i},\textbf{j})$ such that $i^{h}_{c+1}=n_{h}$ and
$j_{h}^{r}>0$. Then there exists $(\tilde{\textbf{i}},\tilde{\textbf{j}})$
such that

$\lim_{\epsilon\to
0}V(\textbf{i},\textbf{j})/V(\tilde{\textbf{i}},\tilde{\textbf{j}})\to 0.$

###### Proof of Lemma S2.6.

Set $(\tilde{\bm{i}},\tilde{\bm{j}})$ equal to $(\bm{i},\bm{j})$, apart from
$\tilde{j}_{h}^{r}=0$ and $\tilde{j}_{k+1}^{r}=m_{r}$. ∎

The three lemmas imply that each relevant $(\bm{i},\bm{j})$ corresponds to an
admissible latent structure, i.e.

$\displaystyle\rho({u})\sim\sum_{\textbf{p}\in\mathcal{P}}(-1)^{n+m}e^{-\psi_{b}(\underline{u})}$
$\displaystyle\prod_{(i,j)\in\Delta_{\textbf{p}}}\left\\{\frac{\partial^{n_{i}+m_{j}}}{\partial
u_{1}^{n_{i}}\partial u_{2}^{m_{j}}}\int_{A_{i}\times
B_{j}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\times$
$\displaystyle\times\prod_{(i,j)\in\Delta^{1}_{\textbf{p}}}\left\\{\frac{\partial^{n_{i}}}{\partial
u_{1}^{n_{i}}}\int_{A_{i}\times
B_{c+1}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\times$
$\displaystyle\times\prod_{(i,j)\in\Delta^{2}_{\textbf{p}}}\left\\{\frac{\partial^{m_{j}}}{\partial
u_{2}^{m_{j}}}\int_{A_{k+1}\times
B_{j}}\int_{\mathbb{R}_{+}^{2}}(1-e^{-u_{1}s_{1}-u_{2}s_{2}})\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}.$

Evaluating the derivatives we have

$\displaystyle\rho({\underline{u}})\sim\sum_{\textbf{p}\in\mathcal{P}}e^{-\psi_{b}({\underline{u}})}$
$\displaystyle\prod_{(i,j)\in\Delta_{\textbf{p}}}\left\\{\int_{A_{i}\times
B_{j}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}^{n_{i}}s_{2}^{m_{j}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\times$
$\displaystyle\times\prod_{(i,j)\in\Delta^{1}_{\textbf{p}}}\left\\{\int_{A_{i}\times
B_{c+1}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}^{n_{i}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\times$
$\displaystyle\times\prod_{(i,j)\in\Delta^{2}_{\textbf{p}}}\left\\{\int_{A_{k+1}\times
B_{j}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{2}^{m_{j}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}.$

Finally, we get

$\displaystyle\mathcal{I}\sim\sum_{\textbf{p}\in\mathcal{P}}\int_{\mathbb{R}_{+}^{2}}\frac{u_{1}^{n-1}u_{2}^{m-1}}{\Gamma(n)\Gamma(m)}e^{-\psi_{b}({\underline{u}})}$
$\displaystyle\prod_{(i,j)\in\Delta_{\textbf{p}}}\left\\{\int_{A_{i}\times
B_{j}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}^{n_{i}}s_{2}^{m_{j}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\times$
$\displaystyle\times\prod_{(i,j)\in\Delta^{1}_{\textbf{p}}}\left\\{\int_{A_{i}\times
B_{c+1}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{1}^{n_{i}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\times$
$\displaystyle\times\prod_{(i,j)\in\Delta^{2}_{\textbf{p}}}\left\\{\int_{A_{k+1}\times
B_{j}}\int_{\mathbb{R}_{+}^{2}}e^{-u_{1}s_{1}-u_{2}s_{2}}s_{2}^{m_{j}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\,\mathrm{d}{u}.$

Evaluating the numerator of (9) the same reasoning yields a formula asymptotic
to

$\displaystyle\sum_{\textbf{p}\in\mathcal{P}}\int_{\mathbb{R}_{+}^{2}}\frac{u_{1}^{n-1}u_{2}^{m-1}}{\Gamma(n)\Gamma(m)}e^{-\psi_{h}({\underline{u}})}$
$\displaystyle\prod_{(i,j)\in\Delta_{\textbf{p}}}\left\\{\int_{A_{i}\times
B_{j}}\int_{\mathbb{R}_{+}^{2}}e^{-(h_{1}(x)+u_{1})s_{1}-(h_{2}(x)+u_{2})s_{2}}s_{1}^{n_{i}}s_{2}^{m_{j}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}$
$\displaystyle\prod_{(i,j)\in\Delta^{1}_{\textbf{p}}}\left\\{\int_{A_{i}\times
B_{c+1}}\int_{\mathbb{R}_{+}^{2}}e^{-(h_{1}(x)+u_{1})s_{1}-(h_{2}(x)+u_{2})s_{2}}s_{1}^{n_{i}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}$
$\displaystyle\prod_{(i,j)\in\Delta^{2}_{\textbf{p}}}\left\\{\int_{A_{k+1}\times
B_{j}}\int_{\mathbb{R}_{+}^{2}}e^{-(h_{1}(x)+u_{1})s_{1}-(h_{2}(x)+u_{2})s_{2}}s_{2}^{m_{j}}\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)\right\\}\,\mathrm{d}{u}.$

where
$\psi_{h}({\underline{u}})=\int_{\mathbb{X}^{2}}\int_{\mathbb{R}_{+}^{2}}\left(1-e^{-(h_{1}(x)+u_{1})s_{1}-(h_{2}(x)+u_{2})s_{2}}\right)\,\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)$.
Note that

$\displaystyle 1-e^{-(h_{1}(x)+u_{1})s_{1}-(h_{2}(x)+u_{2})s_{2}}$
$\displaystyle=e^{-u_{1}s_{1}-u_{2}s_{2}}\left[e^{u_{1}s_{1}+u_{2}s_{2}}-1+1-e^{-h_{1}(x)s_{1}-h_{2}(x)s_{2}}\right]$
$\displaystyle=\left[1-e^{-u_{1}s_{1}-u_{2}s_{2}}\right]+\left[1-e^{-h_{1}(x)s_{1}-h_{2}(x)s_{2}}\right],$

so that

$\displaystyle e^{-\psi_{h}({\underline{u}})}$
$\displaystyle=e^{-\psi_{b}({\underline{u}})}e^{-\int_{\mathbb{X}^{2}}\int_{\mathbb{R}_{+}^{2}}\left[1-e^{-h_{1}(x)s_{1}-h_{2}(x)s_{2}}\right]\rho(\mathrm{d}s)G_{0}(\mathrm{d}x)}$
$\displaystyle=e^{-\psi_{b}({\underline{u}})}\text{E}\left[e^{-\int_{\mathbb{X}^{2}}h_{1}(x)\,\hat{\mu}_{1}(\mathrm{d}x)-\int_{\mathbb{X}^{2}}h_{2}(x)\,\hat{\mu}_{2}(\mathrm{d}x)}\right].$

Furthermore

$G_{0}(A_{h}\times B_{r})=\epsilon\frac{G_{0}(A_{h}\times
B_{r})}{\epsilon}\sim\epsilon g_{h,r},\quad 1\leq i\leq c,1\leq j\leq k,$

and

$G_{0}(A_{h}\times\mathrm{d}x)\sim\epsilon
g_{h,c+1}Q_{X_{h}^{*}}(\mathrm{d}x),\quad G_{0}(\mathrm{d}x\times
B_{r})\sim\epsilon g_{k+1,r}P_{Y_{r}^{*}}(\mathrm{d}x).$

Thus, evaluating the limit in (9) we get

$\displaystyle\text{E}\biggl{[}$ $\displaystyle
e^{-\int_{\mathbb{X}^{2}}h_{1}(x)\,\mu_{1}(\mathrm{d}x)-\int_{\mathbb{X}^{2}}h_{2}(x)\,\mu_{2}(\mathrm{d}x)}\mid(X_{i})_{i\geq
1}^{n},(Y_{j})_{j\geq 1}^{m}\biggr{]}=$
$\displaystyle\times\sum_{p\in\mathcal{P}}\int_{\mathbb{R}_{+}^{2}}\text{E}\left[e^{-\int_{\mathbb{X}^{2}}h_{1}(x)\,\hat{\mu}_{1}(\mathrm{d}x)-\int_{\mathbb{X}^{2}}h_{2}(x)\,\hat{\mu}_{2}(\mathrm{d}x)}\right]\times$
$\displaystyle\times\prod_{(i,j)\in\Delta_{p}}\int_{\mathbb{R}_{+}^{2}}e^{-h_{1}(X_{i}^{*},Y_{j}^{*})s_{1}-h_{2}(X_{i}^{*},Y_{j}^{*})s_{2}}\frac{s_{1}^{n_{i}}s_{2}^{m_{j}}e^{-u_{1}s_{1}-u_{2}s_{2}}\rho(\mathrm{d}s)}{\tau_{n_{i},m_{j}(\underline{u})}}\times$
$\displaystyle\times\prod_{(i,j)\in\Delta^{1}_{p}}\int_{\mathbb{X}}\int_{\mathbb{R}_{+}^{2}}e^{-h_{1}(X_{i}^{*},x_{2})s_{1}-h_{2}(X_{i}^{*},x_{2})s_{2}}\frac{s_{1}^{n_{i}}e^{-u_{1}s_{1}-u_{2}s_{2}}\rho(\mathrm{d}s)}{\tau_{n_{i},0(\underline{u})}}Q_{X_{i}^{*}}(\mathrm{d}x_{2})\times$
$\displaystyle\times\prod_{(i,j)\in\Delta^{2}_{p}}\int_{\mathbb{X}}\int_{\mathbb{R}_{+}^{2}}e^{-h_{1}(x_{1},Y_{2}^{*})s_{1}-h_{2}(x_{1},Y_{2}^{*})s_{2}}\frac{s_{2}^{m_{j}}e^{-u_{1}s_{1}-u_{2}s_{2}}\rho(\mathrm{d}s)}{\tau_{0,m_{j}(\underline{u})}}P_{Y_{j}^{*}}(\mathrm{d}x_{1})\times$
$\displaystyle\times\left(\frac{\int_{\mathbb{R}_{+}^{2}}u_{1}^{n-1}u_{2}^{m-1}\prod_{(i,j)\in
p}g_{i,j}\tau_{n_{i},m_{j}}(\underline{u})e^{-\psi_{b}(\underline{u})}\,\mathrm{d}\underline{u}}{\sum_{\textbf{q}\in\mathcal{P}}\int_{\mathbb{R}_{+}^{2}}u_{1}^{n-1}u_{2}^{m-1}\prod_{(i,j)\in\textbf{q}}g_{i,j}\tau_{n_{i},m_{j}}(\underline{u})e^{-\psi_{b}(\underline{u})}\mathrm{d}\underline{u}}\right)\times$
$\displaystyle\times\frac{u_{1}^{n-1}u_{2}^{m-1}\prod_{(i,j)\in
p}\tau_{n_{i},m_{j}}(\underline{u})e^{-\psi_{b}(\underline{u})}\,\mathrm{d}u}{\int_{\mathbb{R}_{+}^{2}}u_{1}^{n-1}u_{2}^{m-1}\prod_{(i,j)\in
p}\tau_{n_{i},m_{j}}(\underline{u})e^{-\psi_{b}(\underline{u})}\,\mathrm{d}\underline{u}},$

as desired. ∎

###### Proof of Corollary 2.

We use the shorthand notation
$\mu_{1}(f)=\int_{\mathbb{X}}f(x)\,\mu_{1}(\mathrm{d}x)$ for any measurable
function $f\,:\,\mathbb{X}\,\to\,\mathbb{R}$ such that $\mu_{1}(|f|)<\infty$.
Letting ${U}$ be the set of latent variables of Theorem 1, i.e.
${U}=\left(p,U_{1},U_{2},Z^{x},Z^{y}\right)$ for any
$y_{1},\dots,y_{n}\in(0,1)$ and $A_{1},\dots,A_{n}\in\mathcal{X}^{2}$ we get

$\displaystyle\mathbb{P}\left[p_{1}(A_{1})\leq y_{1},\dots,p_{n}(A_{n})\leq
y_{n}\mid{U},(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m}\right]$ $\displaystyle=$
$\displaystyle\mathbb{P}\left[\mu_{1}(\mathbbm{1}_{A_{1}}-y_{1})\leq
0,\dots,\mu_{n}(\mathbbm{1}_{A_{n}}-y_{n})\leq
0\mid{U},(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m}\right].$

The result follows since the finite dimensional distributions of $p_{1}$ given
${U}$, $(X_{i})_{i=1}^{n}$, and $(Y_{j})_{j=1}^{m}$ coincide with the ones of
the normalized posterior distribution of $\mu_{1}$, given ${U}$,
$(X_{i})_{i=1}^{n}$, and $(Y_{j})_{j=1}^{m}$. ∎

###### Proof of Theorem 2.

Set ${\underline{H}}=\left(p,U_{1},U_{2}\right)$ with domain $D$. Then

$\displaystyle\mathbb{P}(X_{n+1}\in\mathrm{d}x\mid(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m})$
$\displaystyle=\text{E}[\tilde{p}_{1}(\mathrm{d}x)\mid(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m}]$
$\displaystyle=\int_{D}\text{E}[\tilde{p}_{1}(\mathrm{d}x)\mid\underline{H}=\underline{h},(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m}]\,F(\mathrm{d}\underline{v}),$

where $F(\cdot)$ is the posterior distribution of $\underline{H}$, with
$\underline{h}=(p,u_{1},u_{2})$. Recalling the notation in Corollary $2$ we
have

$\displaystyle\text{E}[\tilde{p}_{1}(\mathrm{d}x)\mid$
$\displaystyle\underline{H}=\underline{h},(X_{i})_{i=1}^{n},(Y_{j})_{j=1}^{m}]=\text{E}\left[\frac{\hat{\mu}_{1}(\mathrm{d}x\times\mathbb{X})}{R}\right]+\text{E}\left[\frac{\sum_{(i,j)\in\Delta_{p}}J^{1}_{i,j}\delta_{X_{i}^{*}}}{R}\right]+$
$\displaystyle+\text{E}\left[\frac{\sum_{(i,j)\in\Delta^{1}_{p}}J^{1}_{i,c+1}\delta_{X_{i}^{*}}}{R}\right]+\text{E}\left[\frac{\sum_{(i,j)\in\Delta^{2}_{p}}J^{1}_{k+1,j}\delta_{Z_{j}^{y}}}{R}\right]=\sum_{k=1}^{4}I_{k},$

where
$R=T_{1}+\sum_{(i,j)\in\Delta_{p}}J_{i,j}^{1}+\sum_{(i,j)\in\Delta^{1}_{p}}J_{i,c+1}^{1}+\sum_{(i,j)\in\Delta^{2}_{p}}J_{k+1,j}^{1}$.
Set
$S=\sum_{(i,j)\in\Delta_{p}}J_{i,j}^{1}+\sum_{(i,j)\in\Delta^{1}_{p}}J_{i,c+1}^{1}+\sum_{(i,j)\in\Delta^{2}_{p}}J_{k+1,j}^{1}$
and exploit the conditional independence between $J_{ij}^{1}$ and
$\hat{\mu}_{1}$ to obtain

$\displaystyle I_{1}$
$\displaystyle=\int_{\mathbb{R}_{+}}\text{E}\left[e^{-vS}\right]\text{E}\left[\hat{\mu}_{1}(\mathrm{d}x\times\mathbb{X})e^{-vT_{1}}\right]\,\mathrm{d}v$
$\displaystyle=\theta
P_{0}(\mathrm{d}x)\int_{\mathbb{R}_{+}}\left(\prod_{(i,j)\in
p}\frac{\tau_{n_{i},m_{j}}(u_{1}+v,u_{2})}{\tau_{n_{i},m_{j}}(u_{1},u_{2})}\right)\tau_{1,0}(u_{1}+v,u_{2})e^{-\psi^{\underline{u}}_{b}(v,0)}\,\mathrm{d}v,$

where $\psi_{b}^{\underline{u}}(\lambda_{1},\lambda_{2})$ is the Laplace
exponent of $(\hat{\mu}_{1},\hat{\mu}_{2})$ in Theorem 1. Observing that
$\psi_{b}^{\underline{u}}(v,0)+\psi(u_{1},u_{2})=\psi(u_{1}+v,u_{2})$ and
denoting with $L(\cdot)$ the distribution of p, we obtain

$\displaystyle\xi_{0}$
$\displaystyle=\int_{D}I_{1}\,F(\mathrm{d}\underline{u})$
$\displaystyle=\theta
P_{0}(\mathrm{d}x)\int\int_{\mathbb{R}_{+}^{3}}\biggl{\\{}u_{1}^{n-1}u_{2}^{m-1}\left(\prod_{(i,j)\in
p}\tau_{n_{i},m_{j}}(u_{1}+v,u_{2})\right)\tau_{1,0}(u_{1}+v,u_{2})\times$
$\displaystyle\times
e^{-\psi(u_{1}+v,u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}\mathrm{d}vL(\mathrm{d}p)\biggr{\\}}$
$\displaystyle=\frac{\theta
P_{0}(\mathrm{d}x)}{n}\int\int_{\mathbb{R}_{+}^{2}}u_{1}^{n}u_{2}^{m-1}\left(\prod_{(i,j)\in
p}\tau_{n_{i},m_{j}}(u_{1},u_{2})\right)\tau_{1,0}(u_{1},u_{2})e^{-\psi(u_{1},u_{2})}\,\mathrm{d}u_{1}\mathrm{d}u_{2}L(\mathrm{d}p)$
$\displaystyle=\frac{\theta
P_{0}(\mathrm{d}x)}{n}\int_{D}u_{1}\tau_{1,0}(u_{1},u_{2})\,F(\mathrm{d}\underline{u}),$

where the second equality follows from the change of variables
$(w,z)=(u_{1}+v,u_{1})$. The proof for the remaining weights follows along the
same lines and leads to

$\xi_{i}^{x}=\frac{1}{n}\int_{D}u_{1}\left[\frac{\tau_{n_{i}+1,m_{j}}(u_{1},u_{2})}{\tau_{n_{i},m_{j}}(u_{1},u_{2})}+\frac{\tau_{n_{i}+1,0}(u_{1},u_{2})}{\tau_{n_{i},0}(u_{1},u_{2})}\right]\,F(\mathrm{d}\underline{u})$