Enriched Pitman-Yor processes
Tommaso Rigon Bruno Scarpa Sonia Petrone
stat.ME math.ST stat.TH

Mutually Unbiased Unitary Bases of Operators on d-dimensional Hilbert Space
    Stefano Mancini
    December 30, 2023
===========================================================================




In Bayesian nonparametrics there exists a rich variety of discrete priors, including the Dirichlet process and its generalizations, which are nowadays well-established tools. Despite the remarkable advances, few proposals are tailored for modeling observations lying on product spaces, such as R^p. 
Indeed, for multivariate random measures, most available priors lack flexibility and do not allow for separate partition structures among the spaces. We introduce a discrete nonparametric prior, termed enriched Pitman–Yor process (epy), aimed at addressing these issues. Theoretical properties of this novel prior are extensively investigated. We discuss its formal link with the enriched Dirichlet process and normalized random measures, we describe a square-breaking representation and we obtain closed-form expressions for the posterior law and the involved urn schemes. 
In second place, we show that several existing approaches, including Dirichlet processes with a spike and slab base measure and mixture of mixtures models, implicitly rely on special cases of the epy, which therefore constitutes a unified probabilistic framework for many Bayesian nonparametric priors. Interestingly, our unifying formulation will allow us to naturally extend these models while preserving their analytical tractability. As an illustration, we employ the epy for a species sampling problem in ecology and for functional clustering in an e-commerce application. 



§ INTRODUCTION


The Dirichlet process (dp) of <cit.> is a widely employed nonparametric discrete prior which arguably stands because of its analytical tractability. Despite its popularity, the dp may be a restrictive modeling choice as it depends on a single parameter controlling both the variability and the random partition it induces. To overcome these limitations several proposals have been made, including the Pitman–Yor  process (py) discussed in <cit.>, <cit.>, and the very general classes of Gibbs-type priors <cit.>, normalized random measures with independent increments <cit.>, and species sampling models <cit.>. These extensions offer a richer modeling framework while preserving the tractability of the dp. One may refer to <cit.> for an overview.

Unfortunately, none of the aforementioned priors has been specifically designed to model exchangeable observations lying on a product space, such as R^p.  More precisely, let X and Y be two complete and separable Polish spaces and let ((X_n,Y_n))_n ≥ 1 be an infinite sequence of exchangeable random elements taking values in the product space X×Y, for example R^p=R^k ×R^p-k. Then, de Finetti theorem guarantees that conditionally on a random probability measure p̃ the random elements (X_n,Y_n) are independent and identically distributed (iid), namely

    (X_n, Y_n) |p̃   ∼p̃,      n ≥ 1,
    p̃   ∼𝒬,

where 𝒬 is the probability law of the random , namely the prior distribution. In a first motivating application, the random variables X_n and Y_n represent families and species of trees in the Amazonian basin, respectively. Given a collection of observations, we are then interested in predicting the number of novel families and species that one would get in a future sample. However, neither the dp, the py, or other species sampling models are suitable priors 𝒬. Indeed, these choices would associate the discovery of a new species with that of a new family, which is an unrealistic assumption. Instead, we seek priors inducing nested mechanisms, so that the discovery of a new family of trees corresponds to that of a new species, but not vice versa. In the second motivating application, we aim at forming groups of functional observations, each representing the number of searches on a website carried out over time by different customers. The random variables X_n and Y_n represent latent features and latent functions, respectively, which induce a nested and more interpretable clustering mechanism that will turn useful for market segmentation. 

A more suitable prior for the joint random probability measure p̃ in model (<ref>) is the enriched Dirichlet process (edp) of <cit.>. The edp allows for finer control of the dependence structure between the X_n and Y_n and leads to the desired nested clustering. The edp broadens the principles of enriched conjugate priors for natural exponential families <cit.> to the nonparametric setting since it extends the construction of the enriched Dirichlet distribution of <cit.> to random probability measures; refer to <cit.> for further details. These appealing features have been exploited among others in <cit.> for Bayesian nonparametric regression models, in <cit.> for causal inference with missing covariates, and in <cit.> for functional clustering of longitudinal data. 




Here we move away from the original derivation of the edp, considering constructions beyond conjugacy. The edp inherits also some of the drawbacks shared by all Dirichlet-based priors, which indeed motivates our extension. For instance, the probability of observing a new species in the urn scheme of the dp <cit.> solely depends on the sample size and not on the previously observed values. As remarked by <cit.>, this simplifying assumption is particularly problematic in species sampling problems, because the posterior probability of discovering a new species does not depend on the data. In addition, this feature of the dp leads to a logarithmic growth of the number of clusters and to a lack of robustness with respect to miscalibrated prior choices, which might be undesirable in several applied contexts; see e.g. <cit.> and <cit.>. 

We address these issues by proposing a novel discrete prior law 𝒬, that builds upon the edp and the py process. By combining their appealing properties, the proposed enriched Pitman–Yor process (epy) leads to different rates for the number of clusters and more robust Bayesian estimators for species sampling models. Importantly, improved flexibility is attained while preserving analytical and computational tractability. We obtain a simple urn scheme, a tractable posterior characterization, and a so-called square-breaking representation. In addition, we show that the epy can be defined by normalizing a suitable random measure. This alternative definition parallels the construction of <cit.> and has important theoretical implications. 

Beside their key role for the modeling of observations in product spaces, epy priors are provably useful also in other settings. Specifically, consider an exchangeable sequence (Y_n)_n ≥ 1 taking values in Y and let

    Y_n |p̃_Y    ∼p̃_Y,      n ≥ 1, 
    p̃_Y    ∼𝒬_Y.

Hence, one can consider a marginal epy process as the prior 𝒬_Y, which is defined as 

    p̃_Y(B) = p̃(X× B),     p̃∼𝒬,

for any Borel set B ⊆Y, where the prior law 𝒬 is an epy process on X×Y and the space X should be interpreted as a latent dimension that induces enriched specifications.  It will be shown that the marginal epy p̃_Y is an infinite mixture of py processes, which is arguably much more flexible than a single py process. 

We show that model (<ref>) includes a rich variety of prior proposals in the literature as a special case. To the best of our knowledge, their connection with enriched processes has not been previously emphasized. For example, a specific marginal epy process has been implicitly studied in <cit.> and <cit.> for the analysis of functional data. Dirichlet processes with spike and slab base measures <cit.>, or general atomic contaminations <cit.>, are actually special cases of a marginal edp. A mixture of finite-dimensional dps has been employed in <cit.> and <cit.> to perform model-based clustering, while convex combinations of dps have been considered by <cit.>, <cit.> to induce dependence across groups of random variables. These models are also strongly linked to the marginal epy. Through the paper, we will point out the connections between the epy and the aforementioned methods, aiming at providing a unified probabilistic framework for these classes of processes. 

The paper is organized as follows. In Section <ref> we introduce the epy process and we discuss its fundamental probabilistic characterizations, including the square-breaking construction. In Section <ref> we discuss an enriched urn scheme and posterior representations. In Section <ref> we illustrate the marginal epy process and its connection with several existing approaches, and we propose numerous extensions. In Section <ref> we employ the epy  to estimate the number of unobserved species in the Amazonian tree flora.  Finally, in Section <ref> we employ the epy in a mixture model for functional clustering, illustrating its usefulness for market segmentation. Concluding remarks are given in Section <ref>. All the proofs are collected in the Supplementary Material.



§ THE ENRICHED PITMAN–YOR PROCESS


The epy is built upon the dp and the py processes, of which we provide a concise overview that is also useful to set the notation. One can refer to <cit.> and <cit.> for more structured reviews. The epy process is then defined in Section <ref>, together with the so-called square-breaking representation. 



 §.§ Background material


The py process is a probability law on a random discrete distribution, that can be defined through the so-called stick-breaking construction. Let P be a probability measure on X×Y and let (ν_j)_j≥ 1 be a sequence of independent Beta random variables with ν_j ∼beta(1-σ,α + jσ), where either σ∈[0,1) and α> - σ, or σ = -α/H and α > 0 for some integer H ∈{2,3,…}. A discrete random probability measure p̃ follows a py process with parameters (σ, α P), written p̃∼py(σ, α P), if 

    p̃(·) = ∑_h = 1^∞ξ_h δ_(ϕ_h, θ_h)(·),      (ϕ_h,θ_h) ∼ P,

with ξ_h=ν_h∏_j=1^h-1(1-ν_j) for h≥ 1, where we agree that ξ_1 = ν_1. The parameter σ is often called discount or stable  parameter whereas α is termed total mass or precision. If σ = 0, then p̃ in equation (<ref>) defines a dp and we write p̃∼dp(α P). If instead σ = -α/H and α > 0, then the stick-breaking construction is degenerate because ν_H = 1, implying that ∑_h=1^Hξ_h =1 for a finite integer H < ∞. This special case of py, called Dirichlet multinomial process, or Fisher process, admits the following alternative representation

    p̃(·) =∑_h=1^H W_h δ_(ϕ_h, θ_h)(·),      (W_1,…,W_H)∼dirichlet(α/H,…,α/H),

where (ϕ_h,θ_h) ∼ P and with = denoting the equality in distribution. Thus, when σ < 0, the py reduces to a finite-dimensional discrete prior law having symmetric Dirichlet weights. One may refer to the Appendix A.1 of <cit.> for such a distributional equivalence. 

Although the dp is a special case of (<ref>), throughout the paper we will make extensive use of an alternative construction based on completely random measures. This approach is somewhat less straightforward compared to (<ref>), but it has important theoretical implications. Broadly speaking, a dp can be obtained as the normalization of a Gamma process <cit.>, which in turn can be represented as μ̃(·) = ∑_h=1^∞ J_h  δ_(ϕ̃_h,θ̃_h)(·), where J_1 ≥ J_2 ≥… is a collection of ordered positive random jumps whose distribution is given in <cit.>, that are independent on the random locations (ϕ̃_h,θ̃_h) ∼ P. The law of a Gamma random measure μ̃ is uniquely characterized by its Laplace functional, namely

    E{e^- μ̃(f)} = exp{ - α∫_X×Ylog{1 + f(x,y)}  P( x,  y) },

with α >0 and for any positive and measurable function f : X×Y→R^+ such that μ̃(f) = ∫_X×Yf(x,y) μ̃( x,  y) < ∞ almost surely. We write μ̃∼gap(α P).  The dp is then defined as the normalization of μ̃, that is, if μ̃∼gap(α P) then p̃(·) = μ̃(·) / μ̃(X×Y) ∼dp(α P).




 §.§ Definition and alternative representations

Informally, the epy process is obtained by integrating a collection of independent Pitman–Yor processes over a Dirichlet process. The definition has the same rationale that underlies the construction of enriched conjugate priors for natural exponential families <cit.>, which is extended to nonparametric settings in the edp and, to some extent and for univariate random measures, in neutral to the right processes <cit.>. 
Roughly speaking, a multivariate distribution is decomposed in terms of the marginal of and the conditional distributions. Then, an enriched prior law on the parameters of the joint distribution is obtained by assigning independent conjugate priors on the parameters of the marginal and on those of the conditional distributions. In the nonparametric case, the construction is more delicate as the distributions involved are random probability measures. 

 Let P be a probability measure on the product space X×Y, with P(A × B) = ∫_A P_Y | X(B | x) P_X( x) for any Borel sets A ⊆X and B ⊆Y. Moreover, let α > 0 and σ(x), β(x) be functions such that either σ(x) ∈ [0,1) and β(x) > - σ(x) or σ(x) = -β(x)/H(x) and β(x) > 0, where  H(x) ∈{2,3,…}. Define a random probability measure p̃_X on X and a family of random probability measures p̃_Y | X(·| x) on Y, for x ∈X, such that

    p̃_X ∼dp(α P_X),    p̃_Y | X(·| x) ∼py{σ(x), β(x) P_Y | X(·| x) },      x ∈X,

independently among themselves. Then the random probability measure p̃ on the product space X×Y, defined as

    p̃(A × B) = ∫_A p̃_Y | X(B | x) p̃_X( x),       A ⊆X,    B ⊆Y,

is said to be distributed as an enriched Pitman–Yor process (epy) with parameters α P_X, σ(x) and β(x) P_Y | X.  We will write p̃∼epy(α P_X, σ(x), β(x) P_Y | X).


Hence, the edp of <cit.> is a special case of the epy if σ(x) = 0 for any x ∈X. Ensuring that an epy process is a well defined stochastic process, namely that it provides a law of the random joint probability measure p̃, is a quite subtle measure-theoretic issue,  involving random conditional distributions. However, this can be shown along the steps given in <cit.> for the edp, which in turn are based on the results of <cit.>. Indeed, the only conditions there required are that the conditionals p̃_Y | X( ·| x) are independent across x ∈X, the marginal _X is a.s. discrete and _X and _Y | X are independent among themselves; all these properties hold for the epy by construction. A detailed proof is provided in the Supplementary Material for completeness.   


Note that the dp on the product space X×Y is a limiting case of the epy, occurring when β(x) → - σ(x) for any x ∈X. Indeed, at the limit, each conditional law reduces to a point mass, i.e. p̃_Y | X(B| x) = δ_θ_1(x)(B), with θ_1(x) ∼ P_Y | X(·;x), therefore 

    p̃(A × B) = ∑_h=1^∞ξ_h  δ_ϕ_h(A) δ_θ_1(ϕ_h)(B) = ∑_h=1^∞ξ_h  δ_(ϕ_h,θ_h)(A × B),      (ϕ_h,θ_h) ∼ P,

as each β(x) → -σ(x), where (ξ_h)_h ≥ 1 are the stick-breaking weights of p̃_X. Summarizing, an epy process p̃ on X×Y reduces to a dp(α P), as β(x) → - σ(x) for any x ∈X.  In addition, note that the baseline measures P_X, P_Y | X and P can be interpreted as “prior guesses” for the distribution of the observations, because one has

    E{p̃_X(A)} = P_X(A),   E{p̃_Y | X(B | x)} = P_Y | X(B | x),   E{p̃(A × B)} = P(A × B),

for any x ∈X and Borel sets A ⊆X, B ⊆Y, recalling that P(A × B) = ∫_A P_Y | X(B | x)P_X( x). 

The epy can be alternatively defined through a square-breaking representation. Such an equivalent definition, presented in the next proposition, is important especially for computational reasons, as one might truncate the involved series to approximate the infinite-dimensional process; see e.g. <cit.>. In addition, it emphasizes that the epy is a discrete random probability measure. 

 Let the quantities α P_X, σ(x), β(x) P_Y | X be as in Definition <ref> and let p̃ be a random probability measure on X×Y such that p̃∼epy(α P_X, σ(x), β(x) P_Y | X). Then

    p̃(A × B) =∑_ℓ =1^∞∑_h=1^∞ξ_ℓ π_h(ϕ_ℓ) δ_ϕ_ℓ(A) δ_θ_h(ϕ_ℓ)(B),      A⊆X,   B⊆Y,

where ξ_ℓ=ν_ℓ∏_j=1^ℓ-1(1-ν_j) and π_h(x)=η_h(x)∏_j=1^h-1{1-η_j(x)} for any ℓ≥ 1 and h ≥ 1, with

    ν_j    ∼beta(1,α),           ϕ_ℓ∼ P_X, 
    η_j(x)    ∼beta{1-σ(x),β(x) + jσ(x)},           θ_h(x) ∼ P_Y | X(·| x),

independently among themselves for any j ≥ 1, ℓ≥ 1, h ≥ 1 and x ∈X.


When σ(x) = 0, one recovers the square-breaking construction of the edp given in <cit.>. Moreover, if the discount parameter σ(x) is strictly negative for any x ∈X, then the square-breaking representation simplifies, because in this case the  conditional laws p̃_Y | X of the epy process are finite-dimensional; refer to Section <ref>.  

Paralleling the construction of the dp, we present a third alternative definition of the epy process, through the normalization of a random measure that we call Gamma and Pitman–Yor process (ga-py). 

 Let the quantities α P_X, σ(x), β(x) P_Y | X be as in Definition <ref>. Define a random probability measure μ̃_X on X and a family of random probability measures p̃_Y | X(·| x) on Y, for x ∈X, such that μ̃_X ∼gap(α P_X) and p̃_Y | X(·| x) ∼py{σ(x), β(x) P_Y | X(·| x) }, independently among themselves. Then the random probability measure μ̃ on the product space X×Y, defined as

    μ̃(A × B) = ∫_A p̃_Y | X(B | x) μ̃_X( x),       A ⊆X,    B ⊆Y,

is said to be distributed as a Gamma and Pitman–Yor process (ga-py) with parameters α P_X, σ(x) and β(x) P_Y | X.  We will write μ̃∼ga-py(α P_X, σ(x), β(x) P_Y | X).


Let μ̃ be a random measure with μ̃∼ga-py(α P_X,σ(x),  β(x) P_Y | X). Then, the random probability measure p̃ on the product space X×Y

    p̃(·) = μ̃(·)/μ̃(X×Y),

is distributed according to an enriched Pitman–Yor process (epy) with parameters α P_X, σ(x) and β(x) P_Y | X. The normalizing constant μ̃(X×Y) in the above definition is a positive random variable such that μ̃(X×Y)= ∫_Xp̃_Y | X(Y| x)μ̃_X( x) = μ̃_X(X) almost surely (a.s.). Hence, for any Borel sets A ⊆X and B ⊆Y, an epy process can be written as follows

    p̃(A × B)    = ∫_A p̃_Y | X(B| x) μ̃_X/μ̃_X(X)( x) = ∫_A p̃_Y | X(B| x) p̃_X( x),

where p̃_X(·) = μ̃_X(·)/μ̃_X(X) ∼dp(α P_X). 




§ PREDICTIVE RULE AND POSTERIOR LAW


We now turn to the investigation of the conditional properties of the epy. Consider the exchangeable random sequence ((X_n,Y_n))_n≥ 1 such that

    (X_n,Y_n) |p̃∼p̃,     p̃∼epy(α P_X, σ(x), β(x) P_Y | X),      n ≥ 1.

Note that we can sample each (X_n,Y_n) through a two-step mechanism. Specifically, X_n is first drawn from the marginal distribution p̃_X and then, given X_n = x, each Y_n is obtained from the conditional distribution p̃_Y | X(·| x); more precisely,

    X_n |p̃_X  ∼p̃_X,      Y_n | X_n = x,  p̃_Y | X(·| x) ∼p̃_Y | X(·| x),      n ≥ 1.

Recall that, if p̃ is an epy, the marginal law p̃_X is a dp(α P_X) whereas each conditional p̃_Y | X is a py(σ(x), β(x)P_Y | X), independently among themselves and on p̃_X for any x ∈X. Our results will easily follow from (<ref>) and well-known properties of dp and py processes. To facilitate their derivation, we shall assume that each conditional baseline measure P_Y | X(·| x) is a.s. diffuse, that is, it does not have a discrete component. In contrast, the marginal baseline measure P_X may have atoms or even being discrete. Indeed, in the py case the assumption of a diffuse baseline measure is essential to avoid more complicate probabilistic calculations  <cit.>. 




 §.§ Enriched urn scheme


We discuss an enriched urn scheme that extends the predictive mechanism of the edp given in <cit.>. By “urn scheme” we broadly mean an urn-like predictive rule, such as the one described in <cit.>. For the epy, the urn scheme provided in the next theorem sheds light on the underlying random partition and highlights the importance of the additional set of parameters σ(x), which indeed allows for a much finer calibration of the random partition compared to the edp.

The a.s. discreteness of the marginal law p̃_X implies that, in model (<ref>), there will be ties in the realization x^(n) = (x_1,…,x_n) of the random variables X^(n) = (X_1,…,X_n), with positive probability. Let (x_1^*,…,x_k_x^*) denote the k_x distinct values within x^(n), in order of appearance, with associated frequencies (n_1,…,n_k_x), so that n_1 + ⋯ + n_k_x = n. In the sampling mechanism (<ref>), the random variables (Y_1r,…,Y_n_rr) associated to the rth distinct value x_r^* are conditionally iid draws from the discrete distribution p̃_Y | X(·| x_r^*). Thus, with positive probability there will be further ties in each sample (Y_1r,…,Y_n_rr), with distinct values (y_1r^*,…,y_k_rr^*), in order of appearance, and with frequencies (n_1r,…,n_k_r r), so that n_1r + ⋯ + n_k_r r = n_r. Hence, the number of distinct values k_y = k_1 + ⋯ + k_k_x within a realization y^(n) = (y_1,…,y_n) of Y^(n) = (Y_1,…,Y_n) is such that k_x ≤ k_y.  This two-step stochastic mechanism can be described through an urn scheme, as illustrated by the next theorem. 


 Suppose ((X_n,Y_n))_n ≥ 1 is an exchangeable sequence as in equation (<ref>), with p̃∼epy (α P_X, σ(x), β(x) P_Y | X). Moreover, suppose that for each x ∈X the probability measure P_Y | X(·| x) is a.s. diffuse. Then, (X_1,Y_1) ∼ P and for any n ≥ 1

    X_n+1|X^(n) = x^(n)   ∼α/α + nP_X + 1/α+n∑_r=1^k_x n_r δ_x^*_r,

where (x_1^*,…,x_k_x^*) are the k_x distinct values within x^(n) with frequencies (n_1,…,n_k_x). Moreover,  for any r = 1,…,k_x and n ≥ 1

    Y_n+1|X^(n) = x^(n), Y^(n) = y^(n), X_n+1 = x 
    
           ∼β_r + k_r σ_r/β_r + n_rP_Y | X(·| x_r^*) + 1/β_r + n_r∑_j=1^k_r(n_jr - σ_r)  δ_y_jr^*,        x = x_r^*,
           ∼ P_Y | X(·| x),        x ∉{x_1^*,…,x_k_x^*}
    ,

with β_r = β(x_r^*) and σ_r = σ(x_r^*), where (y_1r,…,y_n_rr) are the values of y^(n) associated to x_r^*, whereas (y_1r^*,…,y_k_rr^*) are the corresponding k_r distinct values, with frequencies (n_1r,…,n_k_r r), for r = 1,…,k_x. 


 The system of predictive laws in the above theorem uniquely characterizes the epy process; refer to the Supplementary Material for details. 


The two-stage random partition implied by the predictive rule can be described in terms of a nested Chinese restaurant (crp) metaphor, illustrated in Figure <ref>, which extends the one of <cit.>. Let us assume that P is diffuse almost surely. Suppose there exists a restaurant with a potentially infinite number of tables, representing the X_n, which serves a potentially infinite number of dishes, representing the Y_n. A first customer seats in one of the tables and selects a dish. For n ≥ 1, the (n+1)th customer may either sit in one of the occupied  tables, say the rth, with probability n_r / (α + n) for r=1,…,k_x, or she can seat in a new one with probability α/(α + n). If a new table is chosen, she will get a new dish. Otherwise, she may either select a new dish with probability (β_r + k_r σ_r) / (β_r + n_r) or choose one of dishes previously served at her table, say the jth, with probability (n_jr - σ_r)/(β_r + n_r), for j=1,…,k_r.  In comparison, the classical Chinese restaurant process only partitions customers in tables, disregarding the dish of choice. Such a nested crp describes the random partition implied by the epy. Moreover, if we label the tables with iid draws from P_X, and the dishes at the xth table with iid draws from P_Y | X(·| x), then one obtains the enriched Pólya sequence ((X_n, Y_n))_n ≥ 1 defined by the predictive rule (<ref>)-(<ref>). 





If σ(x) = 0 for any x ∈X, then the sequence ((X_n, Y_n))_n ≥ 1 with predictive rule as in Theorem <ref> corresponds to the enriched Pólya sequence described in <cit.>. Under the hypothesis of Corollary <ref>, ((X_n, Y_n))_n ≥ 1 is a Pólya sequence <cit.>.  Moreover, this also occurs when β(x) → -σ(x) for any x ∈X, consistently with the discussion in Section <ref>. As in the edp, the precision parameter α and the function β(x) regulate the number of distinct values within X^(n) and Y^(n). However, the additional function σ(x) controls the asymptotic clustering behavior of each subsequence (Y_1r,…,Y_n_rr), given X^(n) = x^(n), that is, the growth rate of the number of clusters of (Y_1r,…,Y_n_rr). In addition, the discount parameter σ(x) allows to regulate the variance of the within-group number of clusters k_r, leading to more robust specifications. It is hence clear that the epy allows for a much  greater flexibility compared to the edp. We refer to <cit.> and <cit.> for an extensive discussion about the role of the discount parameter σ(x) and its usefulness both for species sampling and mixture models.

For example, positive values of σ(x) ∈ (0,1) lead to a within-group polynomial growth rate of the number of distinct values, which is much faster than the logarithmic rate of the dp, occurring when σ(x) = 0. Conversely, if the discount parameter is negative, i.e. σ(x) = - β(x)/ H(x)  with β(x) >0, then number of clusters is bounded by H(x). Indeed, in this case for any n ≥ 1

    Y_n+1|X^(n) = x^(n), Y^(n) = y^(n), X_n+1 = x 
           ∼(1 - k_r/H_r)β_r/β_r + n_rP_Y | X(·| x_r^*) + 1/β_r + n_r∑_j=1^k_r(n_jr + β_r / H_r)  δ_y_jr^*,     x = x_r^*,

with β_r = β(x_r^*), σ_r = σ(x_r^*) and H_r = H(x_r^*). The above equation highlights that the within-group number of clusters cannot be greater than H_r, a feature which has been provably useful in several applied contexts, including the application of Section <ref>.  




 §.§ Posterior distribution


We now derive the posterior law of the random probability measure p̃. The epy process is not conjugate but the corresponding posterior is nonetheless analytically tractable. Recall that, by definition, (A × B) = ∫_A p̃_Y | X(B| x) p̃_X( x), for A ⊆X  and B⊆Y. Therefore, its posterior distribution may be obtained, at least in principle, from the posterior laws of the random marginal distribution p̃_X  and of the conditionals p̃_Y | X(·| x). Those posterior distributions are obtained in the following theorem.

 Let (X_i,Y_i) |p̃∼p̃ for i=1,…,n, with p̃∼epy(α P_X, σ(x), β(x) P_Y | X), and suppose that each conditional probability measures P_Y | X(·| x) is a.s. diffuse. Then, under the notation of Theorem <ref>, one has

    p̃_X |X^(n) = x^(n),Y^(n) = y^(n)∼dp(α P_X + ∑_r=1^k_x n_r δ_x^*_r).

Moreover, for any x ∈X one has

    p̃_Y | X(·| x) |X^(n) = x^(n), Y^(n)= y^(n)   = W_0r p̃_Y | X^*(·| x_r^*) + ∑_j=1^k_r W_jr δ_y_jr^*(·),   x = x_r^*, 
    p̃_Y | X(·| x) |X^(n) = x^(n), Y^(n)= y^(n)   ∼py{σ(x), β(x) P_Y | X(·| x)},  x ∉{x_1^*,…,x_k_x^*},

independently on p̃_X and among themselves, where 

    (W_0r,W_1r,…,W_k_rr) ∼dirichlet( β_r + k_rσ_r, n_1r - σ_r,…, n_k_r r - σ_r)
 
and for any r =1,…,k_x

    p̃_Y | X^*(·| x_r^*) ∼py{σ_r, (β_r + k_rσ_r)P_Y | X(·| x_r^*)}.





Note that if σ_r = - β_r/ H_r is strictly negative in the above theorem, then the random probability measure p̃_Y | X^*(·| x_r^*) follows a Dirichlet multinomial process with H_r - k_r components, that is

    p̃_Y | X^*(·| x_r^*) =∑_j=k_r+1^H_rW_jr^* δ_θ_j r(·),

where (W_k_r+1 r^*,…,W_Hr^*) ∼dirichlet(β_r/H_r,…,β_r/H_r) and θ_jr∼ P_Y | X(·| x_r^*) for any j=k_r + 1,…,H. Hence, the posterior law of p̃_Y | X is finite-dimensional, meaning that it is characterized by a finite number of random variables. On the other hand, if we set σ(x) = 0 for all x ∈X, we obtain that, if x = x_r^*,

    p̃_Y | X(·| x) | (X^(n) = x^(n), Y^(n) = y^(n)) ∼dp(β_r P_Y | X(·| x_r^*) + ∑_j=1^k_rn_jr δ_y_jr^*(·) ),

 which implies that p̃| (X^(n) = x^(n), Y^(n) = y^(n)) is an edp with updated parameters, as established in <cit.>. 







§ THE MARGINAL PROCESS WITH A DISCRETE BASELINE MEASURE


Recall from equations (<ref>)-(<ref>) that an exchangeable sequence (Y_n)_n ≥ 1 is directed by a marginal epy process if

    Y_n |p̃_Y    ∼p̃_Y,      n ≥ 1,

where p̃_Y(·) = p̃(X×·) and p̃∼epy(α P_X,σ(x),  β(x) P_Y | X). In this section we focus on a special case of the epy process, arising when the marginal baseline measure P_X is discrete, that is when X_n takes values on a fixed set of points. More precisely, let P_X(·) = ∑_ℓ=1^L α_ℓ/α δ_x_ℓ(·), with x_1,…,x_L ∈X and α = ∑_ℓ = 1^L α_ℓ. For the sake of the exposition we consider L < ∞, although our results may be easily extended to the countable case. Thus, a random probability measure p̃_Y on Y follows a marginal epy process with a discrete baseline measure P_X if

    p̃_Y(·)    =∑_ℓ=1^L Π_ℓ p̃_ℓ(·) = ∑_ℓ=1^L ∑_h=1^∞Π_ℓ π_ℓ hδ_θ_ℓ h(·), 
       (Π_1,…,Π_L)    ∼dirichlet(α_1,…,α_L),     p̃_ℓ∼py(σ_ℓ, β_ℓ P_ℓ),

having set p̃_ℓ(·) = p̃_Y | X(·| x_ℓ), θ_ℓ h = θ_h(x_ℓ), π_ℓ h = π_h(x_ℓ), σ_ℓ = σ(x_ℓ), β_ℓ = β(x_ℓ), and P_ℓ(·) = P_Y | X(·| x_ℓ), for notational convenience. Compared to the general square-breaking representation, in the above equation the X-valued random locations (ϕ_h)_h ≥ 1 are replaced by the fixed values x_1,…,x_L. 

In the following, we discuss further theoretical properties and we clarify the link between the prior (<ref>) and other methods available in the literature. Interestingly, although seemingly unrelated, these proposals arise as special cases of the marginal epy process with a discrete P_X. This highlights the central role of the epy in a variety of contexts. In addition, making such a connection explicit allows us to develop extensions and to obtain novel modeling strategies and computational advances, that naturally arise in our unifying framework.




 §.§ Theoretical characterizations


 Discrete baseline measures may have important and unexpected distributional consequences <cit.>. In our case, such an assumption for P_X leads to remarkable simplifications. To illustrate the effects of this choice, let us consider μ̃_X ∼gap(α P_X) with P_X(·) = ∑_ℓ=1^L α_ℓ/α δ_x_ℓ(·).  Then, it holds

    μ̃_X(·) =∑_ℓ=1^L V_ℓ δ_x_ℓ(·),     V_ℓ∼ga(α_ℓ, 1),    ℓ=1,…,L,

which is arguably a simpler representation than the general <cit.> series discussed in Section <ref>. Indeed, the locations are  deterministic and the jumps are independent.  As a consequence, the distribution of a marginal epy process _Y is in the form of (<ref>).

 We now discuss characterization theorems for the joint epy process  with a discrete P_X, which may be used to study its distributional properties. Indeed, in this case, the Laplace functional characterizing a Gamma and Pitman–Yor random measure admits a simple expression, highlighting important connections with Cauchy-Stieltjes transforms. This is clarified in the following theorem. 

 Let μ̃∼ga-py(α P_X,σ(x),  β(x) P_Y | X) with P_X(·) = ∑_ℓ=1^L α_ℓ/ α δ_x_ℓ(·). Then,

    E{e^-μ̃(f)} = ∏_ℓ=1^L E[{1 + p̃_Y | X(f| x_ℓ) }^-α_ℓ],

where p̃_Y | X(f| x) = ∫_Yf(x,y)p̃_Y | X( y| x) for any x ∈X and for any positive and measurable function f : X×Y→R^+ such that μ̃(f) < ∞ almost surely.


The expectation appearing in the right hand side of the above Laplace functional is called generalized Cauchy-Stieltjes transform and it can be computed in closed form in some special cases. For example, the so-called Cifarelli-Regazzini identity <cit.> implies that if p̃∼dp(α P) and μ̃∼gap(α P) then

    E[{1 + p̃(f)}^-α] = E{e^-μ̃(f)},

for any positive and measurable function f : X×Y→R^+ such that μ̃(f) < ∞ almost surely. 
As an application of the identity in equation (<ref>), we obtain the next Corollary.

 Let μ̃∼ga-py(α P_X, σ(x), β(x) P_Y | X) and let P_X(·) = ∑_ℓ=1^L α_ℓ/ α δ_x_ℓ(·). Moreover, assume that σ(x_ℓ) = 0 and α_ℓ = β(x_ℓ) for any ℓ=1,…,L. Then

    E{e^-μ̃(f)} = exp{ - α∫_X×Ylog{1 + f(x,y)}  P( x,  y) },

for any positive and measurable function f : X×Y→R^+ such that μ̃(f)  < ∞ almost surely.


Hence, under the hypotheses of Corollary <ref>, a ga-py random measure reduces to a Gamma process. In turn, this implies that an edp with discrete baseline measure P_X(·) = ∑_ℓ=1^L α_ℓ / α δ_x_ℓ(·) and whose parameters satisfy the constraint α_ℓ = β(x_ℓ) for any ℓ=1,…,L, is distributed as a dp(α P).  A similar consideration was made by <cit.>, who obtained this result by inspecting the predictive distributions. Instead, our proof relies on the ga-py process. Thus, again under the assumption of Corollary <ref>, the marginal epy p̃_Y of equation (<ref>) becomes a dp, namely 

    p̃_Y(·) =  ∑_ℓ=1^L Π_ℓ p̃_Y | X(·| x_ℓ) ∼dp(∑_ℓ=1^L α_ℓ P_Y | X(·| x_ℓ)).

Therefore, the baseline measure ∑_ℓ=1^L α_ℓ/ α  P_Y | X(·| x_ℓ) is as a mixture.

The equivalent of the Cifarelli-Regazzini identity with p̃∼py(σ,α P), has been obtained by <cit.> for positive σ∈ (0,1) and α > 0. This leads to a second specialization of Theorem <ref>, which is summarized in the following Corollary. 

  Let μ̃∼ga-py(α P_X, σ(x), β(x) P_Y | X) and let P_X(·) = ∑_ℓ=1^L α_ℓ / α δ_x_ℓ(·). Moreover, assume that σ(x_ℓ) > 0 and α_ℓ = β(x_ℓ) for any ℓ=1,…,L. Then

    E{e^-μ̃(f)} = ∏_ℓ=1^L [∫_Y{1 + f(x_ℓ,y)}^σ(x_ℓ) P_Y | X( y| x_ℓ) ]^-β(x_ℓ)/σ(x_ℓ),

for any positive and measurable function f : X×Y→R^+ such that μ̃(f)  < ∞ almost surely. 


Corollary (<ref>) has its own theoretical interests, as it uniquely characterizes a specific ga-py random measure. In addition, it implies that the equivalent of equation (<ref>) does not hold true in the py case. Specifically, under the hypothesis of Corollary <ref>, the marginal epy process p̃_Y(·) =  ∑_ℓ=1^L Π_ℓ p̃_Y | X(·| x_ℓ), is in general not distributed as a py. This may still occur in a very special case, i.e. when the baseline measures and the discount parameters are equal, namely P_Y | X(·| x_1) = ⋯ = P_Y | X(·| x_L) and σ(x_1) = ⋯ = σ(x_L) > 0, which is a known result; see e.g. Proposition 14.35 in <cit.>. 

Few other generalizations of the Cifarelli-Regazzini identity are known. Specifically, consider the transform E[{1 + p̃(f)}^-α], with p̃∼dp(β P), for some positive β≠α. If such a transform were available, this would allow us to obtain the equivalent of Corollary <ref> without imposing the constraint α_ℓ = β(x_ℓ). <cit.> established analytic results for this transform, while <cit.> provided a probabilistic interpretation of such a transform in terms of the Laplace functional of a Beta-Gamma process. However, these findings require more analytical efforts compared to (<ref>).  




 §.§ Discrete priors with atomic contaminations


dp priors with a single atomic component are widely used in Bayesian nonparametrics. 

The common thread is the employment of a discrete random measure p̃_Y on Y having the form

    p̃_Y(·) = Π δ_y_0(·) + (1 - Π) p̃_2(·),     p̃_2  ∼dp(β_2 P_2),

with y_0 ∈Y being a fixed atom, P_2 a diffuse probability measure, and with Π∼beta(α_1,α_2), independently on p̃_2. In some cases <cit.>, we have that y_0 = 0 and therefore p̃_Y may be called a spike and slab dp prior. In contrast, in <cit.> the atom y_0 is allowed to be random. Under the additional constraint α_2 = β_2, the self-similarity property of the dp implies that model (<ref>) reduces to

    p̃_Y ∼dp(α_1 δ_y_0 + α_2 P_2),

which is the specification described in <cit.>, with y_0 = 0. Therefore, models (<ref>) and (<ref>) are closely related and they are sometimes called “outer” and “inner” spike and slab, respectively.

As it may be already clear, model (<ref>) is a marginal epy process with discrete P_X, because it is in the form of equation (<ref>). The equivalence between models (<ref>)-(<ref>) can be viewed as a consequence of Corollary <ref>. This can be better appreciated by noticing that the point mass δ_y_0 is also a trivial dp, namely δ_y_0 = p̃_1 ∼dp(α_1 P_1), with baseline measure P_1 = δ_y_0.  In the epy nested clustering mechanism, the random variables X_n should be interpreted as latent quantities that can only take L = 2 values, and the underlying baseline measure is P_X = α_1 / α δ_x_1 + α_2 / α δ_x_2. In particular, each X_n identifies whether the corresponding Y_n has to be sampled from the  atomic contamination δ_y_0 (i.e. X_n = x_1), or from the  nonparametric component p̃_2 (i.e. X_n = x_2).  




This link between model (<ref>) and epy processes leads to the following natural extension. Indeed, a simple generalization of  (<ref>) accounting for a more flexible clustering mechanism is

    p̃_Y(·) = Π δ_y_0(·) + (1 - Π) p̃_2(·),     p̃_2  ∼py(σ_2, β_2 P_2).

Hence, the resulting p̃_Y is still a marginal epy process and therefore it remains analytically tractable. Motivated by similar considerations, <cit.> studied a py process p̅_y having a contaminated baseline measure, namely p̅_y ∼py(σ, α_1 δ_y_0 + α_2 P_2). Importantly, Corollary <ref> implies that p̅_y cannot be regarded as an epy process, and therefore the equivalence between inner and outer models, as in equations (<ref>)-(<ref>), is lost beyond the dp special case. Thus, the marginal epy process of equation (<ref>) and the prior of <cit.> are different in general, although they may be regarded as closely related  alternatives for the modeling of atomic contaminations. The main distinction is on the computational side: posterior results for p̅_y of <cit.> can be derived, but their practical usage may be complicated due to the presence of cumbersome combinatorial quantities. Conversely, the posterior law of p̃_Y can be readily obtained as a straightforward modification of Theorem <ref>. The spike and slab specification in equation (<ref>) has been independently proposed and carefully investigated by <cit.>. 




 §.§ Mixture of mixtures models




Bayesian nonparametric discrete priors are commonly employed for mixture modeling. In such a setting, the random variables Z_1,…,Z_n are iid draws from a random density f̃, such that

    Z_i |f̃∼f̃,     f̃(z) = ∫_Y𝒦(z | y)p̃_Y( y),      i=1,…,n,

where 𝒦(z | y) is a kernel density and where p̃_Y is a discrete random probability measure. For example, scale-location mixtures of Gaussian kernels are a  popular special case of (<ref>) when Z_i is a vector on R^p. If p̃_Y follows the marginal epy process in equation (<ref>), then model (<ref>) may be termed a “mixture of mixtures” <cit.>, because the random density f̃ becomes 

    f̃(z)=∑_ℓ=1^L Π_ℓ∫_Y𝒦(z | y)p̃_ℓ( y) =  ∑_ℓ=1^L ∑_h=1^∞Π_ℓ π_ℓ h𝒦(z |θ_ℓ h),

which is, indeed, a mixture model whose kernel ∫_Y𝒦(z | y)p̃_ℓ( y) is itself a mixture. Consistently with the nested partition mechanism described in Section <ref>, there will be two levels of clustering, regulated by the latent variables X_1,…,X_n and Y_1,…,Y_n associated to the epy process. The variables X_1,…,X_n control the global clustering, that is, they identify which kernel ∫_Y𝒦(z | y)p̃_ℓ( y) has to be considered. Then, conditionally on the X_1,…,X_n, the variables Y_1,…,Y_n regulate the local clustering occurring within each kernel mixture ∫_Y𝒦(z | y)p̃_ℓ( y). Note that this mechanism is likely to be affected by severe identifiability issues, which may be mitigated by carefully specifying the baseline measures P_ℓ(·) = P_Y | X(·| x_ℓ). 

Motivated by applications to the study of early-birth risk in a population of women, 
<cit.> introduced model (<ref>) to cluster the Y_1,…,Y_n subject-specific latent functions as well as predicting the Y_n+1th curve associated to a new subject. The models of <cit.> and <cit.> are also special cases of (<ref>), when the discount function of the underlying epy process σ(x) is strictly negative for all x ∈X. 












 §.§ Dependent enriched processes


If the random variables Y_i^(j) are structured into groups, for unit i = 1,…,n^(j) and group j=1,…,d, then the conditional iid assumption of model (<ref>) may be inappropriate. In fact, one may assume exchangeability only within the same group, that is

    Y_n^(j)|p̃_j,Y   iid∼p̃_j,Y,      n ≥ 1,

for j=1,…,d, with the vector of random probability measure (p̃_1,Y,…,p̃_d,Y) following some prior law 𝒬_d. Clearly, if 𝒬_d models independence across the p̃_j,Y, then the (Y_n^(j))_n ≥ 1 would reduce to independent exchangeable sequences. 
However, this may cause an important loss of efficiency, failing in borrowing strength across groups.  Hence, it is important to specify a prior law 𝒬_d that induces dependence across the random probability measures p̃_j,Y. Among the available proposals, we focus on a special case of <cit.>. Define

    p̃_j,Y(·)    =  ω_j q̃_j(·) + (1 - ω_j)q̃_0(·),            ω_j ∼beta(α_1, α_2), 
    q̃_0    ∼dp(α_2 P_Y),           q̃_j ∼dp(α_1 P_Y),

independently for j=1,…,d, for some α_1, α_2 > 0, and baseline measure P_Y. This model is closely related to the approach of <cit.>, who assumed a mixture distribution for ω_j with point masses at 0 and 1, and, in addition, that ω_1 = ⋯ = ω_d. Specification (<ref>) induces dependence across groups through the presence of a common random probability measure q̃_0. Moreover, the self-similarity property of the dp implies that the marginals p̃_j,Y are themselves dps, namely p̃_j,Y∼dp{(α_1 + α_2) P_Y}.

In fact, model (<ref>) is of the form (<ref>), thus implying that each p̃_j,Y is distributed as a marginal edp. Hence, extensions leveraging general epy processes can be envisioned. In first place, we could set 

    q̃_0 ∼dp(β_2 P_Y),     q̃_j ∼dp(β_1 P_Y),      j=1,…,d,

with α_1 ≠β_1 and α_2 ≠β_2, therefore allowing for a richer parametrization. This simple modification has quite useful implications. Indeed, it allows two random probability measures p̃_j,Y to be highly correlated (i.e. α_2 →∞) while having few number of clusters (i.e. β_2 ≈ 0). This is not possible in the framework of <cit.>, where high dependence among the p̃_j,Y is necessarily associated with a larger number of clusters. Even more generally, we could let 

    q̃_0 ∼py(σ_2, β_2 P_Y),     q̃_j ∼py(σ_1, β_1 P_Y),      j=1,…,d.

This still implies that p̃_j,Y are marginal epys, and it allows the number of clusters to have different growth rates, as discussed in Section <ref>. 



§ APPLICATIONS


We consider here two different case studies involving the epy process, devoting special emphasis to the elicitation or estimation of the hyperparameters α, σ(x), β(x) and P.  Indeed, as for any discrete nonparametric prior, the effect of the hyperparameters depends on the role of the random probability measure p̃ within the statistical model. For example, in the first case study p̃ is the conditional law of the observable random vectors (X_n, Y_n) and interest is in predicting the number of species within a future sample. In such a setting, the choice of the baseline measure P is irrelevant as it only plays the role of labeling the different species.  In addition, we assume the diffuseness of both P_Y | X and P_X and we make direct usage of the predictive rule obtained in Section <ref>. Instead, in the second case study we consider a mixture of mixtures model described in Section <ref>, where interest is in  clustering subject-specific functions according to local  Y_1,…,Y_n and global X_1,…,X_n latent features. In this case, the choices of the conditional baseline measures P_Y | X and of the discrete marginal measure P_X have a relevant impact on the inferential results, as it commonly occurs in mixture models. 



 §.§ The Amazonian tree flora dataset


The Amazonian flora is the richest assemblage of plant species on Earth. However, the exact number of tree species present in the Amazon basin, or an estimate of the proportion between common and rare species, is still unavailable. The lack of this basic information means that ecologists do not have a clear picture of the world's largest trees community.  Here we analyze the same dataset of <cit.>, which is openly available online. 
Specifically, we aim at predicting the number of new species that researchers are expected to encounter in subsequent surveys, a problem with a very rich statistical literature; see e.g. <cit.>. Refer also to <cit.> for a historical account on the topic. The estimation of the unobserved number of species has been more recently addressed using Bayesian nonparametrics tools, most notably in <cit.> and <cit.>, who spurred an interesting piece of literature focusing on species sampling models beyond the Dirichlet process. However, these approaches focus on the marginal distribution of the species, disregarding the relevant information provided, for instance, by the corresponding families. We will empirically show that enriched specifications outperform marginal models in these specific applications as they can leverage a richer parametrization and a greater amount of information.

A total of n = 553949 trees has been recorded in our dataset, comprising k_y = 4962 different species and k_x = 115 families of trees. The data consists in a collection of frequencies n_jr denoting how many times the jth species of the rth family has been observed. For example, the Euterpe oleracea species belongs to the Arecaceae family and it has been observed 8572 times. We let n_r and  k_r be the number of trees and the number of distinct species 
associated to the rth family, so that ∑_r=1^k_x n_r = ∑_r=1^k_x∑_j=1^k_r n_jr = n. For instance, the Arecaceae family comprises 70 different species for a total of 51862 trees. We aim at providing an accurate estimate of the number of new species k_y(m) ∈{0, 1, …,m} that one would observe if an additional sample of m ≥ 1 trees were collected. The epy model constitutes a natural probabilistic framework for this problem. Let X_i and Y_i be the random variables denoting the family and the species of the ith tree in the sample, respectively. We assume

    (X_i, Y_i) |p̃   ∼p̃,      i ≥ 1,
    p̃   ∼epy(α P_X,σ(x),  β(x) P_Y | X).

The distributions P_X and P_Y | X are diffuse probability measures. However, they only serve as a mathematical tool for identifying “new” or “old” families and species in the urn scheme of Theorem <ref>. Therefore, the baseline measure P can be treated as a nuisance parameter. 
Note that the diffuseness of P_X and P_Y | X does not imply that observations (X_i,Y_i) will be distinct; it only implies that new species are given new labels, as desirable. 
Still, the almost sure discreteness of p̃ will lead to ties among the data with positive probability. 

Recall that we are interested in the posterior distribution of k_y(m) given the data, namely

    k_y(m) |X^(n) = x^(n), Y^(n) = y^(n),      m ≥ 1,

which corresponds to the number of distinct species that one would observe within a future sample Y_n+1,…, Y_n+m and that were not observed among the data Y_1,…,Y_n. Similarly, albeit this is not the focus of this case study, one could be interested in estimating k_x(m), which is the number of new families that were not observed among X_1,…,X_n. The enriched Pólya urn scheme of Theorem <ref> provides a very efficient way to simulate independent samples from the posterior law of (<ref>) by first drawing samples for (X_n+1,Y_n+1),…,(X_n+m,Y_n+m), given the data, and then counting the distinct values among the simulated species Y_n+1,…, Y_n+m that were not previously observed. Beside being useful for simulations, Theorem <ref> and the associated nested Chinese restaurant metaphor shed light on the data generating process. Indeed, the families X_1,…,X_n are generated according to the <cit.> Pólya-urn scheme and then, inside each family, the species are obtained from the urn scheme of a Pitman–Yor prior. 





To perform predictions about k_y(m), we need to either specify or estimate the parameters α, σ(x) and β(x). For any fixed value x ∈X, the parameters σ(x) and β(x) are related to the growth rate and the number of species of that specific family. Since there is no natural ordering among families of trees, we let these values be themselves exchangeable. In particular, we assume that (σ(x), β(x)) 𝒬_σ,β, that is they are independent and identically distributed  – hence, a fortiori, exchangeable – according to some prior law 𝒬_σ,β. In principle, one could rely on a formal Bayesian procedure to obtain the posterior laws of α and each σ(x), β(x). In practice, this would require the implementation of a suitable Markov Chain Monte Carlo algorithm, which would heavily increase the computational burden. 
To simplify computations and following common practice in species sampling applications <cit.>, we rely on a plug-in estimate for those parameters, and in particular we use their posterior mode. In other words, in what is referred to as an empirical Bayes approach in these settings, we plug-in the values α̂, σ̂_r = σ̂(x^*_r) and β̂_r = β̂(x^*_r) maximizing a penalized likelihood. This computationally simpler strategy might be regarded as a variational Bayes approximation of the Bayesian posterior distribution from the prior law  𝒬_σ,β, in the lines envisaged in <cit.>. Proper Bayesian inference would be conceptually straightforward to implement, but the resulting Markov Chain Monte Carlo algorithm would be computationally more demanding than the proposed approximation. 

Let ℒ(x^(n), y^(n)|α, β_1,…,β_k_x,σ_1,…,σ_k_x) be the likelihood function associated to the epy process, written ℒ for notational simplicity. As a consequence of Theorem <ref>, when σ(x) > 0 one can show that

    ℒ∝   [ α^k_x/(α)_n∏_r=1^k_x(n_r - 1)! ] [ ∏_r=1^k_x∏_j=1^k_r-1 (β_r + jσ_r)/(β_r + 1)_n_r-1∏_j=1^k_r(1-σ_r)_n_jr-1],

where (a)_n = a(a+1)⋯(a+n-1) and (a)_0=1 is the Pochammer symbol. Moreover, let f(β_r,σ_r) = f(σ_r)f(β_r) be the densities, for each (β_r,σ_r), that are associated to the prior law Q_σ, β. In addition, let us denote with f(α) the prior density of α. Thus α̂ and the pairs (β̂_r,σ̂_r) for r=1,…,k_x are obtained as follows

    α̂   = max_α f(α) α^k_x/(α)_n∏_r=1^k_x(n_r - 1)!, 
     
    (σ̂_r, β̂_r)    = max_(σ_r,β_r) f(σ_r,β_r) ∏_j=1^k_r-1 (β_r + jσ_r)/(β_r + 1)_n_r-1∏_j=1^k_r(1-σ_r)_n_jr-1,

for r = 1,…,k_x. Hence, the posterior modes α̂, σ̂_r and β̂_r can be easily found via numerical maximization.  As for the prior distributions f(α), f(β_r) and f(σ_r), we let α∼gamma(2,0.01), σ_r ∼beta(1,10) and each β_r ∼gamma(2,1). These prior laws regularize the otherwise ill-behaved estimates occurring when the  functions in equation (<ref>) are unbounded. The hyperparameters are selected so that the estimates α̂, σ̂_r and β̂_r are comparable with the maximum likelihood estimates in the regular cases; see the Supplementary Material for details. These prior choices implicitly assume σ(x) > 0 for any x ∈X, implying that the number of species within each family is allowed to grow indefinitely with the sample size.  Finally, if a new family X^*_k_x+1 is drawn among the sampled X_n+1,…,X_n+m, we then set the corresponding parameters according to the prior information, namely (σ̂_k_x+1, β̂_k_x+1) = max_(σ_k_x+1,β_k_x+1) f(σ_k_x+1)f(β_k_x+1) = (0,1). 

We validate the predictive performance of the proposed model by comparing the epy with a marginal dp on the species Y_i and with an edp. We also considered a py specification; however, the obtained empirical Bayes estimate was σ≈ 0, therefore collapsing to a dp model.  We randomly split the dataset in a training set and a test set having n_ = 250000 and  n_ = n - n_ = 303949 observations, respectively. The hyperparameters of the dp and edp competing models were also estimated from the data; see the Supplementary Material for the details. Conditionally on the training set, we predict the number of unobserved species k_y(m), under the dp, edp and epy models for various choices of m=1,…,n_. We compare the predictions with the actual number of new species present in the full test set, if m = n_, or in a random subsample of it, when m < n_.  As apparent from Figure <ref>, all the competing methods provide a reasonable estimate for the number of unobserved species. However, the additional flexibility and information (i.e. the families) available to the enriched processes leads to more accurate predictions compared to the dp, with the epy being a slight improvement also with respect to the edp. 

Finally, in Figure <ref> we provide the predicted number of new species k_y(m) |X^(n)= x^(n),Y^(n)= y^(n) within m=1,…, 200000 subsequent samples, by employing an epy process and conditioning on the full dataset. The expected value and  90% pointwise credible intervals are provided.  We should remark that in principle one could evaluate k_y(m) for even larger values of m in the attempt of estimating the total number of species present in Amazonia. However, this might lead to misleading inferential conclusions, for example because it heavily relies on the homogeneity assumptions implicit in an exchangeable model. 



 §.§ E-commerce market segmentation


An e-commerce company is interested in understanding the preferences of its customers, to implement targeted marketing strategies. This company operates online and sells flight tickets. Our analysis is based on a dataset that contains the number of times a route has been searched on the company's website, comprising a collection of weekly counts for each flight route. A few examples of these functional observations are depicted in Figure <ref>. The main goal of the analyses is finding a limited number of route groups to target cluster-specific policies. Several Bayesian nonparametric approaches involving functional clustering have been recently proposed, including the works of <cit.>. We here revisit the work of <cit.> under the lenses of general epy processes, which allows us
to clarify several aspects related to the elicitation of its hyperparameters. 


Let us regard the route-specific measurements Z_i(t), for each route i=1,…,n and time value t ∈R^+, as error-prone realizations of unknown functions Y_i(t), that is

    Z_i(t)  = Y_i(t) + ϵ_i(t),      i=1,…,n,

with ϵ_i(t) denoting a Gaussian random noise term, independent over flight routes and time. 
We will write Z_i = {Z_i(t)}_t ∈R^+ and Y_i = {Y_i(t)}_t ∈R^+ for i=1,…,n to denote each functional observation and latent trajectory, respectively. To perform Bayesian clustering of the Z_1,…,Z_n, we let the latent trajectories Y_1,…,Y_n to be conditionally independent draws from the marginal epy process described in equation (<ref>), namely

    Y_i |p̃_Y iid∼p̃_Y,      i=1,…,n,

where p̃_Y(·) = p̃(X×·), p̃∼epy(α P_X,σ(x),  β(x) P_Y | X). Moreover, recall that we are assuming a discrete baseline measure P_X. This is an instance of a mixture of mixtures model, described in equation (<ref>), for the data Z_1,…,Z_n, represented in a hierarchical manner. Note that the atoms θ_ℓ h = θ_h(x_ℓ) = {θ_ℓ h(t)}_t ∈R^+ of p̃_Y are themselves functions and can be interpreted as the latent trajectories of each cluster. In this specific application, a priori information is available about the shapes of these curves. For example, some flight routes are characterized by a strong cyclical component, e.g. the ones depicted in Figure <ref>. The epy allows to naturally incorporate this prior knowledge into the model by assuming the existence of x_1,…,x_L latent feature classes, each describing a specific functional shape (i.e. monotone, cyclical, etc.). The allocation of each function Y_i to the corresponding class is regulated by the random variables X_1,…,X_n. For example, if X_i = x_ℓ then the ith latent trajectory Y_i belongs to the ℓth feature class characterizing, say, cyclical functions. 

Compared to the first case study, the elicitation of the hyperparameters α, σ(x), β(x) and of the baseline measures P_X, P_Y | X is more delicate.  In first place, note that each conditional base measure P_ℓ(·) = P_Y| X(·| x_ℓ) can be regarded as a “functional prior guess”, since for any t ∈R^+ we have that E{Y_i} =   ∑_ℓ=1^L α_ℓ / αE{θ_ℓ 1}, with θ_ℓ 1∼ P_ℓ,  for ℓ = 1,…,L. Note that if inference about the latent classes X_1,…,X_n is of interest, the measures P_ℓ must be “distinguishable”, i.e. they should characterize quite different functional shapes. Otherwise, it might be difficult to infer the global clusters from the data due to severe identifiability issues. In practice, the elicitation of the conditional base measures for functional data may be based on available prior information <cit.>. Consistent with the above discussion, we assume that each θ_ℓ h(t) is linear in the parameters, with a Gaussian prior on the regression coefficients. Moreover, we let L = 2 and we select the conditional base measures P_1 and P_2 so that they have interpretable shapes. The first functional class (ℓ=1) captures yearly cyclical patterns and characterizes the routes having a single peak of web-searches per year. This is the case for example of the Milan–Alghero and Milan–Cagliari routes, as apparent from Figure <ref>. Hence, we let

    θ_1h(t) = ∑_m=1^4 γ_m1h𝒮_m(t) +  γ_51hcos(2π7/365 t) + γ_61hsin(2π7/365 t),

where 𝒮_1(t),…,𝒮_4(t) are deterministic cubic spline basis functions, whose introduction allows moderate deviations from the sinusoidal component. Instead, the second functional class (ℓ=2) characterizes functions having two peaks per year, namely

    θ_2h(t) = ∑_m=1^4 γ_m2h𝒮_m(t) +γ_52hcos(2π14/365 t) + γ_62hsin(2π14/365 t).

The Milan–Rome and Turin–Naples routes, also depicted in Figure <ref>, are potential members of this functional class. As for the prior distributions of the regression coefficients, we let γ_mℓiid∼N(0,1),  for m=1,…,6 and ℓ=1,2, therefore inducing a quite uninformative prior, recalling that the data are standardized. 




We then need to select the functions β(x) and σ(x). In this specific application, this amounts to the specification of four parameters β_1 = β(x_1), β_2 = β(x_2) and σ_1 = σ(x_1), σ_2 = σ(x_2), because the baseline measure P_X is discrete and L = 2. Moreover, in certain clustering applications it is convenient to bound the number of groups to facilitate the practical implementation of marketing campaigns. This leads to the choice of negative values for the discount parameters σ_1 = -β_1 / H_1 and σ_2 = -β_2 / H_2, having set H_1 = 20 and H_2 = 5. In other terms, we allow no more than 20 and 5 clusters within each class, respectively. Moreover, we let β_1 = β_2 = 1 which induces a strong shrinkage on the within-class number of clusters <cit.>, effectively eliminating redundant groups. Hence, the hyperparameters β(x) and σ(x) are selected on the basis of prior information and practical considerations. This is in contrast with the first case study of Section <ref>, in which the hyperparameters α, β(x) and σ(x) were estimated from the data. 

Finally, let us consider the marginal base measure P_X. In mixture of mixture models, the fixed values x_1,…, x_L are arbitrary distinct labels, say the set of numbers {1,…,L}, each with prior proportion α_ℓ/α = E(Π_ℓ) =  P_X(x_ℓ). Moreover, the total mass parameter α = ∑_ℓ = 1^L α_ℓ in this setting can be though as the prior confidence about the latter proportions. We specify a uniform prior for the functional class probabilities (Π_1,Π_2), that is, we let α_1 = α_2 = 1, implying that P_X(x_1) = P_X(x_2) = 1/2. 

Posterior computations are performed using variational Bayes, although Markov Chain Monte Carlo methods may be considered. This results in a point estimate for the nested partition, whose outcomes are briefly summarized in the following. Two representative clusters are shown in Figure <ref>, one for each functional class P_1 (red lines) and P_2 (blue lines). As expected, the functions depicted in red are characterized by single peaked functions, while blue lines display two-peaked functions. The routes depicted in red have a peak of web searches in July, implying that the marketing manager could consider special discounts and offers in that period to incentivize sales. Similarly, marketing actions could be taken in April and subsequently in October for the group of blue routes. Clearly, in Figure <ref> we offer a description of just two specific groups as an exemplification. However, we shall remark that a total of 19 latent trajectories have been found by the aforementioned variational Bayes algorithm, which indeed would require specific analyses and targeted marketing policies. This is carefully detailed in <cit.>, to which we refer for a more structured and detailed explanation in terms of migratory patterns between different geographical locations and further computational details. 




§ CONCLUDING REMARKS


In this paper we extended the notion of enriched nonparametric discrete priors set forth by <cit.>, by considering Pitman–Yor random probability measures in place of Dirichlet processes. The proposed enriched Pitman-Yor process allows for more flexibility and for nested partition structures that can include power-law behaviors; moreover, it provides a unifying probabilistic framework for several existing models within the lively Bayesian nonparametric literature. We illustrated the broad applicability of the epy through two different case studies, both highlighting the usefulness of nested partitions as well as the need for more flexible constructions than the edp. 

Our work underlines a general “enrichment construction” that may lead to further extensions. 
In first place, one could consider a general species sampling model for _X, including the py itself, and the resulting joint random measures  would remain well-defined. Moreover, the urn-scheme presented of Theorem <ref> and the posterior distribution of Theorem <ref> would be straightforward to adapt, as long as P_X is diffuse. The practical advantage of using a dp for _X becomes evident whenever a discrete marginal baseline measure P_X is employed. Indeed, beyond the Dirichlet process case, the discreteness of P_X entails some additional theoretical difficulties, as testified by the work of <cit.>. However, if a py were specified for the marginal random measure _X in place of a dp, the random weights Π_1,…,Π_L of equation (<ref>) would follow the so-called ratio-stable distribution, whose posterior law has been recently obtained by <cit.>. Hence, extensions to general processes for _X would remain tractable even in presence of discrete baseline measures. Moreover, one could also consider more general priors for the random conditional probability measures _Y | X, such as Gibbs-type priors <cit.> and normalized random measures with independent increments <cit.>. As before, the posterior properties derived in Section <ref> may be extended in the aforementioned cases. 

In another research direction, one could consider more complex nested partition mechanisms, aiming at defining random probability measures on product spaces such as X_1 ×⋯×X_d, that is allowing for more than two nesting levels. This idea has been implicitly used in <cit.> for the development of a taxonomic classifier, but a rigorous theoretical study of the involved prior process is currently unavailable. Finally, we shall remark that covariate-dependent extensions of enriched processes can be certainly envisioned, allowing the random variables X_n and Y_n to depend on a set of predictors; refer to <cit.> for a recent and comprehensive overview. Indeed, the exchangeable assumption we made for example in Section <ref> holds only in an approximate sense, because the species occurrence may depend on the spatial location, as pointed out by <cit.>. These considerations encourage further theoretical investigations, for example in the development of space-dependent species sampling models, perhaps building on the contribution of <cit.>.







§ SUPPLEMENTARY MATERIAL







 §.§ Proofs




  §.§.§ Existence of an enriched Pitman–Yor process


Let us denote by B_X and B_Y the Borel σ-fields associated to  X and Y,  respectively. Denote by 𝒬_X the marginal probability law of _X, by 𝒬_Y | X the joint probability law of the conditionals _Y| X(·| x) for any x ∈X, and by 𝒬_X, Y the joint law 𝒬_X, Y = 𝒬_X ×𝒬_Y | X. The proof follows the steps in <cit.>.  First, we obtain that {_Y|X(·| x), x ∈ X} is a set of random conditional probability measures, as we have that  


  (i) each _Y | X(·| x) is a probability measure on X almost surely with respect to 𝒬_Y | X, for any x ∈X;

  (ii) for any Borel set B ∈ B_X, as a function of x,  _Y | X(B | x) is B_X-measurable, a.s. with respect to 𝒬_Y | X.

The property (i) is immediate, whereas (ii) is obtained from results in <cit.> as follows. Let Δ be the set of probability measures on X such that _Y | X(·| x) is measurable as a function of x; results by <cit.> ensure that, if the _Y | X( ·| x) are independent, then the product measure 𝒬_Y | X, obtained via Kolmogorov's existence theorem, assigns outer measure one to Δ. 

Then, let 𝒫_D be the set of discrete probability measures on X.
From the properties of the dp, _X is discrete a.s. with respect to 𝒬_X and by independence of _X and _Y | X it follows that 𝒬_X,Y(𝒫_D ×Δ) = 1. Then, results in  <cit.> ensure that, on 𝒫_D  ×Δ, for any measurable subset A × B of the product space, the map (_X, _Y | X) →∫_A _Y| X(B| x) _X(d x) is jointly measurable in (_X, _Y | X); therefore,  we can define a probability measure 𝒬, on the set of probability measures on the product space X×Y, induced from 𝒬_X,Y restricted to 𝒫_D  ×Δ via the map  (_X, _Y | X) →∫_( ·)_Y | X( ·| x) _X(dx). 






  §.§.§ Proof of Proposition <ref>


The square-breaking representation of the epy follows directly from the stick-breaking representation of the dp and the py that has been recalled in Section <ref>. In particular, note that

    p̃(A × B)    = ∫_A p̃_Y | X(B| x) p̃_X( x) = ∫_A ∑_h=1^∞π_h(x) δ_θ_h(x)(B)p̃_X( x)  
       = ∑_ℓ =1^∞∑_h=1^∞ξ_ℓ π_h(ϕ_ℓ) δ_ϕ_ℓ(A) δ_θ_h(ϕ_ℓ)(B),

for any Borel sets A ⊆X and B ⊆Y. 



  §.§.§ Proof of Theorem <ref>

From equation (<ref>) and recalling Definition <ref> we obtain that the Laplace functional of a Gamma and Pitman–Yor random measure μ̃ can be written as

    E{e^-μ̃(f)} = E[exp{-α∫_Xlog{1 + p̃_Y | X(f| x)}P_X( x)}],

where p̃_Y | X(f| x) = ∫_Yf(x,y)p̃_Y | X( y| x) and for any positive and measurable function f : X×Y→R^+ such that μ̃(f) < ∞ almost surely. The above Laplace functional is fully general and it does not require further restrictions on P_X. Then, by exploiting the discreteness of P_X and the independence among the conditional laws p̃_Y | X we obtain

    E{e^-μ̃(f)}   = E(exp[- ∑_ℓ=1^L α_ℓlog{1 + p̃_Y | X(f| x_ℓ)}]) 
       = ∏_ℓ=1^L E[{1 + p̃_Y | X(f| x_ℓ) }^-α_ℓ],

which concludes the proof for L < ∞. This result can be easily extended to the countable case as an application of dominated convergence theorem. 



  §.§.§ Proof of Theorem <ref>

By definition of the epy and from equation (<ref>), we get that the marginal law p̃_X is independent on Y^(n), given X^(n)= x^(n), so that for any Borel set A ⊆X

    P(X_n+1∈ A |X^(n)   = x^(n), Y^(n)= y^(n)) =
               = E(p̃_X(A) |X^(n)= x^(n), Y^(n)= y^(n))
               = E(p̃_X(A) |X^(n)= x^(n))
               =  P(X_n+1∈ A |X^(n)= x^(n)),

which leads to the <cit.> scheme for the exchangeable sequence (X_n)_n ≥ 1. Note that the (potential) discreteness of P_X poses no issues here. The second part of the theorem is obtained  by exploiting the independence among the py conditional laws p̃_Y | X. Indeed, each subset of observations Y^(n_r) follows the well-known scheme of the py, which is described e.g. in <cit.> when P_Y | X is diffuse. 

We now show that, in fact, the predictive scheme in Theorem <ref> provides a characterization of the epy process, as claimed in Remark <ref>. By Ionescu-Tulcea theorem, the sequence of predictive distributions,  say P(Z_n+1∈·| Z_1=z_1, …, Z_n = z_n), for n ≥ 1 uniquely characterizes the probability law of the stochastic process (Z_n)_n ≥ 1.  Therefore, the sequence of the predictive distributions (<ref>)-(<ref>) characterizes the probability law, say 𝒫, of the stochastic process ( (X_n, Y_n))_n ≥ 1. On the other hand, an exchangeable probability law with an epy directing measure leads to the predictive rule (<ref>)-(<ref>), as shown in Theorem (<ref>). By unicity, the law 𝒫 necessarily coincides with such an exchangeable law.  



  §.§.§ Proof of Theorem <ref>


The independence among p̃_X and Y^(n), given X^(n)= x^(n) immediately leads to the first part of the theorem, thank to conjugacy of the dp <cit.>. Similarly, the posterior distribution of each conditional law p̃_Y | X, thanks to their independence, is obtained as an application of Corollary 20 in <cit.> to each subset of observations Y^(n_r), which requires the diffuseness of P_Y | X. 




 §.§ Choice of the hyperparameters


We discuss here in more detail the prior choices for the three competing models considered in Section <ref>. We compare the estimates α̂, σ̂_r and β̂_r that results from these prior choices with the corresponding maximum likelihood estimates. We show that in most cases the penalization has little effect and therefore it mainly serves as a tool to regularize ill-behaved estimates.










  §.§.§ DP model


In the dp model we need to specify the so-called precision parameter, say β, and we let β∼gamma(2,0.001). The consequent penalized likelihood approach is β̂ = 765.53, whereas the maximum likelihood estimate is β̂_mle = 765.48. Hence, the prior penalty has almost no effect in this specific case.



  §.§.§ EDP and EPY models


The hyperparameter settings for the the edp and the epy models are similar. For the edp, we use α∼gamma(2,0.01), σ_r = 0, and β_r ∼gamma(2,1). For the epy, we let α∼gamma(2,0.01), σ_r ∼beta(1,10), and β_r ∼gamma(2,1). Note that the estimates for α coincide in these two models and one gets α̂ = 11.34 and α̂_mle = 11.24.  Hence, the prior penalty has almost no effect in the estimation of α.

Viceversa, we expect the prior to have an effect in the estimation of β_r and σ_r. As an illustration, we report in Table <ref> the penalized estimates β̂_r, for r=1,…,10,  for the edp model and we compare them with the corresponding maximum likelihood estimates β̂_r,mle. A summary of all the estimates for r=1,…,k_x is given in Table <ref>. Note that the expected value of the gamma(2,1) prior roughly corresponds to the median of the maximum likelihood estimates. As expected, the prior on β_r shrinks the maximum likelihood estimates towards 2 and avoids ill-behaved maximization problems. For example, without the penalty term one would get β̂_3,mle = 0, which is not an admissible value. 

Finally, we report the summary of the estimates (β̂_r,σ̂_r) for the epy model in Table <ref>. The prior on σ_r shrinks the estimates σ̂_r towards 0 and therefore the corresponding β_r parameters are quite similar to those of the edp.

chicago



