1 Introduction↩︎

Historically, Bayesian nonparametric and hierarchical models have addressed data complexity in different ways. Nonparametric models ensure full flexibility to the marginal distribution of the observations by increasing the dimensionality of the parameter space, while hierarchical models focus on the interactions between the observations, grouping them and modeling their dependencies through shared parameters or latent features. The interaction between parameters is recursively modeled in a similar fashion, defining a hierarchical structure that enables pooling of information across different groups while preserving their distinct characteristics – all within the principled framework of Bayesian inference. The hierarchical Dirichlet process [1] represented a significant breakthrough, demonstrating the advantages of combining the two approaches through the sharing of infinite-dimensional parameters. Since then, it has proved effective in an impressive number of contexts, including natural language processing [1], [2], genomics [3]–[5], computer vision [6], [7], music segmentation and speaker diarization [8], [9], cognitive science [10], robotics [11], [12], network analysis [13].

The computational feasibility of the hierarchical Dirichlet process in [1] is strictly linked to a compelling posterior representation via latent variables, often referred to as tables in the restaurant franchise metaphor. The need for tables arises from the nature of the infinite-dimensional parameter, which is an almost surely discrete random probability \(\tilde{P} = \sum_{i \ge 1} J_i \delta_{\theta_i}\), characterized by a countably infinite number of jumps \(J_i\) and atoms \(\theta_i\). In a Bayesian nonparametric setting, a single group of observations is often modeled as conditionally independent and identically distributed from \(\tilde{P}\), where \(\tilde{P}\) is e.g. a Dirichlet process [14] with diffuse mean measure \(P_0\). In this case, since the atoms \(\theta_i\) are independent and identically distributed from \(P_0\), two observations coincide if and only if they share the same atom \(\theta_i\). Conversely, in the hierarchical Dirichlet process, the distribution \(P_0\) of the atoms is itself a latent parameter, which is shared across different groups and likewise modeled as an almost surely discrete random probability. Hence, contrary to the exchangeable case, the atoms display ties with positive probability, and two observations coincide either if they share the same atom or if their respective atoms have identical values. This complicates the posterior representation, unless one keeps track of the atom associated to each observation, which is precisely the role of the latent tables [15]. For \(n\) observations, the tables induce a distribution on the space of partitions of \(n\) elements; the dimensionality of this space increases dramatically with the number of observations, and one typically needs a Gibbs sampling algorithm to sample from it.

The most popular implementation for the hierarchical Dirichlet process is the Gibbs sampler based on the restaurant franchise metaphor [1], often addressed as marginal or collapsed Gibbs sampler. Despite the convenient expressions of its full conditionals and its remarkable flexibility, it is not an exact algorithm, since it samples from a Markov chain and thus approximates the posterior distribution only asymptotically. In addition, there are other well-known drawbacks: it involves a considerable amount of bookkeeping, and scales poorly as the number of observations \(n\) increases, as it relies on a sequential updating scheme for the \(n\) latent tables. Moreover, each table allocation depends on all the other allocations, inducing high autocorrelation in the Markov chain, slow mixing, and preventing parallelization. Such limitations have been recognized by several works; see e.g. [1], [16]–[19]. Accordingly, a plethora of sampling-based strategies have been proposed to reduce computational time and improve mixing and scalability of the standard implementation. These are typically Gibbs samplers adopting a conditional or blocked approach, that is, instantiating the random measures with some finite-dimensional approximation of the posterior. Direct assignment schemes instantiate the jumps of the common random probability, often addressed as global weights, through the stick-breaking construction [1] or a finite-dimensional Dirichlet approximation [9], while variational methods [20]–[22] construct mean field approximations of stick-breaking ratios at both levels of the hierarchy. Nevertheless, their full conditional distributions still depend on the nested partition induced by the tables, and thus do not completely avoid some of its inherent drawbacks. Interestingly, the conditional Gibbs samplers proposed in [18] and [19] obviate the need for tables by considering a finite-dimensional approximation of the model, i.e. by truncating the common random probability a priori.

In this work we pursue a different strategy that eliminates the need for tables while preserving the infinite-dimensionality of the model, thanks to a specific gamma hyperprior for the shared concentration parameter of the Dirichlet processes. Intuitively, the hierarchical Dirichlet process defines a vector of dependent random probabilities with independent jumps and common atoms with ties. Our hyperprior allows for an alternative representation as dependent random probabilities with dependent jumps and common atoms without ties, thereby making the tables superfluous. A fundamental result of this work shows that this representation is the normalization of a specific completely random vector [23], that is, a vector of dependent random measures with jointly independent increments. This is the natural multivariate extension of a completely random measure [24], whose normalization is very popular in the Bayesian nonparametric literature [25]–[27]. In particular, [27] derive an almost conjugate representation of the posterior, conditionally on a real-valued latent variable \(U\) amenable to standard approximate sampling schemes, such as random walk Metropolis-Hastings [28]. In this work, we extend the posterior representation to any normalized completely random vector, which in principle could be applied to many other well-established models, including the normalization of GM-dependent measures [29], Lévy copulas [30], compound random measures [31], and thinned random measures [32].

The characterization as a normalized completely random vector allows to identify a novel posterior representation for the hierarchical Dirichlet process with hyperprior that does not require the latent tables. The price to pay for casting off the tables is the introduction of a latent vector and \(k\) jump vectors, where \(k\) is the total number of distinct observations, typically much smaller than the number of observations \(n\). In principle, the dimension of the latent and jump vectors coincides with the number of groups \(d\), which could be potentially large, thus slowing down the posterior inference algorithm. However, through an in-depth analysis of their distributions, we reduce the non-standard sampling steps to sampling a random vector supported on \([0,\infty)^{2k+1}\). This approach offers three fundamental advantages: (i) the random vector lives in a standard space, and thus can be approximately sampled with standard Markov Chain Monte Carlo (MCMC) techniques, which show good mixing and can be analyzed with a well-established set of diagnostics; (ii) its dimension does not increase with the number of observations \(n\) nor with the number of groups \(d\), which leads to scalable algorithms whenever \(n\) or \(d\) increase but \(k\) remains only moderately large; (iii) when \(n\) is moderately large, we devise an exact i.i.d. sampling algorithm, thus avoiding the approximation error of Gibbs sampling procedures. In summary, our novel posterior representation allows for the development of a scalable algorithm when \(n\) or \(d\) are large compared to \(k\), and of an exact i.i.d. sampling algorithm when \(n\) is only moderately large. Note that \(k\) consistently smaller than \(n\) is the most common setting where the model should be applied, since, e.g., for the hierarchical Dirichlet process, \(k\) is of the same order as \(\log(\log(n))\).

The hierarchical Dirichlet process with hyperprior is only a particular case of the more general class of nonparametric hierarchical models discussed in this work. These models arise as the normalization of a vector \((\tilde{\mu}_1,\dots, \tilde{\mu}_d)\) of conditionally independent completely random measures, given another completely random measure \(\tilde{\mu}_0\). We show that these dependent random measures are completely random vectors and, for this reason, we term them hierarchical completely random vectors. Such hierarchical structures have been used to model dependent hazards [33], [34] and beta processes [35]–[37], but have never been considered for normalization. This novel specification merges two compelling properties for the first time: the naturalness of the hierarchical construction and the analytical tractability of independent multivariate increments. The former leads to simple representations and draws interesting parallels with popular hierarchical models in the literature [1], [38]. The latter enhances the theoretical investigation of the model, leading to the almost-conjugate posterior representation we have already discussed, and playing a crucial role in deriving closed-form expressions for the prior moments and dependence structure.

The work is structured as follows. Section 2 introduces and characterizes hierarchical completely random vectors (hCRVs). Properties of their normalization are studied in Section 3, where we clarify the connections with other hierarchical models in the literature. Section 4 provides a characterization of the posterior distribution for normalized hCRVs. These results are specialized to the gamma-gamma model in Section 5, where we derive our new posterior sampling algorithms and perform relevant comparisons in terms of computational complexity on simulated datasets. Finally, Section 6 sets rules for a fair comparison with the hierarchical Dirichlet process. Using both simulated and real datasets, we argue that a meaningful comparison between dependent priors should entail matching the first two marginal moments and measure of dependence. Proofs and technical details for the sampling algorithms are deferred to the Supplementary Material.

1.0.0.1 Notation

The Cartesian product of \(d\) copies of a set A is denoted as \(A^d = A \times \dots \times A\). The product measure of \(d\) probabilities \(P_1,\dots,P_d\) is \(\prod_{i=1}^d P_i\), or \(P^d\) if \(P_1 = \dots = P_d = P\). Bold symbols indicate vectors, such as \(\boldsymbol{s} = (s_1,\dots,s_d)\), and \(\mathrm{d}\boldsymbol{s} = \mathrm{d}s_1 \cdot \dots \cdot \mathrm{d}s_d\) is the Lebesgue measure on \(\mathbb{R}^d\). The normal distribution with mean \(\mu\) and variance \(\sigma^2\) is denoted by \(N(\mu,\sigma^2)\). If \(f : \mathbb{X} \to \mathbb{Y}\) is a measurable map and \(\rho\) is a measure on \(\mathbb{X}\), then \(f_{\#} \rho\) is the pushforward measure on \(\mathbb{Y}\) defined as \((f_{\#} \rho)(B) = \rho(f^{-1}(B))\), for any measurable set \(B\). The symbol \(\sim\) underlines the randomness of a random measure (\(\tilde{\mu}\)). The abbreviation a.s. stands for almost surely; we repeatedly use \(\Omega_d = [0,+\infty)^d \setminus \{\boldsymbol{0}\}\). Contents in the Supplementary Material are labeled with an ‘S’ prefix, e.g. Section S1.

2 Hierarchical completely random vectors↩︎

In this section, we recall the definition of hierarchical random measures, highlight that they are homogeneous completely random vectors, and recover both their multivariate Laplace exponent and their Lévy measure. The proofs follow a structure similar to the subordination of Lévy processes, first defined in [39], and beautifully described in [40] and [41]. Their joint Lévy measure displays some similarities with that of compound random measures [31], but we show that there is no intersection between the two classes, at least for the same choice of the subordinating random measure. We note that hierarchical completely random vectors can be approximately sampled with standard techniques, exploiting their hierarchical structure. The specification of their law involves an outer and an inner Lévy measure, whose identifiability is studied with relevant examples.

A vector of random measures \(\boldsymbol{\tilde{\mu}} = (\tilde{\mu}_1,\dots,\tilde{\mu}_d)\) is a measurable function on \(M_{\mathbb{X}}^d\), where \(M_{\mathbb{X}}\) denotes the space of boundedly finite measures on a Polish space \(\mathbb{X}\). A completely random vector [23] is the natural multivariate generalization of a completely random measure (CRM), defined in [24]. For a Borel set \(A\) of \(\mathbb{X}\), we use the notation \(\boldsymbol{\tilde{\mu}}(A) = (\tilde{\mu}_1(A),\dots,\tilde{\mu}_d(A))\), which is a random vector in \(\mathbb{R}^d\).

Refer to Section 8 for a brief and self-contained account on completely random measures, Lévy measures, Lévy intensities, Laplace exponents, and their multivariate extension to completely random vectors. Henceforth, \(\tilde{\mu}\sim \mathrm{CRM}(\rho\otimes P_0)\) denotes a CRM with product Lévy intensity \(\mathrm{d}\rho (s) \, \mathrm{d}P_0(x)\), and \(\mathrm{ID}(\rho)\) indicates a pure-jump infinitely divisible distribution with Lévy measure \(\rho\); its expression in integrals is \(\mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(\rho)}\) and its probability density function (p.d.f), if it exists, is denoted by \(f_{\scriptscriptstyle \mathrm{ID}( \rho)}(s)\).

The next theorem shows that the vector of random measures \(\boldsymbol{\tilde{\mu}} \sim \mathrm{hCRV}(\rho,\rho_0,P_0)\) is in fact a completely random vector in the sense of Definition 1. Moreover, the theorem provides the expression of its multivariate Laplace functional through the Laplace exponent and determines its multivariate Lévy intensity; since \(\boldsymbol{\tilde{\mu}}\) is a CRV, such intensity automatically satisfies the integrability conditions (i) and (ii) recalled in Section 8. This construction admits a natural interpretation in terms of subordination. Indeed, the proof follows the same structure as the classical subordination of Lévy processes [39]–[41], and hierarchical CRVs may be regarded as the corresponding random-measure analogue. A related construction, based on the subordination of CRMs by infinitely divisible random measures, is also developed in [42].

In particular, if \(\mathrm{ID}(t \rho)\) has p.d.f. \(f_{ \scriptscriptstyle \mathrm{ID}(t \rho)}\), then \(\rho_{h}\) has a Lévy density \[\rho_{h}(s_1,\dots,s_d) = \int_0^{+\infty} \prod_{i=1}^d f_{\mathrm{ID}(t \rho)}(s_i) \, \mathrm{d}\rho_0(t).\] A sufficient condition for \(\mathrm{ID}(t \rho)\) to have a probability density function, for every \(t > 0\), is that the Lévy measure \(\rho\) has infinite mass and is absolutely continuous [41].

The next result studies the identifiability of the parameters of a hierarchical CRV. For \(c \in \mathbb{R}\), denote by \(c_{\#}\) the pushforward measure of the multiplication map \(s \mapsto cs\).

The last condition can be easily checked on specific classes of models; for example, if we restrict to \(\rho_0 = \alpha \rho\) for some \(\alpha>0\), then \(\rho\) and \(P_0\) are identifiable. To clarify the interpretation of the identifiability conditions, note that \(\rho_0^{(2)}(s) = c \rho_0^{(1)} (cs)\) if and only if \(\tilde{\mu}_0^{(1)} = c \tilde{\mu}_0^{(2)}\).

We conclude this section with two leading examples of hierarchical CRVs.

3 Normalized hierarchical completely random vectors↩︎

One of the most common uses of completely random measures in Bayesian statistics is their normalization [25], which defines random probabilities whose law can act as nonparametric priors. The same construction can be extended to vectors of dependent random measures, such as hierarchical CRVs. In this section, we provide conditions for the normalization to be well-defined and investigate connections with popular models, such as the hierarchical Dirichlet process [1], [49] and the hierarchical normalized \(\sigma\)-stable process [49]. Moreover, we discuss general techniques to measure the dependence of normalized hierarchical CRVs.

For \({\boldsymbol{\tilde{\mu}}} \sim {\rm hCRV}(\rho, \rho_0,P_0)\), we derive a vector of dependent random probabilities as \[\label{eq:norm} \frac{\boldsymbol{\tilde{\mu}}}{\boldsymbol{\tilde{\mu}}(\mathbb{X})} := \bigg(\frac{\tilde{\mu}_1}{\tilde{\mu}_1(\mathbb{X})}, \dots , \frac{\tilde{\mu}_d}{\tilde{\mu}_d(\mathbb{X})}\bigg) ,\tag{1}\] which is well-defined if \(0<\tilde{\mu}_i<+\infty\) a.s., for \(i=1,\dots,d\). The upper bound forces \(P_0\) to be a finite measure; thus, we can assume without loss of generality that \(P_0\) is a probability measure. The lower bound forces each \(\tilde{\mu}_i\) to be infinitely active, that is, the corresponding Lévy measures to have infinite mass.

Therefore, \(\tilde{\mu}_i\) is infinitely active if and only if both \(\tilde{\mu}_0\) and \(\tilde{\mu}_i \mid \tilde{\mu}_0\) are infinitely active.

The construction in 1 is similar to the normalized hierarchical model in [15], [49], where, however, the base random measure is normalized as well, that is \[\label{def:herarchical95structure} \tilde{\mu}_1,\dots,\tilde{\mu}_d \mid \tilde{\mu}_0 \stackrel{\scriptsize{\rm iid}}{\sim}\mathrm{CRM}\bigg(\rho \otimes \frac{\tilde{\mu}_0}{\tilde{\mu}_0(\mathbb{X})} \bigg); \qquad \tilde{\mu}_0 \sim \mathrm{CRM}(\rho_0 \otimes P_0).\tag{2}\] This slight modification has crucial implications on the overall law of \(\boldsymbol{\tilde{\mu}}\), which is no longer a CRV, and whose marginals \(\tilde{\mu}_i\)’s are not CRMs. Interestingly, at least two popular hierarchical specifications can be expressed in terms of a normalized hierarchical CRV, as discussed in the following.

Recall that \(\boldsymbol{\tilde{P}} = (\tilde{P}_1,\dots,\tilde{P}_d) \sim {\rm HDP}(\alpha, \alpha_0,P_0)\) is a hierarchical Dirichlet process [1] with concentration parameters \(\alpha, \alpha_0>0\) and base probability \(P_0\) if \[\label{def:hdp} \tilde{P}_1,\dots,\tilde{P}_d \mid \tilde{P}_0 \stackrel{\scriptsize{\rm iid}}{\sim}{\rm DP}( \alpha \tilde{P}_0); \qquad \tilde{P}_0 \sim {\rm DP}(\alpha_0 P_0),\tag{3}\] where \({\rm DP}(\alpha_0 P_0)\) denotes a Dirichlet process [14] with base measure \(\alpha_0 P_0\). In fact, the normalization of the gamma-gamma hCRV in Example 1 recovers the HDP with a specific gamma prior on the concentration parameter. Here, \({\rm Gamma}(a,b)\) denotes the gamma distribution with shape \(a\) and rate \(b\).

The distribution of a normalized gamma-gamma hCRV depends in fact only on \(\alpha_0\) and \(\alpha/b_0\). The role of the ratio \(\alpha/b_0\) for identifiability is highlighted in Example 1, while \(b\) is a scale parameter for \(\boldsymbol{\tilde{\mu}}\) and disappears with normalization; see also Figure 1 and Section 10.6. In practice, one may restrict to \(b = b_0 = 1\), without loss of generality.

Additionally, when the idiosyncratic component is a stable CRM, the normalized hierarchical model of [49] can be expressed as a normalized hierarchical CRV. This result generalizes, with a different technique, a result of [38] which assumes the base CRM to be stable as well; this same fact is also observed for Lévy processes in [40]. The key property is that, for a stable CRM \(\tilde{\mu}\) with Lévy measure \(\rho\), any proportional measure \(c \tilde{\mu}\) with \(c>0\) is a stable CRM with proportional Lévy measure \(c' \rho\), for some \(c'>0\).

In the following, we provide closed form expression for summary statistics encoding the mean, the variance, and the dependence structure of hCRVs. The simplest and most widely used measure of dependence between random probabilities \(\boldsymbol{\tilde{P}} = \boldsymbol{\tilde{\mu}}/\boldsymbol{\tilde{\mu}}(\mathbb{X})\) is the pairwise linear correlation \(\mathrm{Corr}({\tilde{P}_i}(A),{\tilde{P}_j}(A))\), for a Borel set \(A\). Note that its expression does not depend on \(A\), since normalized hCRVs are normalized homogeneous CRVs and thus belong to the class of multivariate species sampling processes [50]. At the level of the random measures \(\boldsymbol{\tilde{\mu}}\), its computation is a slight modification of [15], where the expressions are derived by leveraging on the conditional independence structure. We report these results in Proposition 18, and present an alternative proof that builds on their joint infinite divisibility. Interestingly, this second technique greatly simplifies the derivation in the normalized case.

4 Posterior representation for normalized hCRVs↩︎

Vectors of dependent random probability measures are commonly employed in Bayesian statistics to model partially exchangeable observations. Indeed, any infinitely active CRV is suitable for this scope through normalization 1 . Many models in the literature fall within this framework, including GM-dependent measures [29], [51], compound random measures [31], Lévy copulas [30], and thinned random measures [32]. In this section, we derive the expression of the posterior distribution for a generic normalized CRV; this can be seen as the multivariate extension of [27] and a special case of FuRBI random measures [52] with shared atoms. These results are later specialized to hierarchical normalized CRVs, and their posterior structure is explored in greater detail.

Let \(\boldsymbol{X}_i = (X_{i1}, \dots, X_{in_i})\) be the \(i\)-th group of observations, for \(i=1,\dots,d\), and consider the model \[\label{eq:model} \boldsymbol{X}_1, \dots, \boldsymbol{X}_d \mid \boldsymbol{\tilde{\mu}} \sim \bigg(\frac{\tilde{\mu}_1}{\tilde{\mu}_1(\mathbb{X})}\bigg)^{n_1} \times \cdots \times \bigg(\frac{\tilde{\mu}_d}{\tilde{\mu}_d(\mathbb{X})}\bigg)^{n_d}; \qquad \boldsymbol{\tilde{\mu}} \sim \text{CRV}(\nu),\tag{4}\] where \(P^m\) denotes the \(m\)-fold product measure and \(\mathrm{d}\nu(\boldsymbol{s}, x) = \mathrm{d}\rho_x(\boldsymbol{s}) \, \mathrm{d}P_0(x)\) is a multivariate Lévy intensity with \(P_0\) a diffuse probability. In the following, we often use the compact notation \(\boldsymbol{X}_{1:d} = (\boldsymbol{X}_1,\dots, \boldsymbol{X}_d)\).

The almost-sure discreteness and the dependence structure of the random measures \(\boldsymbol{\tilde{\mu}}\) imply that the observations \(\boldsymbol{X}_{1:d}\) display tied values with positive probability, both within and across groups. Let \(\boldsymbol{x}^* = (x_1^*, \dots, x_k^*)\) be the distinct values taken by the observed values \(\boldsymbol{x}_{1:d}\), that is, for every \(j=1,\dots,k\), there exist \(x_{ih}\) such that \(x_{ih} = x_j^*\), and \(x_j^* \neq x_{\ell}^*\) for every \(j \neq \ell\). Moreover, for every \(i = 1, \dots, d\) and \(j=1,\dots,k\), denote by \(n_{ij}\) the number of observed values in \(\boldsymbol{x}_i\) equal to \(x_j^*\). This implies that, for every \(i=1,\dots,d\), the decomposition \(n_i = n_{i1} + \cdots + n_{ik}\) holds. Finally, for each \(\boldsymbol{m}\in \mathbb{N}^d\), \(x \in \mathbb{X}\) and \(\boldsymbol{u} \in \Omega_d\), define the cumulant \[\label{eq:cumulants} \tau_{\boldsymbol{m} \mid x}(\boldsymbol{u}) = \int_{\Omega_d} e^{-\boldsymbol{u} \cdot \boldsymbol{s}} \prod_{i=1}^d s_i^{m_i} \mathrm{d}\rho_x(\boldsymbol{s}).\tag{5}\]

The next theorem provides a general expression for the posterior of \(\boldsymbol{\tilde{\mu}} \mid \boldsymbol{X}_{1:d}\), showing that it preserves the CRV property, conditionally on a set of dependent latent variables. For this purpose, define a vector of latent variables \(\boldsymbol{U} = (U_1,\dots, U_d)\) with joint p.d.f. \[\label{eq:latent} f_{\scriptscriptstyle \boldsymbol{U}}(\boldsymbol{u}) \propto e^{- \psi(\boldsymbol{u})} \prod_{i=1}^d u_i^{n_i-1} \prod_{j=1}^k \tau_{n_{1j},\dots, n_{dj}\mid x_j^*}(\boldsymbol{u}),\tag{6}\] where \(\psi\) denotes the Laplace exponent of \(\boldsymbol{\tilde{\mu}}\).

Note that the Lévy measure \(\rho_x\) does not need a density on \(\Omega_d\). This may be seen as a technical detail on \((0,+\infty)\), but it is relevant on \(\Omega_d\), as some popular models, such as GM-dependent measures, do not admit a density; see [53].

We now specialize Theorem 10 to hierarchical CRVs. For simplicity, we assume that \(\rho\) and \(\rho_0\) have Lévy densities on \((0,+\infty)\), denoted with the same notation; this implies that both \(\mathrm{ID}(\rho)\) and \(\mathrm{ID}(\rho_0)\) have a p.d.f. [41]. A first result shows the conditional quasi-conjugacy of the model. Indeed, conditionally on \(\boldsymbol{U}\), we can interpret the vector \(\boldsymbol{\tilde{\mu}}^*\) as a hierarchical CRV with heterogeneous marginal distributions.

In order to sample from the posterior, one has to sample the jumps \(\boldsymbol{J}_j\)’s at fixed locations and the latent variables \(\boldsymbol{U}\), which are \(d\)-dimensional random variables. When \(d\) is small, the sampling task can be performed via \(d\)-dimensional rejection sampling or approximated by a Metropolis-Hastings scheme. However, when \(d\) is moderate or large, reducing the dimension of the proposal becomes essential to preserve efficiency. For hierarchical CRVs, the sampling of jumps can be reduced to the sampling of \(1\)-dimensional random variables. To this end, define the quantity \[\label{eq:tau95bar} \bar \tau_m(u, t) = \int_0^{+\infty} s^m e^{-u s} f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s) \, \mathrm{d}s.\tag{8}\] Remarkably, the quantity above may be computed without integration. Indeed, whenever an explicit expression for the Laplace transform of \({\rm ID}(t\rho)\) is available, one may instead compute its derivatives, since \[\bar \tau_m(u, t) = (-1)^m \frac{\mathrm{d}^m}{\mathrm{d}u^m} \int_0^{+\infty} e^{-u s} f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s) \, \mathrm{d}s.\]

Therefore, Proposition 12 reduces the sampling of each \(d\)-dimensional jumps vector to the easier task of sampling \(d\) conditionally independent \(1\)-dimensional jumps, given an additional variable. These random variables can be sampled exactly in the case of gamma-gamma hCRVs, as detailed in Section 5.

As for the latent variables \({\boldsymbol{U}}\), the evaluation of their p.d.f. 6 up to a normalizing constant requires to compute the cumulants \(\tau_{\boldsymbol{m} \mid x}\) in 5 . Although this task can be tackled on a case-by-case basis, the \(d\)-dimensional integral in their definition can be reduced to a \(1\)-dimensional integral for hierarchical CRVs. Given their homogeneity, the cumulant \(\tau_{\boldsymbol{m} \mid x} = \tau_{\boldsymbol{m}}\) does not depend on \(x \in \mathbb{X}\).

Interestingly, the cumulant \(\tau_{n_{1j},\dots, n_{dj}}(\boldsymbol{U})\) is the normalizing constant for the density of the jump \(J_{0j}\) in Proposition 12. These quantities represent the building blocks for any computational method approximating the \(d\)-dimensional distribution of latent variables \(\boldsymbol{U}\). However, when \(d\) is large, reducing to lower-dimensional distributions becomes essential. Therefore, we propose an alternative representation of these latent variables as conditionally independent gamma random variables with random rate parameters; this result may also be of interest in the univariate case. For \(\boldsymbol{m} \in \mathbb{N}^k\) and \(\boldsymbol{t} \in (0,+\infty)^{k+1}\), set \(m_\bullet = \sum_{j=1}^k m_j\) and define the integral \[\label{eq:norm95constant} C(\boldsymbol{m}; \boldsymbol{t}) = \int_{(0,+\infty)^{k+1}} (s_0 + s_1 + \cdots + s_k)^{-m_{\bullet}} f_{\scriptscriptstyle \mathrm{ID}(t_0 \rho)}(s_{0}) \prod_{j=1}^k s_{j}^{m_{j}} f_{\scriptscriptstyle \mathrm{ID}(t_j \rho)}(s_{j}) \,\mathrm{d}\boldsymbol{s}.\tag{9}\]

Proposition 13 reduces the sampling of the \(d\)-dimensional vector of dependent latent variables \(\boldsymbol{U}\) to the sampling of \(d\) conditionally independent \(1\)-dimensional random variables \(\beta_1,\dots,\beta_d\), which represent the scale parameters of \(U_1, \dots, U_d\). This sampling task can be approached through standard computational methods or further simplified, as described in the following section.

5 The normalized gamma-gamma hCRV↩︎

5.1 Posterior representation↩︎

This section specifies results in Section 4 to the fundamental example of gamma-gamma hCRVs, and provides details for the subsequent implementation of posterior sampling algorithms. In the following, \(((t))_n = \Gamma(t+n)/\Gamma(t)\) denotes the ascending factorial.

Details on sampling algorithms for the hierarchy in (a) are given in Sections 10.4 and 10.5, where we describe an efficient implementation of Newton’s method for inverting the exponential integral, based on a logarithmic transformation that guarantees convergence for each starting point. Notably, this same algorithm can be readily employed to sample a priori from a gamma-gamma hCRV; see Remark 3 and Section 10.4. Moreover, the nontrivial task of sampling the jumps \(\boldsymbol{J}_j\)’s is here reduced to sampling variables \(\alpha J_{0j}\)’s in (c), whose densities can be computed up to normalizing constants. A convenient option is to resort to computational methods that generate approximate samples, such as Metropolis-Hastings algorithms. Details are provided in Section 10.1, where we empirically show that the random-walk Metropolis-Hastings scheme on the log-scale with Gaussian increments outperforms the approach proposed in [28], based on gamma proposals. Remarkably, we also develop an exact sampler in Section 5.2.

The remaining step to sample from the posterior is sampling the latent variables \(\boldsymbol{U}\), whose characterization in Proposition 13 can be simplified through a change of variables.

For implementation convenience, one can directly sample \(U_i/b\). Indeed, these are the relevant quantities in Proposition 14 for sampling the jumps \(J_{i1}, \dots, J_{ik}\) and the random measure \(\tilde{\mu}^*_i\), and computing \(\lambda(\boldsymbol{U})\). The non-standard step for sampling \(\boldsymbol{U}\) is the marginal sampling of the variable \(\alpha T\) in ?? , whose joint density with the auxiliary vector \(\boldsymbol{V}\) is known up to a normalizing constant. In Section 10.1 we describe a Metropolis-Hastings within Gibbs algorithm to sample from the marginal distribution of \(\alpha T\). Alternative procedures to obtain exact samples from \(\alpha T\) are considered in Section 5.2.

5.2 Exact sampling↩︎

As highlighted by Propositions 14-15, posterior sampling from the normalized gamma-gamma hCRV mainly consists in sampling gamma random variables, except for two critical steps, namely the sampling of random variables \(\alpha J_{01}, \dots, \alpha J_{0k}\) in ?? and the marginal sampling of random variable \(\alpha T\) in ?? . In both cases, one may resort to MCMC algorithms, based on Metropolis-Hastings steps, to obtain approximate samples from such distributions. In the following, we propose alternative algorithms based on analytical manipulations of their density functions, which instead allow for exact sampling. For this purpose, introduce the coefficients \(S(q_1,\dots,q_d; h)\), defined by the recursive relation \[\label{eq:multivariate95stirling} S(q_1,\dots,q_\ell+1,\dots,q_d; h) = q_\ell S(q_1,\dots,q_d; h) + S(q_1,\dots,q_d; h-1),\tag{10}\] with boundary conditions \(S(0,\dots,0; 0) = 1\) and \(S(q_1,\dots,q_d; 0) = S(0,\dots,0; h) = 0\) for \(q_{\bullet}> 0\) or \(h > 0\). Remarkably, these coefficients are the natural generalization of the unsigned Stirling numbers of the first kind to blocked permutations. Indeed, \(S(q_1,\dots,q_d; h)\) represents the number of permutations with \(h\) cycles in the Young subgroup of \(\mathcal{S}_{q_\bullet}\) (group of permutations of \(q_\bullet\) elements) induced by the integer partition \((q_1,\dots, q_d)\). Henceforth, we refer to the coefficients \(S(q_1,\dots,q_d; h)\) as multivariate Stirling numbers. For convenience, let \(m_{ij} = \min(1,n_{ij}) \in \{0,1\}\), and define \(m_{\bullet j} = \sum_{i=1}^d m_{ij}\).

Therefore, the random variable \(\alpha J_{0j}\) is in fact a finite mixture of gamma distributions, with mixing weights \(p_{jh_j} \propto \lambda(\boldsymbol{U})^{-h_j} \, \Gamma(h_j) \, S(n_{1j},\dots,n_{dj}; h_j)\), for \(h_j = m_{\bullet j}, \dots, n_{\bullet j}\). The evaluation of the \(S(n_{1j},\dots,n_{dj}; h_j)\)’s via the recursive relation in 10 has computational cost \(\mathcal{O}\big(n_{\bullet j}^2\big)\); see Remark 19.

The evaluation of the \(c_h\)’s may seem computationally expensive when \(k\) is large. Indeed, a naive approach would involve \(k\) nested cycles, with a computational cost of \(\mathcal{O}\big( \prod_{j=1}^k n_{\bullet j}\big)\). However, these coefficients can be obtained through a sequence of discrete convolutions and computed at cost \(\mathcal{O}\big(\sum_{j < \ell} n_{\bullet j} n_{\bullet \ell} \big)\), as detailed in Section 10.2. Exact samples from the distribution of \(\alpha T\) can be obtained resorting to rejection sampling algorithms. Indeed, considering a real parameter \(r\), ?? can be rewritten as \[f_{\scriptscriptstyle \alpha T}(t) \propto t^{\alpha_0+r-1} \, e^{-(b_0/\alpha)\,t} \, (t^{-r} R(t)), \qquad R(t) = \prod_{i=1}^d \frac{1}{((t))_{n_{i}}} \left( \sum_{h=m}^n \frac{c_h}{((\alpha_0))_h} \, t^h \right),\] where \(R(t)\) is a ratio of polynomials, and continuous for \(t \ge 0\). The distribution of \(\alpha T\) is absolutely continuous with respect to a \({\rm Gamma}(\alpha_0 + r, b_0 / \alpha)\), with Radon-Nikodym derivative proportional to \(t^{-r} R(t)\), from which we construct a rejection sampling algorithm. A necessary condition for the Radon-Nikodym derivative is to be bounded above: this is guaranteed when \(0 \le r \le m-d\). Within this interval, we choose the value of \(r\) that maximizes the acceptance probability, which is equivalent to maximizing \[r \log t^*(r) - \log R(t^*(r)) + (\alpha_0 + r) \log(b_0/\alpha) - \log \Gamma(\alpha_0 + r),\] where \(t^*(r)\) is the value of \(t\) that maximizes \(t^{-r} R(t)\). Details on the optimal choice of \(r\) can be found in Section 10.3. Alternatively, the adaptive rejection sampling algorithm of [54] may considerably improve the acceptance rate. Although computationally convenient, theoretical guarantees are only available if the density of \(\alpha T\) is logarithmically concave, which we were unable to prove in general. Therefore, for numerical experiments, we consider the standard rejection sampler outlined above.

5.3 Posterior sampling algorithms↩︎

Collecting the distributional results described above, as well as the details in Section 10, we obtain algorithms to sample jumps \(\boldsymbol{J}\) at fixed locations in the posterior distribution (Proposition 11) induced by the gamma-gamma hCRV. The Julia implementation of these algorithms is available at github.com/claudiodelsole/hCRV.jl, together with an interface to allow their integration within the R environment. The common structure of conditional dependencies within the sampling procedures is summarized in Figure 1, where the computational bottlenecks are highlighted by red circles. Remarkably, the graph of dependencies has no proper cycles, i.e. it is a directed acyclic graph (DAG): variables are sampled in topological order, with each variable depending only on (some of) its ancestors. Such dependence structure may enhance mixing. Furthermore, the proposed algorithms allow for parallelization at different stages. Indeed, after the initial sampling of \(\alpha T\) from its marginal, one can parallelize across groups to sample the latent variables \(\boldsymbol{U}\) (Proposition 15). Likewise, after the computation of \(\lambda(\boldsymbol{U})\), one can parallelize across distinct values to sample the conditionally independent variables \(\alpha J_{01}, \dots, \alpha J_{0k}\) and further parallelize across groups and distinct values to sample jumps \(\boldsymbol{J}\) (Proposition 14).

The MCMC scheme for posterior sampling the jumps \(\boldsymbol{J}\) is summarized in Algorithm 2. The state of the Markov chain consists of the latent variable \(\alpha T\) and jumps \(\alpha J_{01}, \dots, \alpha J_{0k}\), plus auxiliary variables \(\boldsymbol{V}\), thus having dimension \(2k+1\). The exact posterior sampling procedure developed in Section 5.2 is instead summarized in Algorithm 3. This scheme outputs i.i.d. samples from the posterior, but requires a non-trivial initialization. The main computational bottleneck is represented by the rejection sampling step for the latent variable \(\alpha T\), which may suffer from low acceptance rates. Numerical illustrations proving the effectiveness of the proposed algorithms are discussed in Section 10.8.

As highlighted in Proposition 6, the normalized gamma-gamma hCRV model coincides with the hierarchical Dirichlet process (HDP) with a particular gamma prior on the concentration parameter. Therefore, the algorithms described above can be compared with standard samplers for posterior analysis of the HDP. Among the various alternatives and their countless variations, we focus on marginal algorithms, and consider (i) the Gibbs sampler based on the restaurant franchise metaphor [1], in which the allocation to dishes is observed (and thus fixed), and (ii) a collapsed Gibbs sampler for the number of tables serving each dish, derived from the pEPPF; cfr. Section 10.7 for details. The latter algorithm does not appear in the existing literature, and may be of independent interest. In fact, it is related to a sampler considered in [55] for the hierarchical Pitman-Yor process [49], but relies on a further marginalization of the latent tables through the multivariate Stirling numbers in 10 , reducing from \(kd\) to \(k\) latent variables. The marginal sampler (i) involves \(n\) latent variables, one for each observation’s allocation to a table; hence, its computational complexity rapidly increases with the sample size. In this respect, the collapsed Gibbs sampler (ii) provides a substantial improvement, since its state space has dimension \(k\), which is typically much smaller than \(n\). However, updating each of the \(k\) variables (namely the number of tables serving each dish) requires the evaluation of the multivariate Stirling numbers in 10 . These grow extremely fast and often cause numerical overflows even for moderately large sample sizes. For this reason, when \(n\) is (moderately) large, the marginal sampler (i) remains the practical choice.

In fact, our proposals present some interesting advantages. Firstly, in presence of a (moderately) large number of observations, the MCMC approach in Algorithm 2 relies on a state of dimension \(2k+1\), thereby reducing the computational burden of the marginal sampler (i), without requiring the evaluation of multivariate Stirling numbers. Moreover, the structure of conditional (in)dependencies enables parallelization and potentially speeds up mixing. Finally, running a Markov chain on the Euclidean space \((0,\infty)^{k+1} \times \Delta^k\) is more convenient than sampling on the constrained space of partitions: many standard tools for MCMC diagnostics are designed for Euclidean spaces, and acceptance rates of Metropolis-Hastings steps can be optimized by tuning the variance of the proposal; see Section 10.1. On the other hand, when \(n\) is small or moderately large, so that the computation of multivariate Stirling numbers in the collapsed sampler (ii) remains feasible, Algorithm 3 allows for exact i.i.d. sampling from the posterior. This approach avoids any potential issue inherent to MCMC schemes, without increasing the computational complexity.

5.4 Simulation studies↩︎

In the following, we carry out several simulation studies to compare the algorithms in terms of their execution time, as the dimensions of the input data grow. Specifically, we compare four algorithms: (i) the MCMC sampler in Algorithm 2 with symmetric random-walk Metropolis-Hastings steps on the log-scale (mcmc); (ii) the exact sampler in Algorithm 3; (iii) the marginal Gibbs sampler of [1] based on the restaurant franchise metaphor (CRF); (iv) the collapsed Gibbs sampler for the number of tables serving each dish, detailed in Section 10.7. Algorithms are compared in terms of execution time per effective sample, averaging over \(100\) simulated datasets for each experimental setting. The burn-in time for algorithms (i), (iii) and (iv) and the initialization time for algorithm (ii) are deducted from the total execution times. The exact and collapsed Gibbs samplers are stopped when the computations of the multivariate Stirling numbers in 10 encounter a numerical overflow. The corresponding execution times are plotted when at least 75 out of 100 experiments are completed without errors.

The left panel of Figure 4 compares the algorithms for increasing number of groups \(d\), while the number of observations per group \(n_{i} = 50\) is fixed. Similarly, the right panel compares the algorithms for increasing number of observations per group \(n_{i}\), while the number of groups \(d = 20\) is fixed. In both cases, observations are sampled from a hierarchical Dirichlet process with concentration parameters \(\alpha = 5\) and \(\alpha_0 = 3\). Note that parameter \(\alpha_0\) impacts the number of columns in the counts matrix \((n_{ij})_{ij}\), while parameter \(\alpha\) controls its sparsity. For the CRF-based Gibbs sampler, the time per effective sample grows more than linearly in the total number of observations, as expected given its sequential allocation structure. The other algorithms show a linear growth rate in the number of groups (Figure 4, left panel), with the exact and collapsed samplers displaying a smaller slope compared with the MCMC approach. These algorithms appear less affected by the number of observations per group, showing a nearly constant behavior as \(n_i\) grows (Figure 4, right panel). We argue that their execution times may instead mostly depend on the number of distinct values \(k\), which determines the dimension of the state space for both the MCMC and collapsed samplers. In fact, \(k\) is slowly increasing in the experimental settings considered above, as displayed in the bottom panels of Figure 4. Refer to Section 10.9 for a further simulation study supporting this claim.

In conclusion, the time complexity of both proposed sampling schemes scales linearly in the number of groups \(d\) and distinct values \(k\), while it is essentially unaffected by the number of observations. In our implementation, the exact sampler shows better performances, yielding around twice as many independent posterior samples as the MCMC approach, in the same computation time. This is comparable with the collapsed Gibbs sampler for the HDP, which however does not output i.i.d. samples. Despite their efficiency, both algorithms suffer numerical issues as the sample size grows: in our experimental setting, their applicability is limited to \(500\)-\(700\) observations. On the contrary, the MCMC algorithm runs efficiently regardless of the sample size. In particular, it proves much faster than the CRF-based Gibbs sampler for moderate to large sample sizes.

5.5 Extensions to other hCRVs↩︎

The normalized gamma-gamma hCRV represents both the most relevant and the analytically simplest hCRV construction, since Proposition 6 establishes its tight connection with the HDP, arguably the cornerstone of Bayesian nonparametric dependent models, and several latent variables involved in its posterior characterization collapse to gamma distributions. Nevertheless, some analytical tractability is preserved beyond the gamma-gamma case, with potential advantages in terms of model flexibility, e.g. in the asymptotic behavior of the number of distinct values. We argue that the tractability of the gamma-gamma hCRV is primarily due to the idiosyncratic component \(\rho\) being the Lévy measure of a gamma CRM. Indeed, the results in Propositions 12 and 13 heavily depend on the distributional form of \({\rm ID}(t \rho)\), while the choice of \(\rho_0\) only affects the distribution of the \(J_{0j}\)’s and of the vector \(\boldsymbol{T}\). Therefore, one may consider different choices of the base Lévy measure \(\rho_0\) to enhance model flexibility, at a reduced cost in terms of tractability; for example, \(\rho_0\) can be taken to be the Lévy measure of a generalized gamma or stable CRM. The main limitation is the need for pointwise evaluations of the density of \({\rm ID}(\rho_0)\), whose analytical form is rarely available; however, this task can be performed through well-established numerical methods. Instead, the extension of the idiosyncratic measure \(\rho\) to other choices is less straightforward; in particular, within the present posterior representation, the availability of a closed-form expression for the density of \({\rm ID}(t \rho)\) for \(t > 0\) is essential to sample the jumps \(J_{ij}\)’s in Proposition 12, as well as the vectors \(\boldsymbol{S}_i\)’s in Proposition 13. As a result, reasonable fully tractable choices for \(\rho\) are limited to few special cases, namely the inverse Gaussian process and \(\sigma\)-stable process with \(\sigma = 0.5\). Nonetheless, the possibility to directly sample the \(\beta_i\)’s in Proposition 13 is limited to the gamma case. The development of alternative techniques to handle \({\rm ID}(t \rho)\) beyond these specifications would considerably enlarge the class of tractable hCRVs, and represents an interesting direction for future work. In conclusion, we point out that there are no particular reasons to constrain \(\rho\) and \(\rho_0\) to take the same analytical form, neither from the modelling nor the computational perspectives. Instead, their choice should be informed by the desired flexibility in the marginal and dependence structures; see Proposition 9 and the discussion in [15].

6 Application: a fair comparison with the HDP↩︎

Section 11 investigates the impact of the hyperparameters on the borrowing of information and shrinkage induced by the normalized gamma-gamma hCRV, and discusses the proper elicitation of its prior dependence structure. In this section, we compare such model with the hierarchical Dirichlet process. We consider a simulated dataset with \(d = 3\) groups of independent Poisson observations, each of size \(n_{i} = 10\), with means equal to \(2\), \(3\) and \(4\). Proposition 6 shows that the normalized gamma-gamma hCRV is equivalent to the HDP with a particular gamma prior on the concentration parameter \(\alpha\). Hence, it may seem natural to compare the normalized gamma-gamma hCRV with the HDP by fixing the same parameters \((\alpha, \alpha_0, P_0)\), as in Figure 5 (left plot), where \(P_0 = N(8,1)\). However, this approach may alter the comparison, since fixing the same parameters can entail very different marginal variance and correlation structures. In particular, for the HDP\((\alpha,\alpha_0,P_0)\), following [49], \[\begin{gather} E(\tilde{P}_i(A)) = P_0(A), \qquad \mathrm{Var}(\tilde{P}_i(A)) = \frac{1 + \alpha + \alpha_0}{(1+\alpha_0)(1 + \alpha)} \, P_0(A)(1-P_0(A)), \\ \mathrm{Corr}(\tilde{P}_i(A), \tilde{P}_j(A)) = \frac{1 + \alpha}{1 + \alpha + \alpha_0}. \end{gather}\] Table ¿tbl:table:limits? reports some limiting behaviors, which notably differ from those of the normalized gamma-gamma hCRV, derived from Example 3. Therefore, rather than considering the same values of \((\alpha, \alpha_0)\), a fair comparison should consider the same values of the mean, variance and correlation. Specifically, one should choose \(\sigma^2, \rho \in (0,1)\) such that \(\mathrm{Var}(\tilde{P}_i(A)) = \sigma^2 \, P_0(A)(1-P_0(A))\) and \(\mathrm{Corr}(\tilde{P}_i(A), \tilde{P}_j(A)) = \rho\), and set the hyperparameters of the models accordingly. For the normalized gamma-gamma hCRV, the corresponding values of \(\alpha\) and \(\alpha_0\) are obtained by solving the system of non-linear equations \[\rho\, \big(1+\alpha_0/\alpha \, e^{1/\alpha} E_{\alpha_0}(1/\alpha)\big) - 1 = 0, \qquad \sigma^2 (1 + \alpha_0) - 1/\rho = 0,\] using standard numerical methods. The same task is achieved for the HDP by setting \[\alpha = \frac{1}{1-\rho} \bigg( \frac{1}{\sigma^2}-1 \bigg), \qquad \alpha_0 = \frac{1}{\rho \sigma^2} -1.\] Figure 15 displays the values of \((\alpha, \alpha_0)\) for the two models as functions of \((\sigma^2, \rho)\). The right plot of Figure 5 compares the normalized gamma-gamma hCRV with the HDP fixing the same mean, variance, and correlation. In contrast, the middle plot illustrates a common approach in which the correlation is fixed without accounting for variance. In conclusion, when fixing the same parameters \((\alpha, \alpha_0)\), one might mistakenly conclude that the gamma-gamma hCRV borrows much more information than the HDP. A similar misinterpretation arises when fixing the correlation without adjusting for variance. Ultimately, by matching variance and correlation, it becomes clear that the gamma-gamma hCRV exhibits only slightly more borrowing of information. This increase can be attributed to the hierarchical CRV leveraging not only the common base random probability but also information contained in the random concentration parameter.

A similar analysis is reproduced on the Palmer Archipelago (Antarctica) penguin data [56], with the goal of predicting the flipper length for female penguins. The dataset includes three species of penguins (Adelie, Chinstrap, Gentoo) with many ties (\(n = 165\), \(k = 41\)). The prior guess is \(P_0 = N(100, 10)\), which provides a significant bias towards lower values of the flipper length. As shown in Figure 6, fixing the same value of the hyperparameters \((\alpha, \alpha_0)\) could incorrectly suggest that the gamma-gamma hCRV gives too much weight to the prior with respect to the HDP (left plot). This effect is mitigated by fixing the same correlation but different variance (middle plot) and disappears by fixing the same variance and correlation (right plot). Under this fair model comparison, the gamma-gamma hCRV borrows slightly more information than the HDP. This behavior is the same observed for the simulated data, and can be ascribed to the information contained in the common concentration parameter.

7 Discussion↩︎

This work introduces normalized hierarchical completely random vectors (hCRVs), a new way of building dependent priors that combines the naturalness of hierarchical structures with a convenient posterior representation due to multivariate infinite divisibility. On the one hand, they provide an alternative sampling strategy for existing models, such as the hierarchical Dirichlet process [1]; on the other hand, they represent a general recipe for constructing dependent priors with a convenient posterior representation. As such, normalized hCRVs enter the flourishing context of partially exchangeable models, which includes, among others, dependent stick-breaking constructions [57]–[60], additive random measures [29], [61], nested structures [62]–[65], Lévy copulas [30], [66], compound random measures [31], dependent Pólya trees [67], FuRBI priors [52]; see [68] for a recent review.

Future work will explore additional theoretical properties and methodological advantages of normalized hCRVs in different settings and for different families of CRMs. A first natural development is the derivation of the induced random partition through the pEPPF and of the predictive distribution for each group of observations, in the spirit of [49]. Moreover, in the context of multigroup data, the combination of normalized hCRVs with probability kernels in mixture models could lead to a novel class of priors for density estimation and clustering [69]. This raises interesting questions concerning the contraction rate to the true densities for partially exchangeable models [70], the behavior of the mixing measure [71], and that of the clustering structure [72]. Additionally, normalized hCRVs could be employed in stochastic block models for multi-layer network data [13], and, in the one-dimensional setting (\(d = 1\)), they could represent a viable alternative as priors beyond Gibbs type [38] or as building blocks for multiview clustering [73]. Finally, from a numerical perspective, it would be of interest to evaluate the performance of hCRVs in more structured application domains, such as topic modelling, genomics, or sequential data, in order to assess the regimes in which the table-free representation yields the largest practical gains in terms of mixing, scalability, and computational cost.

Funding↩︎

Claudio Del Sole was partially supported by the European Union - Next Generation EU PRIN-PNRR (project P2022H5WZ9).

8 Background on completely random measures↩︎

This appendix contains a concise account of completely random measures [24] and their multivariate extension to vectors of random measures. In particular, we introduce the key notions of Lévy measure and Laplace exponent.

Let \((\mathbb{X}, d_{\mathbb{X}})\) be a Polish space endowed with a distance \(d_{\mathbb{X}}\). The space \(M_{\mathbb{X}}\) of boundedly finite measures on \(\mathbb{X}\) is a Borel space with the topology of weak\(^\sharp\) convergence [74]. A random measure is a measurable function \(\tilde{\mu}: \Omega \to M_{\mathbb{X}}\) from some probability space \(\Omega\).

[24] provides a remarkable representation theorem that decomposes any CRM \(\tilde{\mu}\) into a unique sum \(\tilde{\mu}= \mu_0 + \tilde{\mu}_f + \tilde{\mu}_c\) of three independent components: (i) a deterministic measure \(\mu_0\), (ii) a random measure with fixed atoms \(\tilde{\mu}_f\), and (iii) an a.s. discrete random measure without fixed atoms \(\tilde{\mu}_c\). The use of CRMs as priors in Bayesian nonparametric models usually restricts to the third component. For this reason, unless differently specified, we only focus on this class of CRMs and henceforth assume \(\tilde{\mu}= \tilde{\mu}_c\). Moreover, [24] shows that for any CRM \(\tilde{\mu}\) there exists a Poisson random measure \(\tilde{N}\) on \((0,+\infty) \times \mathbb{X}\) such that, for any Borel set \(A\) of \(\mathbb{X}\), \[\label{eq:poisson} \tilde{\mu}(A) = \iint_{(0,\infty) \times A} s \, \mathrm{d}\tilde{N} (s, x).\tag{11}\] The mean measure \(\nu\) of the Poisson random measure \(\tilde{N}\) characterizes the law of \(\tilde{\mu}\) and is termed the Lévy intensity of the CRM, justifying the notation \(\tilde{\mu}\sim {\rm CRM}(\nu)\). This measure on \((0,+\infty) \times \mathbb{X}\) can have infinite mass on sets of the form \((0,\varepsilon) \times A\). However, for every \(\varepsilon>0\) and for every bounded Borel set \(A\), it must satisfy the following constraints:

1. the jump component is bounded out of the origin, \(\nu((\varepsilon,+\infty) \times A) <+\infty\);

2. the jump component is integrable in the origin, \(\iint_{(0,\varepsilon)\times A} s \, \mathrm{d}\nu(s,x) <+\infty\);

3. the atom component has no point masses: for every \(x \in \mathbb{X}\), \(\nu((0,+\infty) \times \{x\}) = 0\).

Details can be found in [75]. The identity 11 is used in [24] to derive a Lévy-Khintchine representation of the Laplace transform \[\label{eq:levy-khintchine} \log \mathbb{E}(e^{-\lambda \tilde{\mu}(A)}) = - \iint_{(0,+\infty) \times A} (1-e^{-\lambda s}) \, \mathrm{d}\nu(s,x),\tag{12}\] for every \(\lambda \ge 0\). Therefore \(\tilde{\mu}(A)\) has a pure-jump infinitely divisible distribution with Lévy measure \(\mathrm{d}\rho_A(s) = \int_{A} \mathrm{d}\nu(s, x)\), which we compactly write as \(\tilde{\mu}(A) \sim {\rm ID}(\rho_A)\); refer to [41] for further details. For \(A = \mathbb{X}\), the expression in 12 provides the Laplace exponent of the CRM, namely the function \(\psi: [0,+\infty) \to [0,+\infty)\) defined as \[\psi(\lambda) = - \log \mathbb{E}(e^{-\lambda \tilde{\mu}(\mathbb{X})}).\] Remarkably, thanks to the independence of the evaluations on disjoint sets, the Laplace transform characterizes the law of the CRM. The usual specification of Lévy intensities is through disintegration \[\label{eq:disintegration} \mathrm{d}\nu(s,x) = \mathrm{d}\rho_x(s) \,\mathrm{d}P_0(x),\tag{13}\] where \(P_0\) is a \(\sigma\)-finite atomless measure on \(\mathbb{X}\) for condition (c), and \(\rho_x\) is a \(\sigma\)-finite measure on \((0,+\infty)\), \(P_0\)-a.e., such that, for every \(\varepsilon>0\), \[\label{eq:conditions95levy} \text{(a)} \quad \rho_x((\varepsilon,+\infty)) <+\infty, \qquad \text{(b)} \quad \int_{(0,\varepsilon)} s \, \mathrm{d}\rho_x(s) <+\infty.\tag{14}\] Note that \(P_0\) does not have to be a probability measure nor a measure with finite mass, unless we extend condition (a) to unbounded Borel sets. In this case, we can assume \(P_0\) to be a probability measure, and the disintegration 13 is unique. Conditions 14 imply that \(\rho_x\) is a Lévy measure on \((0,+\infty)\); when \(\rho_x\) is absolutely continuous with respect to the Lebesgue measure, we term its Radon-Nikodym derivative Lévy density. In our context, we are interested in two additional conditions. Firstly, we consider Lévy measures \(\rho_x\) that have infinite mass near the origin, that is, such that \(\rho_x((0,+\infty)) = \infty\); we term such Lévy measures infinitely active. Secondly, we consider disintegration in product form, which leads to the definition of homogeneous CRM.

For a homogeneous \(\tilde{\mu}\sim {\rm CRM}(\rho \otimes P_0)\) the Lévy-Khintchine representation in 12 simplifies consistently. Indeed, under homogeneity and assuming \(P_0\) to be a probability measure, the Laplace exponent is \[\psi(\lambda) = \int_{(0,+\infty)} (1-e^{-\lambda s}) \, \mathrm{d}\rho(s),\] and \(\log \mathbb{E}(e^{-\lambda \tilde{\mu}(A)}) = - P_0(A) \, \psi(\lambda)\), for \(\lambda \ge 0\). We observe that \(\psi\) is a non-negative and infinitely differentiable function whose derivative is completely monotone. Moreover, it vanishes at 0 and its derivative vanishes at \(+\infty\). Thanks to the Lévy-Khintchine representation of Bernstein functions (see, e.g. [76], Theorem 3.2) any such function is the Laplace exponent of a CRM.

The two fundamental examples are the gamma CRM and the stable CRM. We recall their definitions, which are used as building blocks throughout the manuscript.

The notion of completely random measure can be naturally extended to vectors \(\boldsymbol{\tilde{\mu}} = (\tilde{\mu}_1,\dots,\tilde{\mu}_d)\) of random measures, as in Definition 1. Following [23], we term them completely random vectors (CRVs). The representation in 11 can be extended to CRVs by considering a multivariate Lévy intensity \(\nu\) on \(\Omega_d \times \mathbb{X}\), where \(\Omega_d = [0,+\infty)^d \setminus \{\boldsymbol{0}\}\). Similarly, the Laplace transform of the random vector \(\boldsymbol{\tilde{\mu}}(A)\) is characterized by \[\label{eq:levy-khintchine95crv} \log \mathbb{E}(e^{-\boldsymbol{\lambda} \cdot \boldsymbol{\tilde{\mu}}(A)}) = - \iint_{\Omega_d \times A} (1-e^{-\boldsymbol{\lambda} \cdot \boldsymbol{s}}) \, \mathrm{d}\nu(\boldsymbol{s},x),\tag{15}\] where \(\boldsymbol{\lambda}= (\lambda_1,\dots,\lambda_d) \in [0,+\infty)^d, \boldsymbol{s}= (s_1,\dots,s_d) \in \Omega_d\) and \(\cdot\) denotes the scalar product. The multivariate Laplace exponent is defined as \[\label{eq:laplace95exponent95crv} \psi(\boldsymbol{\lambda}) = - \log \mathbb{E}(e^{-\boldsymbol{\lambda} \cdot \boldsymbol{\tilde{\mu}}(\mathbb{X})}) = \iint_{\Omega_d \times \mathbb{X}} (1-e^{-\boldsymbol{\lambda} \cdot \boldsymbol{s}}) \, \mathrm{d}\nu(\boldsymbol{s},x),\tag{16}\] Consistently with the univariate case, a CRV is homogeneous if \(\nu = \rho \otimes P_0\), where \(P_0\) is an atomless measure on \(\mathbb{X}\) and \(\rho\) is \(d\)-dimensional Lévy measure on \(\Omega_d\) such that \[\label{eq:conditions95levy95crv} \text{(i)} \quad \rho_x((\varepsilon,+\infty)^d) <+\infty, \qquad \text{(ii)} \quad \int_{[0,\varepsilon)^d \setminus \{\boldsymbol{0}\}} \|\boldsymbol{s}\| \, \mathrm{d}\rho_x(\boldsymbol{s}) <+\infty.\tag{17}\] Under homogeneity, if \(P_0\) is a probability measure, the multivariate Laplace exponent is \[\psi(\boldsymbol{\lambda}) = \int_{\Omega_d} (1-e^{-\boldsymbol{\lambda} \cdot \boldsymbol{s}}) \, \mathrm{d}\rho(\boldsymbol{s}),\] and \(\log \mathbb{E}(e^{-\boldsymbol{\lambda} \cdot \boldsymbol{\tilde{\mu}}(A)}) = - P_0(A) \, \psi(\boldsymbol{\lambda})\), for \(\boldsymbol{\lambda} \in [0,+\infty)^d\). Finally, if \(\boldsymbol{\tilde{\mu}} \sim {\rm CRV}(\nu)\), it easily follows from the definition that the marginal random measures \(\tilde{\mu}_i\) are CRMs. Their Lévy intensities can be obtained by marginalization of \(\nu\) as \(\mathrm{d}\nu_i(s_i,x) = \int_{ \boldsymbol{s}_{-i} \in [0,+\infty)^{d-1}} \, \mathrm{d}\nu(\boldsymbol{s},x)\), where \(\boldsymbol{s}_{-i} = (s_1,\dots,s_{i-1},s_{i+1},\dots,s_d)\).

9 Proofs↩︎

9.1 Proof of Theorem 1↩︎

Step 1. Show that \(\boldsymbol{\tilde{\mu}}\) is a CRV. Let \(\boldsymbol{\tilde{\mu}}(A) = ( \tilde{\mu}_1(A),\dots, \tilde{\mu}_d(A))\) be the random vector of evaluations of \(\boldsymbol{\tilde{\mu}}\) on the Borel set \(A\) of \(\mathbb{X}\). First, we observe that \(\boldsymbol{\tilde{\mu}}\) takes values in the space of boundedly finite measures. Indeed, by Definition 2, for every bounded Borel set \(A\) and conditionally on \(\tilde{\mu}_0\), \[\tilde{\mu}_i(A)\mid \tilde{\mu}_0 \sim \mathrm{ID}(\tilde{\mu}_0(A)\rho),\] where \(\mathrm{ID}(\eta)\) denotes the infinitely divisible distribution with Lévy measure \(\eta\). Since \(\tilde{\mu}_0\) is a completely random measure, \(\tilde{\mu}_0(A)<+\infty\) almost surely; hence, \(\tilde{\mu}_0(A)\rho\) is a Lévy measure, and \(\tilde{\mu}_i(A)<+\infty\) a.s. for every \(i=1,\ldots,d\). At this point, we prove that the random vectors \(\boldsymbol{\tilde{\mu}}(A_1), \dots, \boldsymbol{\tilde{\mu}}(A_k)\) are mutually independent, for every \(A_1, \cdots, A_k\) mutually disjoint sets of \(\mathbb{X}\). Specifically, we show that all linear combinations are mutually independent, that is, for coefficients \(\lambda_{ij}>0\), with \(i=1,\dots,d\) and \(j=1,\dots, k\), \[\mathbb{E}\Big(e^{- \sum_{j=1}^k \sum_{i=1}^d \lambda_{ij} \tilde{\mu}_i(A_j)}\Big) = \prod_{j=1}^k \mathbb{E}\Big(e^{- \sum_{i=1}^d \lambda_{ij} \tilde{\mu}_i(A_j)}\Big).\] This identity is proved exploiting the following properties:

• \(\tilde{\mu}_1,\dots,\tilde{\mu}_d\) are conditionally independent given \(\tilde{\mu}_0\);

• since \(\tilde{\mu}_i \mid \tilde{\mu}_0\) is a CRM, its evaluations on disjoint sets are independent;

• Since \(\tilde{\mu}_i(A_j) \mid \tilde{\mu}_0 \sim \mathrm{ID}(\tilde{\mu}_0(A_j) \rho)\), then \(\mathcal{L}(\tilde{\mu}_i(A_j) \mid \tilde{\mu}_0) = \mathcal{L}(\tilde{\mu}_i(A_j) \mid \tilde{\mu}_0(A_j))\), where \(\mathcal{L}(X)\) is the probability law of the random object \(X\);

• since \(\tilde{\mu}_0\) is a CRM, \(\tilde{\mu}_0(A_1),\dots, \tilde{\mu}_0(A_k)\) are independent;

• again, \(\tilde{\mu}_1,\dots,\tilde{\mu}_d\) are conditionally independent given \(\tilde{\mu}_0\).

This entails \[\begin{align} \mathbb{E}\Big(e^{- \sum_{j=1}^k \sum_{i=1}^d \lambda_{ij} \tilde{\mu}_i(A_j)}\Big) & \overset{\text{(i)}}{=} \mathbb{E}\bigg(\prod_{i=1}^d \mathbb{E}\Big(e^{- \sum_{j=1}^k \lambda_{ij} \tilde{\mu}_i(A_j)} \mid \tilde{\mu}_0 \Big) \bigg) \\ & \overset{\text{(ii)}}{=} \mathbb{E}\bigg(\prod_{i=1}^d \prod_{j=1}^k \mathbb{E}\Big(e^{- \lambda_{ij} \tilde{\mu}_i(A_j)}\mid \tilde{\mu}_0\Big)\bigg) \\ & \overset{\text{(iii)}}{=} \mathbb{E}\bigg(\prod_{i=1}^d \prod_{j=1}^k \mathbb{E}\Big(e^{- \lambda_{ij} \tilde{\mu}_i(A_j)} \mid \tilde{\mu}_0(A_j)\Big)\bigg) \\ & \overset{\text{(iv)}}{=} \prod_{j=1}^k \mathbb{E}\bigg(\prod_{i=1}^d \mathbb{E}\Big(e^{- \lambda_{ij} \tilde{\mu}_i(A_j)} \mid \tilde{\mu}_0(A_j)\Big)\bigg) \\ & \overset{\text{(v)}}{=} \prod_{j=1}^k \mathbb{E}\Big(e^{- \sum_{i=1}^d \lambda_{ij} \tilde{\mu}_i(A_j)}\Big). \end{align}\]

Step 2. Derive the Laplace transform. From the specification in Definition 2, it follows that \(\tilde{\mu}_i(A) \mid \tilde{\mu}_0 \stackrel{\scriptsize{\rm iid}}{\sim}\mathrm{ID}(\tilde{\mu}_0(A) \rho)\). Therefore, \[\begin{align} \log \mathbb{E}\Big(e^{- \boldsymbol{\lambda} \cdot \boldsymbol{\tilde{\mu}}(A)}\Big) & = \log \mathbb{E}\Big(\mathbb{E}\Big(e^{- \boldsymbol{\lambda} \cdot \boldsymbol{\tilde{\mu}}(A)} \mid \tilde{\mu}_0\Big)\Big) = \log \mathbb{E}\bigg( \prod_{i=1}^d e^{- \psi(\lambda_i) \tilde{\mu}_0(A)} \bigg) \\ & = \log \mathbb{E}\Big( e^{- \tilde{\mu}_0(A) \sum_{i=1}^d \psi(\lambda_i)} \Big) = - P_0(A) \, \psi_0\bigg(\sum_{i=1}^d \psi(\lambda_i)\bigg). \end{align}\]

Step 3. Determine homogeneity and Laplace exponent. The multivariate Lévy intensity \(\nu_h\) and the Laplace exponent \(\psi_h\) of \(\boldsymbol{\tilde{\mu}}\) are uniquely characterized by the Laplace transform. Given the product form of Step 2, it follows that \(\nu_h = \rho_h \otimes P_0\) for some multivariate \(\rho_h\) and the Laplace exponent satisfies \[\psi_h (\boldsymbol{\lambda}) = \psi_0 \bigg(\sum_{i=1}^d \psi(\lambda_i) \bigg).\]

Step 4. Determine the Lévy measure. The Lévy measure \(\rho_h\) is characterized by \[\psi_h (\boldsymbol{\lambda}) = \int_{\Omega_d} (1-e^{-\boldsymbol{\lambda} \cdot \boldsymbol{s}}) \, \mathrm{d}\rho_h (\boldsymbol{s}).\] We use the notation \(\mathscr{L}_{\rho}(\lambda) = \mathbb{E}(e^{-\lambda X}) = e^{-\psi(\lambda)}\), where \(\lambda\in[0,+\infty)\) and \(X \sim \mathrm{ID}(\rho)\), to denote the Laplace transform of an infinitely divisible random variable with Lévy measure \(\rho\) and Laplace exponent \(\psi\); note that \(\mathscr{L}_{\rho}(\lambda)^t = \mathscr{L}_{t\rho}(\lambda)\). Defining \(X_i \stackrel{\scriptsize{\rm iid}}{\sim}\mathrm{ID}(t\rho)\), one has \[\begin{align} \psi_h (\boldsymbol{\lambda}) & = \psi_0 \bigg(\sum_{i=1}^d \psi(\lambda_i) \bigg) \\ & = \int_0^{+\infty} \Big(1-e^{-\sum_{i=1}^d \psi(\lambda_i)\,t}\Big) \, \mathrm{d}\rho_0(t) = \int_0^{+\infty} \bigg(1-\prod_{i=1}^d e^{- \psi(\lambda_i)\,t} \bigg) \, \mathrm{d}\rho_0(t) \\ & = \int_0^{+\infty} \bigg(1-\prod_{i=1}^d \mathscr{L}_{\rho}(\lambda_i)^t \bigg) \, \mathrm{d}\rho_0(t) = \int_0^{+\infty} \bigg(1-\prod_{i=1}^d \mathscr{L}_{t\rho}(\lambda_i) \bigg) \, \mathrm{d}\rho_0(t) \\ & = \int_0^{+\infty} \bigg(1-\prod_{i=1}^d \mathbb{E}_{X_i}\Big(e^{-\lambda_i X_i}\Big)\bigg) \, \mathrm{d}\rho_0(t) = \int_0^{+\infty} \mathbb{E}_{ \boldsymbol{X} }\Big(1-e^{-\sum_{i=1}^d \lambda_i X_i}\Big) \, \mathrm{d}\rho_0(t) \\ & = \int_0^{+\infty} \int_{[0,+\infty)^d} \Big(1-e^{-\sum_{i=1}^d \lambda_i s_i}\Big) \prod_{i=1}^d \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)} (s_i) \, \mathrm{d}\rho_0(t) \\ & = \int_0^{+\infty} \int_{\Omega_d} \Big(1-e^{-\sum_{i=1}^d \lambda_i s_i}\Big) \prod_{i=1}^d \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)} (s_i) \, \mathrm{d}\rho_0(t) \\ & = \int_{\Omega_d} ( 1- e^{-\boldsymbol{\lambda}\cdot \boldsymbol{s}}) \int_0^{+\infty} \prod_{i=1}^d \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)} (s_i) \, \mathrm{d}\rho_0(t) , \end{align}\] by exchanging the order of the integrals thanks to Fubini-Tonelli theorem, and observing that the integrand vanishes in \(\boldsymbol{s} = \boldsymbol{0}\). If \(\mathrm{ID}(t \rho)\) has density \(f_{\scriptscriptstyle \mathrm{ID}(t \rho)}\) with respect to the Lebesgue measure, then this is equal to \[\int_{\Omega_d} ( 1- e^{-\boldsymbol{\lambda}\cdot \boldsymbol{s}}) \int_0^{+\infty} \prod_{i=1}^d f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s_i) \, \mathrm{d}\rho_0(t) \, \mathrm{d}\boldsymbol{s}.\]

9.2 Proof of Remark 2↩︎

Assume that the Lévy measure \(\rho_{\rm co}\) of a compound random measure coincides with \(\rho_h\) in Theorem 1 for some base Lévy measure \(\rho_0\). Since \(\rho_{\rm co}\), \(\rho_h\) and \(\rho_0\) are \(\sigma\)-finite measures, by the uniqueness of the disintegration of measures \(\rho_0\)-a.e. it holds that \[\frac{1}{t^d} \, H \bigg(\frac{s_1}{t},\dots,\frac{s_d}{t}\bigg) \mathrm{d}\boldsymbol{s} = \prod_{i=1}^d \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s_i).\] This forces \(H\) to be in product form and \(\mathrm{ID}(t\rho)\) to have a density. Therefore, there exists a density \(H_1\) on \([0,+\infty)\) such that \[\frac{1}{t} \, H_1\bigg(\frac{s}{t}\bigg) = f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s), \qquad \text{\rho_0-a.e.}\] Let \(X\) be a random variable with density \(H_1\); then \(tX\) has density \(f(s) = H_1(s/t)/t\). From the identity above, \(tX\) is an infinitely divisible distribution with Lévy density \(t \rho(s)\). Therefore, \(X\) has (strictly) stable distribution [41]; in particular, it has stability index \(1\), corresponding to a Lévy density \(\rho(s) \propto s^{-2}\). However, this choice of \(\rho(s)\) does not satisfy the integrability condition (b) in 14 , and thus is not a valid Lévy density for a CRM.

9.3 Proof of Theorem 4↩︎

Step 1. Express the condition in terms of Laplace exponents. Since hierarchical CRVs are CRVs, their law is uniquely determined by their multivariate Laplace exponent. The \(1\)-dimensional Laplace exponents \(\psi_0^{(1)}\) and \(\psi^{(2)}\), associated to \(\rho_0^{(1)}\) and \(\rho^{(2)}\) respectively, are strictly increasing by definition, which implies that they are invertible. Therefore, from the expression of Theorem 1, \(\boldsymbol{\tilde{\mu}}_1\) and \(\boldsymbol{\tilde{\mu}}_2\) have the same multivariate Laplace exponent if and only if \[\sum_{i=1}^d \psi^{(1)} \circ (\psi^{(2)})^{-1} (s_i) = (\psi_0^{(1)})^{-1}\circ \psi_0^{(2)} \bigg(\sum_{i=1}^d s_i \bigg), \qquad P_0^{(1)} = P_0^{(2)},\] for every \((s_1,\dots,s_d) \in \Omega_d\). Define \(f = \psi^{(1)} \circ (\psi^{(2)})^{-1}\) and \(g = (\psi_0^{(1)})^{-1}\circ \psi_0^{(2)}\). Since the inverse of a continuous function is continuous, \(f,g: (0,+\infty) \to (0,+\infty)\) are continuous functions such that \(f(0) = g(0) = 0\) and \[\label{eq:add} f (s_1)+\dots + f(s_d) = g (s_1+\dots + s_d).\tag{18}\]

Step 2. Show that \(f = g\) and are linear functions. Let \(s_1 = s\) and \(s_2 = \dots = s_d = 0\). Then, since \(f(0) = 0\), we have \(g = f\). To prove linearity, we show that \(f(s_1 + s_2) = f(s_1) + f(s_2)\) for every \(s_1,s_2 >0\) and \(f(c s) = c f(s)\) for every \(c,s>0\). The first property follows by taking \(s_3 = \dots = s_d = 0\) in 18 . As for the second, we prove it first for \(c\) a natural number and then for \(c\) a rational number. We conclude by density thanks to the continuity of \(f\). For \(n \in \mathbb{N}\) natural number, by 18 , \(f(ns) = f(s + \dots + s) = n f(s)\), while for \(q = n/m\) rational number, with \(n,m \in \mathbb{N}\), \(nf(s) = f(ns) = f(qms) = m f(qs)\), which implies \(f(qs) = nf(s)/m = q f(s)\).

Step 3. Conclusion. We have proved that there exists \(c>0\) such that, for every \(s>0\), \[\psi^{(1)} \circ (\psi^{(2)})^{-1} (s)= (\psi_0^{(1)})^{-1}\circ \psi_0^{(2)}(s) = c s.\] Thus for every \(s>0\), \[\psi_0^{(2)}(s) = \psi_0^{(1)} (c s), \qquad \psi^{(2)}(s) = \frac{1}{c} \, \psi^{(1)} (s), \qquad P_0^{(1)} = P_0^{(2)},\] Using a change of variable, this condition can be equivalently expressed in terms of the Lévy measures as \[\rho_0^{(2)} = c_{\#} \rho_0^{(1)}; \qquad \rho^{(2)}(s) = \frac{1}{c} \, \rho^{(1)} (s); \qquad P_0^{(1)} = P_0^{(2)}.\] If the Lévy measures have densities, this is equivalent to \[\rho_0^{(2)}(s) = \frac{1}{c} \, \rho_0^{(1)} \bigg(\frac{s}{c} \bigg);\qquad \rho^{(2)}(s) = \frac{1}{c} \, \rho^{(1)} (s); \qquad P_0^{(1)} = P_0^{(2)},\] which coincides with the statement of the proposition by substituting \(c\) with \(1/c\).

9.4 Proof of Lemma 1↩︎

Step 1. Find the Lévy measure \(\rho_1\) of \(\tilde{\mu}_i\). By Theorem 1, the multivariate Lévy measure of \(\boldsymbol{\tilde{\mu}}\) is the measure on \(\Omega_d\) satisfying \[\mathrm{d}\rho_{h}(s_1,\dots, s_d) = \int_0^{+\infty} \prod_{i=1}^d \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)} (s_i) \, \mathrm{d}\rho_0(t).\] The marginal Lévy measure of each \(\tilde{\mu}_i\) coincides with the one of \(\tilde{\mu}_1\), that is, \[\mathrm{d}\rho_1(s) = \int_{[0,+\infty)^{d-1}}\int_0^{+\infty} \mathrm{d}P_{ \scriptscriptstyle \mathrm{ID}(t \rho)}(s) \, \prod_{i=2}^d \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)} (s_i) \, \mathrm{d}\rho_0(t) = \int_0^{+\infty} \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s) \, \mathrm{d}\rho_0(t),\] where we have exchanged the order of integration thanks to Fubini-Tonelli’s theorem, since measures are positive and \(\sigma\)-finite.

Step 2. Show that if \(\rho\) and \(\rho_0\) are infinitely active, then \(\rho_1\) is infinitely active. If \(\rho\) is infinitely active, then \(\mathrm{ID}(t \rho)\) gives zero probability to the origin, thus \[\int_{\Omega_1}\mathrm{d}\rho_1(s) = \int_{\Omega_1}\int_0^{+\infty} \mathrm{d}P_{ \scriptscriptstyle \mathrm{ID}(t \rho)}(s) \, \mathrm{d}\rho_0(t) = \int_0^{+\infty} \, \mathrm{d}\rho_0(t),\] where again we have used Fubini-Tonelli’s theorem to exchange the order of the integrals. Since \(\rho_0\) is infinitely active, the last integral is \(+\infty\).

Step 3. If \(\rho_0\) is finitely active, then \(\rho_1\) is finitely active. Since \(\mathrm{ID}(t \rho)\) is a probability measure on \([0,+\infty)\), its mass on \((0,+\infty)\) is smaller than or equal to 1. Thus, by Fubini-Tonelli, \[\int_{\Omega_1}\mathrm{d}\rho_1(s) = \int_0^{+\infty }\bigg( \int_{\Omega_1} \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)} (s) \bigg) \, \mathrm{d}\rho_0(t) \le \int_0^{+\infty} \, \mathrm{d}\rho_0(t) <+\infty,\] where the last inequality is due to the finite activity of \(\rho_0\).

Step 4: If \(\rho\) is finitely active, then \(\rho_1\) is finitely active. If \(\rho\) is finitely active with total mass \(a\), then \(Y \sim \mathrm{ID}(t \rho)\) has positive mass in 0, namely \(\mathbb{P}(Y=0) = e^{-t a}\). This follows by observing that \(Y\) is a compound Poisson distribution, and \(\mathbb{P}(Y=0)\) coincides with the probability mass in \(0\) of a Poisson random variable with rate \(\lambda = a t\). Thus, by Fubini-Tonelli’s theorem, \[\int_{\Omega_1}\mathrm{d}\rho_1(s) = \int_0^{+\infty} \int_{\Omega_1} \mathrm{d}P_{\scriptscriptstyle \mathrm{ID}(t \rho)} (s) \;\mathrm{d}\rho_0(t) = \int_0^{+\infty} (1-e^{-at}) \, \mathrm{d}\rho_0(t) = \psi_0(at) <+\infty.\]

9.5 Proof of Proposition 6↩︎

The Dirichlet process arises as a normalization of a gamma CRM [14], thus \[\frac{\tilde{\mu}_1}{\tilde{\mu}_1(\mathbb{X})},\dots,\frac{\tilde{\mu}_d}{\tilde{\mu}_d(\mathbb{X})} \Bigm| \tilde{\mu}_0 \, \stackrel{\scriptsize{\rm iid}}{\sim}\, {\rm DP}( \alpha \tilde{\mu}_0) = {\rm DP}\bigg( \alpha \tilde{\mu}_0(\mathbb{X}) \,\frac{\tilde{\mu}_0}{\tilde{\mu}_0(\mathbb{X})}\bigg).\] For a gamma CRM, the total mass \(\tilde{\mu}_0(\mathbb{X})\) is independent from the normalization [77]. Moreover, \(\tilde{\mu}_0(\mathbb{X}) \sim \text{Gamma}(\alpha_0, b_0)\) and thus \(\alpha \tilde{\mu}_0(\mathbb{X}) \sim \text{Gamma}(\alpha_0, b_0/\alpha)\) is independent of \(\tilde{P}_0 = \tilde{\mu}_0 / \tilde{\mu}_0(\mathbb{X})\).

9.6 Proof of Proposition 7↩︎

We prove that, conditionally on \(\tilde{\mu}_0\), it holds \(\tilde{\mu}^{(1)} \stackrel{\scriptsize{d}}{=}c(\tilde{\mu}_0) \, \tilde{\mu}^{(2)}\) for some positive value \(c(\tilde{\mu}_0)\) depending on \(\tilde{\mu}_0\). This implies that the normalizations of \(\tilde{\mu}^{(1)}\) and \(\tilde{\mu}^{(2)}\) are equal in distribution. Conditionally on \(\tilde{\mu}_0\), both \(\tilde{\mu}^{(1)}\) and \(\tilde{\mu}^{(2)}\) are CRMs, and we can express the condition \(\tilde{\mu}^{(1)} \stackrel{\scriptsize{d}}{=}c(\tilde{\mu}_0) \, \tilde{\mu}^{(2)}\) in terms of their conditional Lévy measures \(\rho^{(1)}\) and \(\rho^{(2)}\) as \[\rho^{(2)}(s) \otimes \frac{\tilde{\mu}_0}{\tilde{\mu}_0(\mathbb{X})} \stackrel{\scriptsize{d}}{=}c(\tilde{\mu}_0) \, \rho^{(1)}(c(\tilde{\mu}_0) s) \otimes \tilde{\mu}_0.\] Plugging in the expressions for the stable Lévy measures, we need to find \(c(\tilde{\mu}_0)\) such that \[\frac{\alpha \, \sigma}{\Gamma(1-\sigma)} \frac{1}{s^{\sigma+1}} \otimes \frac{\tilde{\mu}_0}{\tilde{\mu}_0(\mathbb{X})} \stackrel{\scriptsize{d}}{=}\frac{\alpha \, \sigma}{\Gamma(1-\sigma)} \frac{1}{c(\tilde{\mu}_0)^{\sigma}} \frac{1}{s^{\sigma+1}} \otimes \tilde{\mu}_0.\] The proof is concluded by choosing \(c(\tilde{\mu}_0) = \tilde{\mu}_0(\mathbb{X})^{1/\sigma}\).

9.7 Proof of Proposition 9↩︎

Firstly, we state and prove a related result concerning the unnormalized random measures. For this purpose, denote by \(M_r(\rho)\) the \(r\)-th moment of a Lévy measure \(\rho\), \[M_r(\rho) = \int_{[0,\infty)} s^r\, \mathrm{d}\rho(s).\]

We now prove Proposition 9. Similarly to the proof of Proposition 18, we could derive the desired results in two ways: (i) exploiting the hierarchical structure, or (ii) using the fact that we are normalizing a CRV. When focusing on the normalization, exploiting the CRV structure is particularly convenient for deriving the mean and variance, whereas using the hierarchical structure brings to a straightforward calculation for the covariance. The univariate results in [26] show that if \(\tilde{\mu}\sim \mathrm{CRM}(\rho \otimes P_0)\), where \(P_0\) is a probability measure and \(\psi\) denotes its Laplace exponent, the normalization \(\tilde{P} = \tilde{\mu}/\tilde{\mu}(\mathbb{X})\) satisfies \[\mathbb{E}(\tilde{P}(A)) = P_0(A), \qquad \mathrm{Var}(\tilde{P}(A)) = - P_0(A)(1-P_0(A)) \int_0^{+\infty} u e^{-\psi(u)} \psi'' (u) \, \mathrm{d}u.\] Since \(\tilde{P}_i(A)\) is the normalization of \(\tilde{\mu}_i \sim \mathrm{CRM}( \rho_h \otimes P_0)\), this immediately implies that \(\mathbb{E}(\tilde{P}_i(A)) = P_0(A)\). Moreover, by considering the Laplace exponent of \(\tilde{\mu}_i\) in Theorem 1, \[\mathrm{Var}(\tilde{P}_i(A)) = - P_0(A)(1-P_0(A)) \int_0^{+\infty} u e^{-\psi_0(\psi(u))} (\psi_0 \circ \psi)'' (u) \, \mathrm{d}u.\] By the law of total covariance and the conditional independence of the random probability measures \(\tilde{P}_i\) and \(\tilde{P}_j\) for \(i\neq j\), given \(\tilde{\mu}_0\), \[\begin{align} \mathrm{Cov}(\tilde{P}_i(A), \tilde{P}_j(A)) & = \mathrm{Cov}( \mathbb{E}(\tilde{P}_i(A) \mid \tilde{\mu}_0), \mathbb{E}(\tilde{P}_j(A) \mid \tilde{\mu}_0)) \\ & = \mathrm{Var}\bigg(\frac{\tilde{\mu}_0(A)}{\tilde{\mu}_0(\mathbb{X})}\bigg) = - P_0(A)(1-P_0(A)) \int_0^{+\infty} u e^{-\psi_0(u)} \psi_0'' (u) \, \mathrm{d}u. \end{align}\]

9.8 Proof of Theorem 10↩︎

Since the random probabilities \(\tilde{\mu}_i/\tilde{\mu}_i(\mathbb{X})\) are a.s. discrete with random atoms from a continuous distribution, the model for the observations \(\boldsymbol{X}_{1:d} \mid \boldsymbol{\tilde{\mu}}\) is non-dominated. Therefore, we cannot rely on Bayes theorem to find the posterior distribution, and should use an alternative strategy, based on the multivariate Laplace functional. This proof can be seen as the multivariate extension of the proof in [27]. The multivariate Laplace functional \[\mathbb{E}\exp \left( - \sum_{i=1}^d \int f_i \,\mathrm{d}\tilde{\mu}_i \right)\] is defined for any non-negative measurable functions \(f_1,\dots f_d\) and characterizes the law of any vector of random measures. In particular, when \(\boldsymbol{\tilde{\mu}} \sim {\rm CRV}(\nu)\), it satisfies \[\label{laplace95functional95crv} \log \mathbb{E}\exp \left( - \sum_{i=1}^d \int f_i \,\mathrm{d}\tilde{\mu}_i \right) = - \int_{\Omega_d \times \mathbb{X}} (1-e^{- \sum_{i=1}^d s_i f_i(x)}) \, \mathrm{d}\nu(\boldsymbol{s},x).\tag{19}\] Our goal is to find an appropriate expression of \(\mathbb{E}( e^{- \sum_{i=1}^d \int f_i \,\mathrm{d}\tilde{\mu}_i} \mid \boldsymbol{X}_{1:d} = \boldsymbol{x}_{1:d})\), for any non-negative measurable functions \(f_1,\dots f_d\). Throughout the proof, we exploit the following properties of the conditional expectation, which hold for any random variable \(X,Y\) and for any Borel set \(A\) such that \(\mathbb{P}(Y \in A) \ge 0\):

• conditional expectation w.r.t. events: \(\mathbb{E}(X \mid Y \in A) = \mathbb{E}(X \mathbb{1}_A(Y))/\mathbb{P}(Y \in A)\);

• tower property: \(\mathbb{E}(X) = \mathbb{E}( \mathbb{E}(X\mid Y))\) and \(\mathbb{P}(A) = \mathbb{E}(\mathbb{1}_A) = \mathbb{E}( \mathbb{P}(A \mid Y))\).

Step 1. Express the conditional expectation in terms of events. By the dominated convergence theorem and the exchangeability of the observations \[\mathbb{E}\Big( e^{- \sum_{i=1}^d \int f_i \mathrm{d}\tilde{\mu}_i } \mid \boldsymbol{X}_{1:d} = \boldsymbol{x}_{1:d}\Big) = \lim_{\varepsilon \to 0} \, \mathbb{E}\bigg( e^{- \sum_{i=1}^d \int f_i \mathrm{d}\tilde{\mu}_i } \Bigm| \boldsymbol{X}_{1:d} \in \prod_{i=1}^d \prod_{j=1}^k B_{\varepsilon}^{n_{ij}}(x_j^*) \bigg),\] where \(\boldsymbol{x}^* = (x_1^*,\dots,x_k^*)\) are the unique values in the observations, with multiplicities \(n_{i1},\dots,n_{ik}\) for each group \(i = 1, \dots, d\); see Section 4. Here, \(B_{\varepsilon}(x) = \{ \omega: d(x,\omega) \le \varepsilon\}\) denotes the ball of radius \(\varepsilon > 0\) centered in \(x\), and we use \(B^{m}_{\varepsilon}(x) = B_{\varepsilon}(x) \times \dots \times B_{\varepsilon}(x)\) for their \(m\)-cartesian product. Without loss of generality, we always consider \(\varepsilon\) sufficiently small for the balls \(\{B_{\varepsilon}(x_j^*)\}_j\) to be pairwise disjoint.

Step 2. Condition with respect to events. Using property (a), \[\begin{gather} \mathbb{E}\Big( e^{- \sum_{i=1}^d \int f_i \,\mathrm{d}\tilde{\mu}_i } \mid \boldsymbol{X}_{1:d} = \boldsymbol{x}_{1:d}\Big) \\ = \lim_{\varepsilon \to 0} \, \frac{\displaystyle \mathbb{E}\Big(e^{- \sum_{i=1}^d \int f_i \,\mathrm{d}\tilde{\mu}_i } \, \mathbb{1}_{\prod_{i=1}^d \prod_{j=1}^k B_{\varepsilon}^{n_{ij}}(x_j^*)} (\boldsymbol{X}_{1:d}) \Big)}{\mathbb{P}\big(\boldsymbol{X}_{1:d} \in \prod_{i=1}^d \prod_{j=1}^k B_{\varepsilon}^{n_{ij}}(x_j^*)\big)} =: \lim_{\varepsilon \to 0} \frac{N_\varepsilon(\boldsymbol{x}^*)}{D_\varepsilon(\boldsymbol{x}^*)}, \end{gather}\] where we have introduced the notation \(N_\varepsilon(\boldsymbol{x}^*)\) for the numerator and \(D_\varepsilon(\boldsymbol{x}^*)\) for the denominator. In the next steps (3–6) we show that both numerator and denominator decrease at the same speed as \(\varepsilon \to 0\), namely \[\begin{align} N_\varepsilon(\boldsymbol{x}^*) &= C_N \prod_{j=1}^k P_0 (B_{\varepsilon}(x_j^*)) + o\bigg(\prod_{j=1}^k P_0 (B_{\varepsilon}(x_j^*)) \bigg), \\ D_\varepsilon(\boldsymbol{x}^*) &= C_D \prod_{j=1}^k P_0 (B_{\varepsilon}(x_j^*)) + o\bigg(\prod_{j=1}^k P_0 (B_{\varepsilon}(x_j^*)) \bigg), \end{align}\] for some constants \(C_N,C_D >0\), where \(P_0\) is the diffuse base probability of the CRV. It then follows that the limit above coincides with \(C_N/C_D\). Since the denominator is a special case of the numerator when \(f_1=\dots = f_d =0\), we focus on finding \(C_N\) and then specialize the result to \(C_D\).

Step 3. Express the numerator in terms of \(\boldsymbol{\tilde{\mu}}\). By the tower property (b), \[\begin{align} N_\varepsilon(\boldsymbol{x}^*) & = \mathbb{E}\bigg( e^{- \sum_{i=1}^d \int f_i \, \mathrm{d}\tilde{\mu}_i } \, \mathbb{P}\bigg( \boldsymbol{X}_{1:d} \in \prod_{i=1}^d \prod_{j=1}^k B_{\varepsilon}^{n_{ij}}(x_j^*) \Bigm| \boldsymbol{\tilde{\mu}} \bigg) \bigg) \\ & = \mathbb{E}\bigg( e^{- \sum_{i=1}^d \int f_i \, \mathrm{d}\tilde{\mu}_i } \prod_{i=1}^d \prod_{j=1}^k \frac{\tilde{\mu}_i(B_{\varepsilon}(x_j^*))^{n_{ij}}}{\tilde{\mu}(\mathbb{X})^{n_{ij}}} \bigg) \\ & = \mathbb{E}\bigg( \prod_{i=1}^d \frac{1}{\tilde{\mu}_i(\mathbb{X})^{n_i}} \, e^{-\int f_i \, \mathrm{d}\tilde{\mu}_i} \prod_{j=1}^k \tilde{\mu}_i(B_{\varepsilon}(x_j^*))^{n_{ij}} \bigg). \end{align}\]

Step 4. Use the gamma trick to separate the integrand into independent components. Using the density of a gamma with shape \(n_i\) and rate \(\tilde{\mu}_i(\mathbb{X})\), we rewrite \[\frac{1}{\tilde{\mu}_i(\mathbb{X})^{n_i}} = \frac{1}{\Gamma(n_i)} \int_0^{+\infty} u_i^{n_i-1}e^{-\tilde{\mu}_i(\mathbb{X}) u_i} \, \mathrm{d}u_i.\] Henceforth, we adopt the compact notation \(B_j = B_\varepsilon(x_j^*)\) and \(B_0 = \mathbb{X}\setminus \{B_1 \sqcup \dots \sqcup B_k\}\), with the convention \(n_0 = n_{i0} = 0\) for \(i=1,\dots,d\). Using Fubini-Tonelli’s theorem and the independence property of a CRV on disjoint set-wise evaluations, \[\begin{align} N_\varepsilon(\boldsymbol{x}^*) & = \mathbb{E}\bigg( \prod_{i=1}^d \frac{1}{\Gamma(n_i)} \int_0^{+\infty} u_i^{n_i-1} \prod_{j=0}^k e^{- \tilde{\mu}_i(B_j) u_i} \, \mathrm{d}u_i \prod_{j=0}^k e^{-\int_{B_j} f_i \, \mathrm{d}\tilde{\mu}_i} \tilde{\mu}_i (B_j)^{n_{ij}} \bigg) \\ & = \frac{1}{ \prod_{i=1}^d\Gamma(n_i)} \int_{(0,+\infty)^d} \bigg(\prod_{i=1}^d u_i^{-1} \bigg) \, n_\varepsilon(u;\boldsymbol{x}^*) \, \mathrm{d}\boldsymbol{u}, \end{align}\] where \[n_\varepsilon(u;\boldsymbol{x}^*) = \prod_{j=0}^k \mathbb{E}\bigg( \prod_{i=1}^d e^{-\int_{B_j} (f_i(x) +u_i) \, \mathrm{d}\tilde{\mu}_i(x)} (u_i\,\tilde{\mu}_i(B_j))^{n_{ij}} \bigg).\] We now study the asymptotic behaviour of the quantity \(n_\varepsilon(u;\boldsymbol{x}^*)\).

Step 5. Express the integrand \(n_\varepsilon(\boldsymbol{x}^*)\) in terms of the derivative of the multivariate Laplace functional. Let \(\eta_1, \dots, \eta_d \ge 1\) be auxiliary quantities, such that \(n_\varepsilon(u;\boldsymbol{x}^*)\) can be written as \[\begin{align} n_\varepsilon(u;\boldsymbol{x}^*) & = \prod_{j=0}^k \mathbb{E}\bigg( \prod_{i=1}^d \lim_{\eta_i \to 1^+} e^{-\int_{B_j} (f_i(x) + \eta_i u_i) \, \mathrm{d}\tilde{\mu}_i(x)} (u_i\,\tilde{\mu}_i(B_j))^{n_{ij}} \bigg) \\ & = \lim_{\boldsymbol{\eta} \to 1^+} \prod_{j=0}^k \mathbb{E}\bigg( \prod_{i=1}^d e^{-\int_{B_j} (f_i(x) + \eta_i u_i) \, \mathrm{d}\tilde{\mu}_i(x)} (u_i\,\tilde{\mu}_i(B_j))^{n_{ij}} \bigg), \end{align}\] where \(\boldsymbol{\eta} \to 1^+\) is a compact notation for \(\eta_i \to 1^+\) for each \(i = 1, \dots, d\), and we have exchanged limit and expectation by monotone convergence theorem. We observe that, for some \(\ell \in \{1, \dots, d\}\), \[\frac{\partial}{\partial \eta_\ell} \prod_{i=1}^d e^{-\int_{B_j} (f_i(x) + \eta_i u_i) \, \mathrm{d}\tilde{\mu}_i(x)} = -u_\ell\,\tilde{\mu}_\ell(B_j) \prod_{i=1}^d e^{-\int_{B_j} (f_i(x) + \eta_i u_i) \, \mathrm{d}\tilde{\mu}_i(x)}.\] This formula can be applied recursively, using the convention \(d^0/du^0 = \text{Id}\), for \(\text{Id}\) the identity function. Specifically, \[n_\varepsilon(u;\boldsymbol{x}^*) = \lim_{\boldsymbol{\eta} \to 1^+} \prod_{j=0}^k (-1)^{n_{\bullet j}} \, \mathbb{E}\bigg( \frac{\partial^{n_{\bullet j}}}{\partial \eta_1^{n_{1j}}\dots \partial \eta_d^{n_{dj}}} \prod_{i=1}^d e^{-\int_{B_j} (f_i(x) + \eta_i u_i) \, \mathrm{d}\tilde{\mu}_i(x)} \bigg),\] where \(n_{\bullet j} = n_{1j} + \dots + n_{dj}\). Since \(\eta_i \ge 1\) and \(f_i \ge 0\), the derivative is bounded above by the product \[\prod_{i=1}^d (u_i \, \tilde{\mu}_i(B_j))^{n_{ij}} \exp (-u_i \tilde{\mu}_i(B_j)),\] which has finite mean since the exponential decay of \(\exp (- \sum_{i=1}^d u_i \tilde{\mu}_i(B_j))\) is not compromised by the slower polynomial growth of \(\prod_{i=1}^d (u_i \, \tilde{\mu}_i(B_j))^{n_{ij}}\). Therefore, the derivative is uniformly integrable in \(\eta_1, \dots, \eta_d\) and we can exchange derivative and expectation. Remarkably, we need to introduce \(\boldsymbol{\eta}\) and cannot derive the expressions directly with respect to \(\boldsymbol{u}\): otherwise, the derivative would not be uniformly integrable for CRMs with unbounded moment measures, such as the \(\sigma\)-stable CRM. Using the expression 16 of the multivariate Laplace functional of a CRV, we obtain \[\begin{align} n_\varepsilon(u;\boldsymbol{x}^*) & = \lim_{\boldsymbol{\eta} \to 1^+} \prod_{j=0}^k (-1)^{n_{\bullet j}} \, \frac{\partial^{n_{\bullet j}}}{\partial \eta_1^{n_{1j}} \dots \partial \eta_d^{n_{dj}}} \, \mathbb{E}\bigg( e^{- \sum_{i=1}^d \int_{B_j} (f_i(x) + \eta_i u_i) \,\mathrm{d}\tilde{\mu}_i(x)} \bigg) \\ & = \lim_{\boldsymbol{\eta} \to 1^+} \prod_{j=0}^k (-1)^{n_{\bullet j}} \, \frac{\partial ^{n_{\bullet j}}}{\partial \eta_1^{n_{1j}}\dots \partial \eta_d^{n_{dj}}} \, \exp\bigg(- \int_{\Omega_d \times B_j} \big(1-e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x) + \boldsymbol{\eta u})}\big) \, \mathrm{d}\nu(\boldsymbol{s}, x)\bigg), \end{align}\] where we have used the notation \(\boldsymbol{f}(x) = (f_1(x),\dots, f_d(x))\) and \(\boldsymbol{\eta u} = (\eta_1 u_1, \dots, \eta_d u_d)\). For \(j=0,\dots, k\), define the function \[g_j(\boldsymbol{\eta}, \boldsymbol{u}) = - \int_{\Omega_d \times B_j} \big(1-e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x) + \boldsymbol{\eta u})} \big) \, \mathrm{d}\nu(\boldsymbol{s}, x).\]

Step 6. Determine the asymptotic behaviour of the partial derivatives of \(e^{-g_j}\). For \(j=1,\dots, k\), the partial derivatives of \(g_j\) satisfy, as \(\varepsilon \to 0\), \[\begin{align} \frac{\partial ^{n_{\bullet j}}}{\partial \eta_1^{n_{1j}} \dots \partial \eta_d^{n_{dj}}} \,g_j(\boldsymbol{\eta}, \boldsymbol{u})= \\& = (-1)^{n_{\bullet j}} \bigg(\prod_{i=1}^d u_i^{n_{ij}} \bigg) \int_{\Omega_d \times B_j} \prod_{i=1}^d s_i^{n_{ij}} e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x) + \boldsymbol{\eta u}) } \, \mathrm{d}\nu(\boldsymbol{s}, x) \\ & = (-1)^{n_{\bullet j}} \bigg(\prod_{i=1}^d u_i^{n_{ij}} \bigg) \int_{\Omega_d} \prod_{i=1}^d s_i^{n_{ij}} e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x_j^*) + \boldsymbol{\eta u}) } \mathrm{d}\rho_{x_j^*}(\boldsymbol{s}) \, P_0(B_j) + o(P_0(B_j)).\end{align}\] The first line comes from uniform integrability of the derivative, which allows exchanging the integral and derivative operators; the second line follows from Lebesgue differentiation theorem for atomless measures on Polish spaces.

The multivariate Faà di Bruno formula [78] allows to express multiple partial derivatives of a function. We introduce the notation as in [79], who observed that partial derivatives of the type \(\partial^{n_{\bullet}}/\partial \eta_1^{n_{1}}\dots \partial \eta_d^{n_{d}}\) can be treated as maximally mixed partial derivatives of the type \(\partial^{n_{\bullet}}/ \partial v_1 \cdot \dots \partial v_{n_{\bullet}}\) by allowing for ties among the variables, which brings to more compact expressions. In particular, define \(v_1 = \dots = v_{n_{1j}} = \eta_1\), and \(v_{\sum_{i'=1}^{i-1} n_{i'j}+1} = \dots = v_{\sum_{i'=1}^{i} n_{i'j}} = \eta_i\) for \(i=2,\dots,d\), so that, \[\label{eq:faa} (-1)^{n_{\bullet j}} \frac{\partial^{n_{\bullet j}}}{\partial \eta_1^{n_{1j}}\dots \partial \eta_d^{n_{dj}}} \, e^{-g_j(\boldsymbol{\eta}, \boldsymbol{u})} = e^{-g(\boldsymbol{\eta}, \boldsymbol{u})} \sum_{\pi} \prod _{A \in \pi} \frac{\partial^{|A|}}{ \prod_{i \in A} \partial v_i} g_j(\boldsymbol{\eta}, \boldsymbol{u}),\tag{20}\] where the sum is over all partitions \(\pi\) of the numbers \(\{1,\dots,n_{\bullet j} \}\). From this expression, one can derive a formula in terms of \(\eta_1,\dots,\eta_d\) through appropriate combinatorial coefficients, retrieving the one in [78]. However, in our case, the combinatorial formulation above is not necessary, as we are only interested in the asymptotic behaviour as \(\varepsilon \to 0\). From the previous discussion, all terms \(\partial^{|A|}/ (\prod_{i \in A} \partial v_i) g(\boldsymbol{\eta}, \boldsymbol{u})\) are asymptotically equivalent up to a constant. Hence, from 20 , the asymptotically slowest term is the summand corresponding to the partition \(\pi\) with a minimal number of sets, that is, \(\pi = \{1,\dots, n_{\bullet j}\}\). Therefore, as \(\varepsilon \to 0\), 20 reads \[\begin{gather} (-1)^{n_{\bullet j}} \frac{\partial^{n_{\bullet j}}}{\partial \eta_1^{n_{1j}}\dots \partial \eta_d^{n_{dj}}} \, e^{-g_j(\boldsymbol{\eta}, \boldsymbol{u})} = \exp\bigg(- \int_{\Omega_d \times B_j} (1- e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x) + \boldsymbol{\eta u}) }) \,\mathrm{d}\nu(\boldsymbol{s},x) \bigg) \\ \times \bigg(\prod_{i=1}^d u_i^{n_{ij}} \bigg) \int_{\Omega_d} \prod_{i=1}^d s_i^{n_{ij}} e^{- \boldsymbol{s}\cdot (\boldsymbol{f}(x_j^*) + \boldsymbol{\eta u}) } \mathrm{d}\rho_{x_j^*}(\boldsymbol{s}) \, P_0(B_j) + o(P_0(B_j)). \end{gather}\]

Step 7. Determine the asymptotic behaviour of \(n_\varepsilon(\boldsymbol{x}^*)\). Considering the term \(e^{-g_0}\), which does not vanish as \(\varepsilon \to 0\), we obtain that \(n_\varepsilon(u; \boldsymbol{x}^*)\) is asymptotically equal to \[\begin{gather} n_\varepsilon(\boldsymbol{u}; \boldsymbol{x}^*) = \exp\bigg(- \int_{\Omega_d \times \mathbb{X}} (1- e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x) + \boldsymbol{u}) }) \,\mathrm{d}\nu(\boldsymbol{s},x) \bigg) \\ \times \bigg(\prod_{i=1}^d u_i^{n_{i}} \bigg) \prod_{j=1}^k \int_{\Omega_d} \prod_{i=1}^d s_i^{n_{ij}} e^{- \boldsymbol{s}\cdot (\boldsymbol{f}(x_j^*) + \boldsymbol{u}) } \mathrm{d}\rho_{x_j^*}(\boldsymbol{s}) \prod_{j=1}^k P_0(B_j) + o \bigg( \prod_{j=1}^k P_0(B_j) \bigg), \end{gather}\] where we have computed the limit for \(\boldsymbol{\eta} \to 1^+\) under the integration using monotone convergence.

Step 8. Determine the asymptotic behaviour of the numerator \(N_\varepsilon(\boldsymbol{x}^*)\). From the relation between \(n_\varepsilon(u; \boldsymbol{x}^*)\) and \(N_\varepsilon(\boldsymbol{x}^*)\) in Step 4, by monotone convergence theorem, we have that \(N_\varepsilon(\boldsymbol{x}^*) = C_N \prod_{j=1}^kP_0(B_j) + o( \prod_{j=1}^k P_0(B_j))\), where \(C_N\) equals \[\begin{gather} \frac{1}{ \prod_{i=1}^d\Gamma(n_i)} \int_{(0,+\infty)^d} \bigg(\prod_{i=1}^d u_i^{n_{i}-1} \bigg) \, \exp\bigg(- \int_{\Omega_d \times \mathbb{X}} (1- e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x) + \boldsymbol{u}) }) \,\mathrm{d}\nu(\boldsymbol{s},x) \bigg) \\ \times \prod_{j=1}^k \bigg( \int_{\Omega_d} \prod_{i=1}^d s_i^{n_{ij}} e^{- \boldsymbol{s}\cdot (\boldsymbol{f}(x_j^*) + \boldsymbol{u}) } \mathrm{d}\rho_{x_j^*}(\boldsymbol{s}) \bigg) \, \mathrm{d}\boldsymbol{u}. \end{gather}\]

Step 9. Determine the expression of the Laplace functional a posteriori. By specializing the formula in Step 8 for \(C_N\) to \(f = 0\), we determine the value of \(C_D\), which can be conveniently expressed in terms of the multivariate Laplace exponent \(\psi\) and cumulants \(\tau_{n_{1j},\dots, n_{dj} \mid x_j^*}\) as \[C_D = \frac{1}{ \prod_{i=1}^d\Gamma(n_i)} \int_{(0,+\infty)^d} \bigg(\prod_{i=1}^d u_i^{n_{i}-1} \bigg) \, e^{- \psi(\boldsymbol{u})} \prod_{j=1}^k \tau_{n_{1j},\dots, n_{dj}|x_j^*}(\boldsymbol{u}) \, \mathrm{d}\boldsymbol{u}.\] It follows that the Laplace functional a posteriori is equal to \[\begin{gather} \qquad \mathbb{E}\Big( e^{- \sum_{i=1}^d \int f_i \mathrm{d}\tilde{\mu}_i} \mid \boldsymbol{X}_{1:d} = \boldsymbol{x}_{1:d}\Big) \\ = \frac{ \int_{(0,+\infty)^d} \big(\prod_{i=1}^d u_i^{n_{i}-1} \big) \, e^{- \int_\mathbb{X}\psi(\boldsymbol{f}(x) + \boldsymbol{u})\,\mathrm{d}P_0(x)} \, \mathrm{d}\boldsymbol{u}}{ \int_{(0,+\infty)^d} \big(\prod_{i=1}^d u_i^{n_{i}-1} \big) \, e^{- \psi(\boldsymbol{u})} \prod_{j=1}^k \tau_{n_{1j},\dots, n_{dj}|x_j^*}(\boldsymbol{u}) \, \mathrm{d}\boldsymbol{u}}. \qquad \end{gather}\]

Step 10. Interpret the denominator as the normalizing constant of latent random variables \(U_1,\dots, U_d\). The multivariate Laplace exponent \(\psi(\boldsymbol{u})\) can be retrieve at the numerator by multiplying and dividing by \(e^{-\psi(\boldsymbol{u})}\). Indeed, \[\begin{align} \int_{\Omega_d \times \mathbb{X}} (1-e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x) + \boldsymbol{u}) }) \, \mathrm{d}\nu(\boldsymbol{s}, x) + \psi(\boldsymbol{u}) & = \int_{\Omega_d \times \mathbb{X}} (e^{- \boldsymbol{s} \cdot \boldsymbol{u}} - e^{-\boldsymbol{s} \cdot (\boldsymbol{f}(x) + \boldsymbol{u}) }) \, \mathrm{d}\nu(\boldsymbol{s}, x) \\ & = \int_{\Omega_d \times \mathbb{X}} (1-e^{- \boldsymbol{s} \boldsymbol{f}(x) }) e^{-\boldsymbol{s} \cdot \boldsymbol{u}} \, \mathrm{d}\nu(\boldsymbol{s}, x). \end{align}\] Therefore, the posterior Laplace functional is rewritten as \[\begin{gather} \mathbb{E}\Big( e^{- \sum_{i=1}^d \int f_i \mathrm{d}\tilde{\mu}_i} \mid \boldsymbol{X}_{1:d} = \boldsymbol{x}_{1:d}\Big) \\ = \mathbb{E}\bigg( e^{- \int_{\Omega_d \times \mathbb{X}} (1-e^{- \boldsymbol{s} \boldsymbol{f}(x) }) e^{-\boldsymbol{s} \cdot \boldsymbol{U}} \mathrm{d}\nu(\boldsymbol{s}, x)} \prod_{j=1}^k \frac{ \int_{\Omega_d} \prod_{i=1}^d s_i^{n_{ij}} e^{- \boldsymbol{s} \cdot (\boldsymbol{f}(x_j^*) + \boldsymbol{U}) } \mathrm{d}\rho_{x_j^*}(\boldsymbol{s})}{\tau_{n_{1j},\dots, n_{dj}|x_j^*}(\boldsymbol{U})} \bigg), \end{gather}\] where \(\boldsymbol{U} = (U_1,\dots,U_d)\) is a vector of random variables with joint p.d.f. \[f_{\scriptscriptstyle \boldsymbol{U}}(\boldsymbol{u}) \propto \prod_{i=1}^d u_i^{n_i-1} e^{- \psi(\boldsymbol{u})} \prod_{j=1}^k \tau_{n_{1j},\dots, n_{dj}|x_j^*}(\boldsymbol{u}).\] This implies that there exist latent variables \(\boldsymbol{U} = (U_1,\dots,U_d)\) such that the posterior Laplace functional of \(\boldsymbol{\tilde{\mu}}\), conditionally on \(\boldsymbol{U}\), satisfies \[\mathbb{E}\Big( e^{- \sum_{i=1}^d \int f_i \mathrm{d}\tilde{\mu}_i } \mid \boldsymbol{X}_{1:d} = \boldsymbol{x}_{1:d}, \boldsymbol{U} \Big) = \mathbb{E}\Big(e^{- \sum_{i=1}^d \int f_i \mathrm{d}\tilde{\mu}_i^* } \mid \boldsymbol{U} \Big) \prod_{j=1}^k\mathbb{E}\Big(e^{-\boldsymbol{f}(x_j^*) \cdot \boldsymbol{J}_j} \mid \boldsymbol{U}\Big),\] where, conditionally on \(\boldsymbol{U}\), the measure \(\boldsymbol{\tilde{\mu}}^*\) is a CRV with Lévy intensity \(e^{-\boldsymbol{s} \cdot \boldsymbol{U}} \mathrm{d}\nu(\boldsymbol{s}, x)\) and \(\boldsymbol{J}_j\) is a vector of jumps with distribution \[\mathrm{d}P_{\boldsymbol{J}_j \mid \boldsymbol{U}}(\boldsymbol{s}) \propto \prod_{i=1}^d s_i^{n_{ij}} e^{- \boldsymbol{s}\cdot \boldsymbol{U} } \mathrm{d}\rho_{x_j^*}(\boldsymbol{s}).\] If we define the random elements \(\boldsymbol{\tilde{\mu}}^*\) and \(\boldsymbol{J}_j\) to be conditionally independent given \(\boldsymbol{U}\), the Laplace functional of their sum is the product of their Laplace functions, and thus \[\mathbb{E}\Big( e^{- \sum_{i=1}^d \int f_i \mathrm{d}\tilde{\mu}_i } \mid \boldsymbol{X}_{1:d} = \boldsymbol{x}_{1:d}, \boldsymbol{U} \Big) = \mathbb{E}\Big( e^{ - \sum_{i=1}^d \int f_i \mathrm{d}(\tilde{\mu}_i^* + \sum_{j=1}^{k} J_{ij} \delta_{x_j^*})} \mid \boldsymbol{U} \Big).\] By uniqueness of the Laplace functional, this implies that \[\mathcal{L}( \boldsymbol{\tilde{\mu}} \mid \boldsymbol{X}_{1:d} = \boldsymbol{x}_{1:d}) = \mathcal{L}\bigg(\boldsymbol{\tilde{\mu}}^* + \sum_{j=1}^k \boldsymbol{J}_j \delta_{x^*_j}\bigg).\]

9.9 Proof of Proposition 11↩︎

We first prove a preliminary Lemma on the exponential tilting of a Lévy measure.

Considering \(\boldsymbol{\tilde{\mu}} \sim {\rm hCRV}(\rho, \rho_0, P_0)\), the expression of the Lévy intensity in Theorem 10(i) is given by \[\mathrm{d}\nu^*_{\scriptscriptstyle U}(\boldsymbol{s}, x) = e^{-\boldsymbol{U} \cdot \boldsymbol{s}} \mathrm{d}\nu_h(\boldsymbol{s},x) = e^{-\boldsymbol{U} \cdot \boldsymbol{s}} \rho_h(\boldsymbol{s}) \, \mathrm{d}\boldsymbol{s} \, \mathrm{d}P_0(x),\] where we have used that ID\((t\rho)\) has a p.d.f. on \((0,+\infty)\). Exploiting Lemma 3 for the exponential tilting, we obtain \[\begin{align} e^{-\boldsymbol{U} \cdot \boldsymbol{s}} \rho_h(\boldsymbol{s}) & = \int_0^{+\infty} \prod_{i=1}^d \big(e^{-U_i s_i} f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s_i)\big) \, \rho_0(t)\, \mathrm{d}t \\ & = \int_0^{+\infty} \prod_{i=1}^d \big( e^{-t \psi(U_i)} f_{\scriptscriptstyle \mathrm{ID}(t \rho_{U_i})}(s_i) \big) \, \rho_0(t) \, \mathrm{d}t \\ & = \int_0^{+\infty} \prod_{i=1}^d f_{\scriptscriptstyle \mathrm{ID}(t \rho_{U_i} )}(s_i) \, e^{-t \sum_{i=1}^d \psi(U_i)} \rho_0(t) \, \mathrm{d}t. \end{align}\] In analogy with Theorem 1, this expression can be interpreted as the Laplace exponent of a hierarchical CRV with heterogeneous marginal distributions, characterized by the Lévy measures \[\mathrm{d}\rho_i^*(s) = \mathrm{d}\rho_{\scriptscriptstyle U_i}(s) = e^{-U_i\,s} \rho(s) \,\mathrm{d}s, \qquad \mathrm{d}\rho_0^*(t) = e^{-t \sum_{i=1}^d \psi(U_i)} \rho_0(t) \,\mathrm{d}t.\]

9.10 Proof of Proposition 12↩︎

Since \(\rho\) and \(\rho_0\) have Lévy densities, the distribution of the vector of jumps \(\boldsymbol{J}_j\) in 7 , given the latent variables \(\boldsymbol{U}\), has p.d.f. proportional to \[\begin{align} \prod_{i=1}^d s_i^{n_{ij}} \, e^{-\boldsymbol{U} \cdot \boldsymbol{s}} \rho_{h}(\boldsymbol{s}) & = \prod_{i=1}^d s_i^{n_{ij}} e^{-U_i\,s_i} \int_0^{+\infty} \prod_{i=1}^d f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s_i) \, \rho_0(t) \, \mathrm{d}t \\ & = \int_0^{+\infty} \prod_{i=1}^d s_i^{n_{ij}} e^{-U_i s_i} f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s_i) \, \rho_0(t) \, \mathrm{d}t \\ & = \int_0^{+\infty} \prod_{i=1}^d \frac{s_i^{n_{ij}} e^{-U_i s_i} f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s_i)}{\bar \tau_{n_{ij}} (U_i, t)} \prod_{i=1}^d \bar \tau_{n_{ij}} (U_i, t) \, \rho_0(t) \, \mathrm{d}t, \end{align}\] where \(\bar \tau_{m} (u, t)\) is defined in 8 . Therefore, the distribution of \(\boldsymbol{J}_j \mid \boldsymbol{U}\) is a mixture of conditionally independent random variables \(J_{1j}, \dots, J_{dj}\) with densities \[f_{\scriptscriptstyle J_{ij} \mid U_i, J_{0j}}(s) = \frac{s^{n_{ij}} e^{-U_i s} f_{\scriptscriptstyle \mathrm{ID}(J_{0j} \rho)}(s)}{\bar \tau_{n_{ij}} (U_i, J_{0j})}, \qquad i = 1, \dots, d,\] given a mixing random variable \(J_{0j}\) having p.d.f. proportional to \(\prod_{i=1}^d \bar \tau_{n_{ij}} (U_i, t) \, \rho_0(t)\).

9.11 Proof of Lemma 2↩︎

Substituting the expression for \(\rho_h\) as in the proof of Proposition 12, we use Fubini-Tonelli Theorem as in Remark 5 to obtain \[\begin{align} \tau_{\boldsymbol{m}}(\boldsymbol{u}) & = \int_{\Omega_d} \prod_{i=1}^d s_i^{m_i} e^{-\boldsymbol{u} \cdot \boldsymbol{s}} \rho_{h}(\boldsymbol{s}) \, \mathrm{d}{\boldsymbol{s}} = \int_0^{+\infty} \int_{\Omega_d} \prod_{i=1}^d s_i^{m_i} e^{-u_i s_i} f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s_i) \, \mathrm{d}{\boldsymbol{s}} \, \rho_0(t) \, \mathrm{d}t \\ & = \int_0^{+\infty} \prod_{i=1}^d \int_0^{+\infty} s_i^{m_i} e^{-u_i s_i} f_{\scriptscriptstyle \mathrm{ID}(t \rho)}(s_i) \, \mathrm{d}s_i \, \rho_0(t) \, \mathrm{d}t = \int_0^{+\infty} \prod_{i=1}^d \bar \tau_{m_i} (u_i, t) \rho_0(t) \, \mathrm{d}t. \end{align}\]

9.12 Proof of Proposition 13↩︎

Starting from the expression of the distribution of the latent variables \(\boldsymbol{U}\) in 6 , specialized for \(\boldsymbol{\tilde{\mu}} \sim {\rm hCRV}(\rho, \rho_0, P_0)\), \[\begin{align} f_{\scriptscriptstyle \boldsymbol{U}}(\boldsymbol{u}) & \propto \prod_{i=1}^d u_i^{n_i-1} e^{- \psi_0(\sum_{i=1}^d \psi(u_i))} \prod_{j=1}^k \tau_{n_{1j},\dots, n_{dj}}(\boldsymbol{u}) \\ & = \prod_{i=1}^d u_i^{n_i-1} e^{- \psi_0(\sum_{i=1}^d \psi(u_i))} \prod_{j=1}^k \int_{0}^{+\infty} \prod_{i=1}^d \bar \tau_{m_i} (u_i, t_j) \rho_0(t_j) \, \mathrm{d}t_j, \end{align}\] where we have substituted the expression of \(\tau_{\boldsymbol{m}}(\boldsymbol{u})\) obtained in Lemma 2. For each \(j = 1, \dots, k\), exploiting the definition of \(\bar \tau_{m} (u, t)\) in 8 , we obtain \[\begin{align} \int_{0}^{+\infty} \prod_{i=1}^d \bar \tau_{m_i} (u_i, t_j) \rho_0(t_j) \, \mathrm{d}t_j & = \int_{0}^{+\infty} \prod_{i=1}^d \int_0^\infty s_{ij}^{n_{ij}} e^{-u_i\,s_{ij}} f_{\scriptscriptstyle {\rm ID}(t_j\rho)}(s_{ij}) \, \mathrm{d}s_{ij} \, \rho_0(t_j) \, \mathrm{d}t_j \\ & = \int_{\Omega_d} \prod_{i=1}^d e^{-u_i\,s_{ij}} \int_{0}^{+\infty} \prod_{i=1}^d s_{ij}^{n_{ij}} f_{\scriptscriptstyle {\rm ID}(t_j\rho)}(s_{ij}) \, \rho_0(t_j) \, \mathrm{d}t_j \, \mathrm{d}\boldsymbol{s}_{\bullet j}, \end{align}\] where we have exchanged the integrals thanks to Fubini-Tonelli’s theorem. Recall that, conditionally on \(\tilde{\mu}_0\), the distribution of \(\tilde{\mu}_i(\mathbb{X}) \sim {\rm ID}(\tilde{\mu}_0(\mathbb{X}) \rho)\), while at the root of the hierarchy \(\tilde{\mu}_0(\mathbb{X}) \sim {\rm ID}(\rho_0)\). Therefore, exploiting the definition of multivariate Laplace exponent, \[\begin{align} \prod_{i=1}^d u_i^{n_i-1} e^{- \psi_0(\sum_{i=1}^d \psi(u_i))} & = \prod_{i=1}^d u_i^{n_i-1} \mathbb{E}\big( e^{- {\boldsymbol{u}} \cdot {\boldsymbol{\tilde{\mu}}}(\mathbb{X})} \big) \\ & = \int_{\Omega_d} \prod_{i=1}^d u_i^{n_i-1} e^{- u_i y_i} \int_{0}^\infty \prod_{i=1}^d f_{\scriptscriptstyle {\rm ID}(z \rho)}(y_i) \, f_{\scriptscriptstyle {\rm ID}(\rho_0)}(z) \, \mathrm{d}z. \end{align}\] In conclusion, the density of the vector of latent variables \(\boldsymbol{U}\) is proportional to \[\begin{gather} f_{\boldsymbol{U}}(\boldsymbol{u}) \propto \int_{\Omega_d^{k+1}} \prod_{i=1}^d u_i^{n_i-1} e^{- u_i \, (s_{i0} + s_{i1} + \cdots + s_{ik})} \\ \times \int_{(0,+\infty)^{k+1}} \prod_{i=1}^d f_{\scriptscriptstyle {\rm ID}(t_0 \rho)}(s_{i0}) \, f_{\scriptscriptstyle {\rm ID}(\rho_0)}(t_0) \prod_{j=1}^k s_{ij}^{n_{ij}} f_{\scriptscriptstyle {\rm ID}(t_j\rho)}(s_{ij}) \, \rho_0(t_j) \, \mathrm{d}\boldsymbol{t} \, \mathrm{d}\boldsymbol{s}, \end{gather}\] where we have renamed integration variables \(y_i\) as \(s_{i0}\) and \(z\) as \(t_0\). Therefore, conditionally on a vector \(\boldsymbol{\beta} = (\beta_1, \dots, \beta_d)\) of dependent random variables, the latent variables \(U_1, \dots, U_d\) are independent and gamma distributed, with \(U_i \sim {\rm Gamma}(n_i, \beta_i)\), for \(i = 1, \dots, d\). Moreover, each \(\beta_i = S_{i0} + S_{i1} + \cdots + S_{ik}\), where the joint density of \(\boldsymbol{S} = (S_{ij})_{ij}\) is proportional to \[\begin{gather} f_{\scriptscriptstyle \boldsymbol{S}}(\boldsymbol{s}) \propto \prod_{i=1}^d (s_{i0} + s_{i1} + \cdots + s_{ik})^{-n_i} \\ \times \int_{(0,+\infty)^{k+1}} \prod_{i=1}^d f_{\scriptscriptstyle {\rm ID}(t_0 \rho)}(s_{i0}) \, f_{\scriptscriptstyle {\rm ID}(\rho_0)}(t_0) \prod_{j=1}^k s_{ij}^{n_{ij}} f_{\scriptscriptstyle {\rm ID}(t_j \rho)}(s_{ij}) \, \rho_0(t_j) \, \mathrm{d}\boldsymbol{t} \, \mathrm{d}\boldsymbol{s}. \end{gather}\] This entails that, conditionally on \(\boldsymbol{T} = (T_0, \dots, T_k)\), \(\boldsymbol{S}_i = (S_{i0}, S_{i1}, \cdots, S_{ik})\) are independent for each \(i = 1, \dots, d\), with density proportional to \[f_{\scriptscriptstyle \boldsymbol{S}_i \mid \boldsymbol{T}}(\boldsymbol{s}_i) \propto (s_{i0} + s_{i1} + \cdots + s_{ik})^{-n_i} f_{\scriptscriptstyle {\rm ID}(T_0 \rho)}(s_{i0})\, \prod_{j=1}^k s_{ij}^{n_{ij}} f_{\scriptscriptstyle {\rm ID}(T_j \rho)}(s_{ij}).\] Finally, the density of vector \(\boldsymbol{T}\) is proportional to \[f_{\scriptscriptstyle \boldsymbol{T}}(\boldsymbol{t}) \propto \prod_{i=1}^d C(n_{i1}, \dots, n_{ik}; \boldsymbol{t}) \, f_{\scriptscriptstyle {\rm ID}(\rho_0)}(t_0) \prod_{j=1}^k \rho_0(t_j),\] where \(C\), defined in 9 , represents the normalizing constant for the distribution of \(\boldsymbol{S}_i \mid \boldsymbol{T}\).

9.13 Proof of Proposition 14↩︎

From Example 1, the gamma-gamma hCRV is characterized by \[\rho(s) = \alpha \, \frac{e^{-bs}}{s}, \qquad \rho_0(s) = \alpha_0 \, \frac{e^{-b_0s}}{s},\] where \(\alpha,\alpha_0>0\) are shape parameters and \(b,b_0>0\) rate parameters. The Laplace exponent in Definition 5 is \(\psi(\lambda) = \alpha \log (1+\lambda/b)\). The rest of the proof follows from Propositions 11 and 12.

• From Proposition 11, and substituting the expressions for \(\rho\), \(\rho_0\) and \(\psi\), \[\begin{gather} e^{-U_i s} \rho(s) = \alpha s^{-1} e^{-bs-U_i s} = \alpha s^{-1} e^{-b(1+U_i/b)s}, \\ e^{-\sum_{i=1}^d \psi(U_i)s} \rho_0(s) = \alpha_0 s^{-1} e^{-b_0s-\alpha \sum_{i=1}^d \log (1+U_i/b) s} = \alpha_0 s^{-1} e^{-\alpha \lambda(\boldsymbol{U})s}, \end{gather}\] where \(\displaystyle \lambda(\boldsymbol{U}) = b_0/\alpha + \sum_{i=1}^d \log (1+U_i/b)\).

• From Proposition 12, for each \(j=1,\dots,k\), the jumps \(J_{1j},\dots,J_{dj}\) are conditionally independent, given \(\boldsymbol{U}\) and \(J_{0j}\). Moreover, from the specification of \(\rho\) above, the random variable \(\text{ID}(t\rho)\) has gamma distribution with shape parameter \(\alpha t\) and rate parameter \(b\). Hence, for each \(i=1,\dots,d\), the jump \(J_{ij}\) has density proportional to \[s^{n_{ij}} e^{-U_is} f_{\scriptscriptstyle \text{ID}(J_{0j} \rho)}(s) \propto s^{\alpha J_{0j} + n_{ij} - 1} e^{-bs - U_i s},\] which is the density of a gamma random variable with shape \(\alpha J_{0j} + n_{ij}\) and rate \(b+U_i\). The normalizing constant \(\bar \tau_{n_{ij}}(U_i,t)\) in 8 is \[\bar \tau_{n_{ij}}(U_i,t) = \int_0^{+\infty} s^{n_{ij}} e^{-U_is} f_{\scriptscriptstyle \text{ID}(t \rho)}(s) \, \mathrm{d}s = \frac{b^{\alpha t}}{(b + U_i)^{n_{ij} + \alpha t}} \frac{\Gamma(n_{ij} + \alpha t)}{\Gamma(\alpha t)}.\]

• Again from Proposition 12, for \(j=1,\dots,k\), the density of \(J_{0j}\), given \(\boldsymbol{U}\), is proportional to \[\begin{align} f_{\scriptscriptstyle J_{0j} \mid \boldsymbol{U}}(t) & \propto \prod_{i=1}^d \bar \tau_{n_{ij}}(U_i,t) \,\rho_0(t) \propto \prod_{i=1}^d \left( \left(\frac{1}{1 + U_i/b} \right)^{\alpha t} ((\alpha t))_{n_{ij}} \right) \alpha_0 \, \frac{e^{-b_0\,t}}{t} \, \mathrm{d}t \\ & \propto t^{-1} \, e^{-b_0 t - \alpha t \sum_{i=1}^d \log(1 + U_i/b)} \prod_{i=1}^d ((\alpha t))_{n_{ij}} \propto t^{-1} \, e^{-\alpha \lambda(\boldsymbol{U}) t} \prod_{i=1}^d ((\alpha t))_{n_{ij}}, \end{align}\] where \(((\alpha t))_n = \Gamma(\alpha t+n)/\Gamma(\alpha t)\) denotes the ascending factorial. The result is obtained computing the density of the linear transformation \(\alpha J_{0j}\).

9.14 Proof of Proposition 15↩︎

Recall that, from the specification of \(\rho\) in Example 1, the random variable \(\text{ID}(t\rho)\) has gamma distribution with shape parameter \(\alpha t\) and rate parameter \(b\). The rest of the proof follows from Proposition 13.

• For each \(i=1,\dots,d\), the density of \(\boldsymbol{S}_i = (S_{i0},\dots, S_{ik})\), given \(\boldsymbol{T} = (T_0,\dots,T_k)\), is proportional to \[\begin{align} f_{\scriptscriptstyle \boldsymbol{S}_i \mid {\boldsymbol{T}}}(s_{i0},\dots,s_{ik}) & \propto (s_{i0} + s_{i1} + \dots + s_{ik})^{-n_i} f_{\scriptscriptstyle \mathrm{ID}(T_0 \rho)}(s_{i0}) \prod_{j=1}^k s_{ij}^{n_{ij}} f_{\scriptscriptstyle \mathrm{ID}(T_j \rho)}(s_{ij}) \\ & \propto (s_{i0} + s_{i1} + \dots + s_{ik})^{-n_i} s_{i0}^{\alpha T_0-1} e^{-b s_{i0}} \prod_{j=1}^k s_{ij}^{\alpha T_j + n_{ij}-1} e^{-b s_{ij}} \\ & \propto (s_{i0} + s_{i1} + \dots + s_{ik})^{-n_i} e^{-b (s_{i0} + s_{i1} + \dots + s_{ik})} s_{i0}^{\alpha T_0-1} \prod_{j=1}^k s_{ij}^{\alpha T_j + n_{ij}-1}. \end{align}\] Applying the change of variables \(\beta_i = S_{i0}+\cdots+S_{ik}\) and \(W_{ij} = S_{ij}/\beta_i\), for \(j=0,\dots,k\), the joint density of \(\beta_i \ge 0\) and \(\boldsymbol{W}_i = (W_{i0},\dots,W_{ik}) \in \Delta^k\), where \(\Delta^k\) is the \(k\)-dimensional unit simplex, is \[f_{\scriptscriptstyle \beta_i,\boldsymbol{W}_i\mid {\boldsymbol{T}}}(z_i,\boldsymbol{w}_i) \propto z_i^{\alpha (T_0+T_1+\cdots+T_k)-1} e^{-b z_i} \, w_{i0}^{\alpha T_0-1} \prod_{j=1}^k w_{ij}^{\alpha T_j + n_{ij}-1}.\] Therefore \(\beta_i = S_{i0}+\cdots+S_{ik}\) is independent from \(\boldsymbol{W}_i\) and has gamma distribution with shape \(\alpha (T_0+\cdots+T_k)\) and rate \(b\). Moreover, the quantity \(C(\boldsymbol{m}; \boldsymbol{t})\) in 9 is given by \[\begin{align} C(\boldsymbol{m}; & \boldsymbol{t}) = \int_{(0,+\infty)^{k+1}} (s_{0} + \dots + s_{k})^{-m_{\bullet}} f_{\scriptscriptstyle \mathrm{ID}(t_0 \rho)}(s_{0}) \prod_{j=1}^k s_{j}^{m_{j}} f_{\scriptscriptstyle \mathrm{ID}(t_j \rho)}(s_{j}) \mathrm{d}\boldsymbol{s} \\ & = \frac{b^{\alpha \sum_{j=0}^k t_j}}{\Gamma(\alpha t_0) \prod_{j=1}^k \Gamma(\alpha t_j)} \int_{(0,+\infty)} z^{\alpha \sum_{j=0}^k t_j -1} e^{-b z} \mathrm{d}z \, \int_{\Delta^k} w_{0}^{\alpha t_0-1} \prod_{j=1}^k w_{j}^{\alpha t_j + m_{j}-1} \mathrm{d}\boldsymbol{w} \\ & = \frac{\Gamma(\alpha (t_0+\cdots+t_k))}{\Gamma(\alpha (t_0+\cdots+t_k) + m_{\bullet})} \prod_{j=1}^k \frac{\Gamma(\alpha t_j+m_j)}{\Gamma(\alpha t_j)}. \end{align}\]

• From the specification of \(\rho_0\) in Example 1, the random variable \(\text{ID}(\rho_0)\) has gamma distribution with shape parameter \(\alpha_0\) and rate parameter \(b_0\). Therefore, the density of \(\boldsymbol{T} = (T_0,\dots,T_k)\) is proportional to \[\begin{align} f_{\scriptscriptstyle {\boldsymbol{T}}}(\boldsymbol{t}) & \propto \prod_{i=1}^d C(n_{i1},\dots,n_{ik}; \boldsymbol{t}) f_{\scriptscriptstyle \mathrm{ID}(\rho_0)}(t_0) \prod_{j=1}^k \rho_0(t_j) \\ & \propto \prod_{i=1}^d \left( \frac{\Gamma(\alpha (t_0+\cdots+t_k))}{\Gamma(\alpha (t_0+\cdots+t_k) + n_i)} \prod_{j=1}^k \frac{\Gamma(\alpha t_j+n_{ij})}{\Gamma(\alpha t_j)} \right) t_0^{\alpha_0-1} \, e^{-b_0 t_0} \prod_{j=1}^k t_j^{-1} \, e^{-b_0 t_j} \\ & \propto \prod_{i=1}^d \left( \frac{1}{((\alpha(t_0+\cdots+t_k)))_{n_i}} \right) t_0^{\alpha_0-1} \, e^{-b_0 t_0} \prod_{j=1}^k t_j^{-1} \, e^{-b_0 t_j} \left( \prod_{i=1}^d ((\alpha t_j))_{n_{ij}} \right), \end{align}\] where \(((\alpha t))_n = \Gamma(\alpha t+n)/\Gamma(\alpha t)\) is the ascending factorial. Since we only need the distribution of \(\alpha T = \alpha(T_0 +\dots + T_k)\) to sample the random variables \(\beta_1, \dots, \beta_d\), we could apply the change of variables \[\alpha T = \alpha(T_0 +\dots + T_k), \qquad V_{j} = T_{j}/T, \quad j = 0,\dots,k,\] and obtain the joint density of \(\alpha T\) and the vector \(\boldsymbol{V} = (V_0,\dots,V_k) \in \Delta^k\) of auxiliary latent variables, supported on the \(k\)-dimensional unit simplex \(\Delta^k\): \[f_{\scriptscriptstyle \alpha T,\boldsymbol{V}}(t,\boldsymbol{v}) \propto t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} \, \prod_{i=1}^d \frac{1}{((t))_{n_i}} \;v_0^{\alpha_0-1} \prod_{j=1}^k \left( v_j^{-1}\prod_{i=1}^d ((t v_j))_{n_{ij}} \right).\] Therefore, the marginal density of \(\alpha T\) is \[f_{\scriptscriptstyle \alpha T}(t) \propto t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} \, \prod_{i=1}^d \frac{1}{((t))_{n_i}} \int_{\Delta^k} v_0^{\alpha_0-1} \, \prod_{j=1}^k \left( v_j^{-1}\prod_{i=1}^d ((t v_j))_{n_{ij}} \right) \mathrm{d}\boldsymbol{v}.\] Note that this function is integrable in \(t\); indeed, it holds that \(((s))_q \sim s\) for \(s \to 0\) when \(q\) is a positive integer. Therefore, for \(t \to 0\), \[\begin{align} f_{\scriptscriptstyle \alpha T}(t) & \sim t^{\alpha_0-1} \prod_{i=1}^d \frac{1}{((t))_{n_i}} \int_{\Delta^k} v_0^{\alpha_0-1} \, \prod_{j=1}^k \left( v_j^{-1}\prod_{i=1}^d ((t v_j))_{n_{ij}} \right) \mathrm{d}\boldsymbol{v} \\ & \sim t^{\alpha_0-1} \prod_{i=1}^d \left( t^{-1} \prod_{j=1}^k t^{m_{ij}} \right) \sim t^{\alpha_0-1} \prod_{i=1}^d t^{m_{i \bullet}-1} \sim t^{\alpha_0 + m-d-1}, \end{align}\] and since we are assuming \(n_i > 0\), then \(m_{i \bullet} \ge 1\) and therefore \(m \ge d\). See also the proof of Proposition 17.

9.15 Proof of Proposition 16↩︎

Recall that the ascending factorial \(((s))_q\), for integer \(q\), can be written as \[\label{eq:ascending95factorial} ((s))_{q} = \sum_{h=0}^{q} S(q,h)\,s^{h},\tag{21}\] where \(S(q,h)\) are the unsigned Stirling numbers of the first kind, defined through the recursive relation \(S(q+1,h) = q\,S(q,h) + S(q,h-1)\), with initial conditions \(S(0,0) = 1\) and \(S(q,0) = S(0,h) = 0\) for \(q > 0\) or \(h > 0\). Since \(S(q,0) = 0\) whenever \(q>0\), the summation above can start from \(1\) if \(q\) is strictly positive.

From Proposition 14, the density of \(\alpha J_{0j}\), for each \(j = 1, \dots, k\), can be rewritten as \[\begin{align} f_{\scriptscriptstyle \alpha J_{0j} \mid \boldsymbol{U}}(t) & \propto t^{-1} \, e^{-\lambda(\boldsymbol{U})\,t} \prod_{i=1}^d \left( \sum_{h_{ij}=m_{ij}}^{n_{ij}} S(n_{ij},h_{ij})\,t^{h_{ij}} \right) \\ & \propto t^{-1} \, e^{-\lambda(\boldsymbol{U})\,t} \sum_{h_j=m_{\bullet j}}^{n_{\bullet j}} \sum_{\substack{h_{1j}+\dots+h_{dj} = h_j \\ m_{ij} \le h_{ij} \le n_{ij}}} \prod_{i=1}^d S(n_{ij},h_{ij}) \, t^{h_j} \\ & \propto \sum_{h_j=m_{\bullet j}}^{n_{\bullet j}} \Bigg( \sum_{\substack{h_{1j}+\dots+h_{dj} = h_j \\ m_{ij} \le h_{ij} \le n_{ij}}} \prod_{i=1}^d S(n_{ij},h_{ij}) \Bigg) \, t^{h_j-1} e^{-\lambda(\boldsymbol{U})\,t} \end{align}\] where \(m_{ij} \in \{0,1\}\) is the indicator for \(n_{ij} > 0\), that is, \(m_{ij} = \min(1,n_{ij})\), and \(m_{\bullet j} = \sum_{i=1}^d m_{ij}\). Note that since \(n_{\bullet j} > 0\), then \(m_{\bullet j} > 0\), and the density of \(\alpha J_{0j}\) is properly defined. Moreover, for each \(j = 1, \dots, k\), define the coefficients \(S(n_{1j},\dots,n_{dj}; h_j)\) where \[\label{eq:multivariate95stirling95supp} S(q_{1},\dots,q_{d}; h) = \sum_{\substack{h_{1}+\dots+h_{d} = h \\ 0 \le h_{i} \le q_{i}}} \prod_{i=1}^d S(q_{i},h_{i}).\tag{22}\] Exploiting the recursive relation for unsigned Stirling numbers of the first kind, we obtain \[\begin{align} S(q_{1},\dots, & q_{\ell}+1, \dots,q_{d}; h) = \\ & = \sum_{\substack{h_{1}+\dots+h_{d} = h \\ 0 \le h_{i} \le q_{i} \\ 0 \le h_\ell \le q_\ell+1}} S(q_\ell+1,h_\ell) \prod_{\substack{i=1 \\ i \neq \ell}}^d S(q_{i},h_{i}) \\ & = \sum_{\substack{h_{1}+\dots+h_{d} = h \\ 0 \le h_{i} \le q_{i} \\ 0 \le h_\ell \le q_\ell+1}} q_\ell S(q_\ell,h_\ell) \prod_{\substack{i=1 \\ i \neq \ell}}^d S(q_{i},h_{i}) + \sum_{\substack{h_{1}+\dots+h_{d} = h \\ 0 \le h_{i} \le q_{i} \\ 0 \le h_\ell \le q_\ell+1}} S(q_\ell,h_\ell-1) \prod_{\substack{i=1 \\ i \neq \ell}}^d S(q_{i},h_{i}) \\ & = q_\ell \sum_{\substack{h_{1}+\dots+h_{d} = h \\ 0 \le h_{i} \le q_{i} \\ 0 \le h_\ell \le q_\ell}} S(q_\ell,h_\ell) \prod_{\substack{i=1 \\ i \neq \ell}}^d S(q_{i},h_{i}) + \sum_{\substack{h_{1}+\dots+h_{d} = h-1 \\ 0 \le h_{i} \le q_{i} \\ 0 \le h_\ell \le q_\ell}} S(q_\ell,h_\ell) \prod_{\substack{i=1 \\ i \neq \ell}}^d S(q_{i},h_{i}) \\ & = q_\ell \,S(q_{1},\dots,q_{d}; h) + S(q_{1},\dots,q_{d}; h-1). \end{align}\] Notably, in the third line we have used the fact that \(S(q_\ell,q_\ell+1) = 0\) (first term) and applied the change of index \(h_\ell \mapsto h_\ell+1\) (second term). Therefore, the coefficients above satisfy the recurrence relation in 10 defining the multivariate Stirling numbers.

9.16 Proof of Proposition 17↩︎

The density of \(\alpha T\) in Proposition 15 can be rewritten using the equality in 21 and the definition of multivariate Stirling numbers in 22 : \[\begin{align} f_{\scriptscriptstyle \alpha T}&(t) \propto \\ & \propto t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} \, \prod_{i=1}^d \frac{1}{((t))_{n_i}} \int_{\Delta^k}v_0^{\alpha_0-1} \prod_{j=1}^k \left( v_j^{-1}\prod_{i=1}^d \left( \sum_{h_{ij}=m_{ij}}^{n_{ij}} S(n_{ij},h_{ij})\,v_j^{h_{ij}} t^{h_{ij}} \right) \right) \mathrm{d}\boldsymbol{v} \\ & \propto t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} \,\prod_{i=1}^d \frac{1}{((t))_{n_i}} \int_{\Delta^k}v_0^{\alpha_0-1} \prod_{j=1}^k \left( \sum_{h_j=m_{\bullet j}}^{n_{\bullet j}} S(\boldsymbol{n}_j; h_j) \,v_j^{h_j-1} t^{h_j} \right) \mathrm{d}\boldsymbol{v} \\ & \propto t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} \, \prod_{i=1}^d \frac{1}{((t))_{n_i}} \int_{\Delta^k}v_0^{\alpha_0-1} \sum_{h=m}^{n} \Bigg(\sum_{\substack{h_{1}+\dots+h_{k} = h \\ m_{\bullet j} \le h_{j} \le n_{\bullet j}}} \prod_{j=1}^k S(\boldsymbol{n}_j; h_j) \, v_j^{h_j-1} t^{h_j} \Bigg)\, \mathrm{d}\boldsymbol{v}, \end{align}\] where \(\boldsymbol{n}_j = (n_{1j},\dots,n_{dj})\). The density of \(\alpha T\) is then obtained computing the integral over the \(k\)-dimensional simplex: \[\begin{align} f_{\scriptscriptstyle \alpha T}&(t) \propto \\ & \propto t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} \prod_{i=1}^d \frac{1}{((t))_{n_i}} \Bigg(\sum_{h=m}^{n} t^h\sum_{\substack{h_{1}+\dots+h_{k} = h \\ m_{\bullet j} \le h_{j} \le n_{\bullet j}}} \prod_{j=1}^k S(\boldsymbol{n}_j; h_j) \int_{\Delta^k}v_0^{\alpha_0-1} \prod_{j=1}^k v_j^{h_j-1} \, \mathrm{d}\boldsymbol{v}\Bigg) \\ & \propto t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} \prod_{i=1}^d \frac{1}{((t))_{n_i}} \Bigg(\sum_{h=m}^{n} \frac{t^h}{((\alpha_0))_h} \sum_{\substack{h_{1}+\dots+h_{k} = h \\ m_{\bullet j} \le h_{j} \le n_{\bullet j}}} \prod_{j=1}^k \Gamma(h_j)\,S(\boldsymbol{n}_j; h_j)\Bigg) \\ & \propto t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} \prod_{i=1}^d \frac{1}{((t))_{n_i}} \left(\sum_{h=m}^{n} \frac{c_h }{((\alpha_0))_h}\, t^{h} \right), \end{align}\] where the coefficients \(c_h\) for \(h = m,\dots,n\) are defined as \[c_h = \sum_{\substack{h_{1}+\dots+h_{k} = h \\ m_{\bullet j} \le h_{j} \le n_{\bullet j}}} \prod_{j=1}^k \Gamma(h_j)\,S(n_{1j},\dots,n_{dj}; h_j).\]

10 Posterior sampling for the gamma-gamma hCRV↩︎

This section provides additional details on the practical implementation of posterior sampling algorithms for the normalized gamma-gamma hCRV described in Section 5, and contains numerical illustrations supporting their effectiveness. Section 10.1 analyzes the Metropolis-Hastings steps of Algorithm 2. Specifically, a random walk Metropolis-Hastings scheme on the log-scale with Gaussian increments is proposed and compared with the approach of [28], based on gamma proposals. The initialization of Algorithm 3 is detailed in Sections 10.2 and 10.3, which respectively discuss the computation of coefficients \(c_h\)’s in Proposition 17 and the optimization of parameter \(r\) within the rejection sampling scheme. Section 10.4 describes the inverse Lévy measure algorithm to obtain a truncated sample from the hierarchy of gamma CRMs in Proposition 14(a). For this purpose, an efficient procedure to sequentially invert the exponential integral function is outlined in Section 10.5. This same simulation strategy remains valid a priori to sample from gamma-gamma hCRVs. Section 10.6 characterizes the posterior random probabilities \(\tilde{p}_i = \tilde{\mu}_i / \tilde{\mu}_i(\mathbb{X}) \mid \boldsymbol{X}_{1:d}\) as conditionally normalized gamma CRMs, i.e. conditionally Dirichlet processes, and thus provides an alternative approach to directly sample their posterior normalized jumps. This allows for a straightforward comparison with standard samplers for the HDP with gamma prior on the concentration parameter, which target the same posterior distributions (Proposition 6). The marginal Gibbs sampler of [1], based on the restaurant franchise metaphor, is tailored to our setting in Section 10.7. This section also proposes an alternative collapsed Gibbs sampler for the HDP, which directly samples the number of tables serving each dish, and may be of independent interest. Numerical illustrations of the effectiveness of the proposed algorithms in terms of mixing properties and posterior accuracy are provided in Section 10.8. Finally, Section 10.9 extends the simulation study of Section 5.4 comparing the different algorithms in terms of execution time per effective sample, as the number \(k\) of distinct values increases.

10.1 Metropolis-Hastings steps in Algorithm 2↩︎

The non-standard steps in the posterior sampling algorithms of Section 5.3 are:

• the marginal sampling of random variable \(\alpha T\) in ?? , whose joint density with the auxiliary vector \(\boldsymbol{V}\) is known up to a normalizing constant;

• the sampling of random variables \(\alpha J_{01}, \dots, \alpha J_{0k}\) in ?? , whose densities are again known up to normalizing constants.

A natural approach to obtain samples from non-standard densities is resorting to MCMC schemes. In Algorithm 2, we consider a blocked Gibbs sampler with Metropolis-Hastings steps. The full conditional distributions are derived from Propositions 14 and 15, namely \[\begin{align} (V_j,V_\ell) \mid \boldsymbol{V}^{-j,\ell}, (\alpha T) = t & \, \sim \, v_j^{-1} \,v_\ell^{-1} \prod_{i=1}^d ((t v_j))_{n_{ij}} ((t v_\ell))_{n_{i\ell}} \qquad (j, \ell = 1, \dots, k), \\ (V_j,V_0) \mid \boldsymbol{V}^{-j,0}, (\alpha T) = t & \, \sim \, v_0^{\alpha_0-1} \,v_j^{-1}\prod_{i=1}^d ((t v_j))_{n_{ij}} \qquad (j = 1, \dots, k), \\ (\alpha T) \mid \boldsymbol{V} & \, \sim \, t^{\alpha_0-1} \, e^{-(b_0/\alpha)\,t} \prod_{i=1}^d \frac{1}{((t))_{n_i}} \prod_{j=1}^k \prod_{i=1}^d ((t v_j))_{n_{ij}}, \\ (\alpha J_{0j}) \mid \boldsymbol{U} & \, \sim \, t^{-1} \, e^{-\lambda(\boldsymbol{U})\,t} \,\prod_{i=1}^d ((t))_{n_{ij}} \qquad (j = 1, \dots, k), \end{align}\] where \(\boldsymbol{V}^{-j,\ell}\) denotes the vector \(\boldsymbol{V} \in \Delta^k\) with components \(V_j\) and \(V_\ell\) removed. We note that variables in \(\boldsymbol{V} \in \Delta^k\) are subject to the constraint \(\sum_{j=0}^k V_j = 1\), and one cannot update variables independently; on the other hand, variables \(\alpha T\) and \(\alpha J_{0j}\) can take every positive value. Throughout the Gibbs sampling procedure, we propose new values for each variable, and accept or reject the proposal according to the Metropolis-Hastings ratio. The rest of the section outlines some of the possible proposals, focusing on rather straightforward and widely known options. Clearly, more sophisticated alternatives can be considered without undermining the validity of our results.

A simple and effective symmetric proposal for the pair of auxiliary variables \((V_j,V_\ell)\), for each \(j,\ell = 0, \dots, k\), is \[v_{j}^* = \varepsilon \, (v_j + v_\ell), \qquad v_{\ell}^* = (1-\varepsilon) \, (v_j + v_\ell),\] where \(v_j\) and \(v_\ell\) are the current values of \(V_j\) and \(V_\ell\) and \(\varepsilon \sim \mathcal{U}(0,1)\); see e.g. [80]. The proposal is accepted with log-probability \[\log(r) = \min \left\{0, \;q_j(v^*_j) - q_j(v_j) + q_\ell(v^*_\ell) - q_\ell(v_\ell) \right\},\] where \(q_0(v) = (\alpha_0 - 1) \, v\), and \(q_j(v) = \sum_{i=1}^d \log((t v))_{n_{ij}} - v\), for \(j = 1, \dots, k\). At each iteration of the Gibbs sampler, we randomly select \(k\) pairs of indexes \((j,\ell)\) and implement the Metropolis-Hastings step detailed above. Note that selecting \(k\) pairs is a good compromise between choosing a single pair and scanning each of the \(k(k-1)/2\) possible pairs.

For the positive variables \(\alpha T\) and \(\alpha J_{0j}\)’s, we consider random walks with two alternative proposals, which we first describe for a generic positive random variable with density \(f(x)\).

• Gamma proposal. Following the approach discussed in [28], we propose a new value \(x^*\) from a gamma distribution centered at the current value \(x\), \[x^* \sim \text{Gamma}(\delta, \delta / x),\] where \(\delta > 0\) controls the variance of the proposal. The proposed value \(x^*\) is accepted with log-probability \[\begin{align} \log(r) & = \min \big\{0, \, \log f(x^*) - \log f(x) + \log g(x; \delta, \delta/x^*) - \log g(x^*; \delta, \delta/x) \big\} \\ & = \min \left\{ 0, \, \log f(x^*) - \log f(x) + (2\delta-1)(\log x - \log x^*) + \delta \left(\frac{x^*}{x} - \frac{x}{x^*} \right) \right\}, \end{align}\] where \(x \mapsto \log g(x; a, b)\) is the log-density of a gamma distribution with shape \(a > 0\) and rate \(b > 0\). For the practical implementation, [28] suggest to restrict to \(\delta \ge 1\).

• Random walk on log-transform. For a positive random variable, we can target the density of its log-transform, which is \(t \mapsto f(e^t)\,e^t\). In this case, we resort to a random walk on the log-scale, and propose a new value \(\log x^*\) from a normal distribution centered at the current value \(\log x\), \[\log x^* \sim \mathcal{N}(\log x, \sigma^2),\] where \(\sigma^2 > 0\) controls the variance of the proposal. The proposed value \(x^*\) is accepted with log-probability \[\log(r) = \min \big\{0, \, \log f(x^*) - \log f(x) + \log x^* - \log x \big\}.\]

Interestingly, both proposals can be interpreted as multiplicative perturbations of the current value, with different distributions for the perturbation factor. Indeed, they both propose a new value \(x^* = x\,\varepsilon\), where \(x\) is the current value and \(\varepsilon\) is distributed as either (a) a gamma random variable \(\varepsilon \sim \text{Gamma}(\delta, \delta)\) with variance \(1/\delta\), or (b) a log-normal random variable \(\log \varepsilon \sim \mathcal{N}(0, \sigma^2)\) with variance \(\sigma^2\) on the log-scale. Therefore, the acceptance rates can be rewritten as \[\begin{align} & {\rm (a)} \qquad \log(r) = \min \left\{ 0, \, \log g(x\,\varepsilon) - \log g(x) + \delta (\varepsilon - 1/\varepsilon) - 2\delta\, \log \varepsilon \right\}, \\ & {\rm (b)} \qquad \log(r) = \min \left\{ 0, \, \log g(x\,\varepsilon) - \log g(x) \right\}, \end{align}\] where \(\log g(x) = \log f(x) + \log x\). The practical implementation of these approaches requires the evaluation of the logarithm of target density at the current and proposed values, up to normalizing constants. For our purposes, we have \[\begin{gather} \log f_{ \scriptscriptstyle \alpha T \mid \boldsymbol{V}}(t) = (\alpha_0-1) \log t - (b_0/\alpha)\,t + \sum_{j=1}^k \sum_{i=1}^d \log((tv_j))_{n_{ij}} - \sum_{i=1}^d \log((t))_{n_i}, \\ \log f_{\scriptscriptstyle \alpha J_{0j} \mid \boldsymbol{U}}(t) = \sum_{i=1}^d \log((t))_{n_{ij}} - \lambda(\boldsymbol{U}) \, t - \log t \qquad (j = 1, \dots, k). \end{gather}\] Note that the random variables \(\alpha J_{0j}\)’s are gamma distributed whenever \(n_{1j}, \dots, n_{dj} \le 1\), and can thus be sampled exactly as \[\label{eq:exact95jumps} \alpha J_{0j} \mid \boldsymbol{U} \sim \text{Gamma}\left( \sum_{i=1}^d n_{ij}, \, \lambda(\boldsymbol{U}) \right).\tag{23}\] Integrating this feature into the algorithm leads to a substantial reduction in computational times in many scenarios. A crucial aspect of random walk Metropolis-Hastings algorithms is the choice of a suitable variance parameter for the proposal, in order for the resulting Markov chain to mix efficiently [81]. For this purpose, we resort to adaptive MCMC techniques, which allow to automatically tune the proposals to optimize their performances. Specifically, we implement a simple but effective updating rule, based on the Robbins-Monro recursion [82]. The explicit updates for the two proposals discussed above are \[{\rm (a)} \quad \log \delta_{s+1} \leftarrow \log \delta_s - \gamma_s (r_s - r^*), \qquad {\rm (b)} \quad \log \sigma^2_{s+1} \leftarrow \log \sigma^2_s + \gamma_s (r_s - r^*),\] where \(\delta_s\) and \(\sigma^2_s\) are the parameters of the proposals at step \(s\), \(r_s\) is the probability of acceptance at step \(s\), \(r^* = 0.44\) is the target acceptance rate [83] and \((\gamma_s)_{s \ge 1}\) is a sequence of positive and non-increasing numbers. In our implementation, we employ the deterministic sequence \(\gamma_s = (10 + s)^{-1/2}\), for \(s \ge 1\); see [82] for a discussion on the desirable properties of such sequence. This adaptive procedure is employed at each iteration during the burn-in phase, and the values of variance parameters reached at the end of the burn-in are kept fixed for the actual posterior sampling. Our implementation allows for different variance parameters for \(\alpha T\) and each of the \(\alpha J_{0j}\)’s.

We conclude this section with a comparison of the gamma and log-normal proposals, carried out through a simple numerical experiment. We consider \(1,000\) synthetic datasets, each consisting of a random number of groups \(d \sim \text{Poisson}(5)\), with \(n_i = 100\) observations per group. Observations are sampled from a hierarchical Dirichlet process with random concentration parameters \(\alpha \sim \text{Gamma}(5)\) and \(\alpha_0 \sim \text{Gamma}(3)\). We draw \(10,000\) posterior samples for each dataset, with a burn-in of \(1,000\) steps, during which variance parameters are tuned. Table 1 reports the effective sample size (ESS) and acceptance rate for each target variable and proposal. Figures for the \(\alpha J_{0j}\)’s are the means computed over \(j\), excluding those indexes for which the exact sampling in 23 is available. The log-normal proposal outperforms the gamma proposal in terms of ESS for both \(\alpha T\) and the \(\alpha J_{0j}\)’s; the difference is particularly evident for the latter. Moreover, the adaptive MCMC procedure systematically reaches the target acceptance rate of \(0.44\). Clearly, the evidence provided by this experiment is not conclusive in any respect, as we are restricting ourselves to a particular algorithm and sampling scenario. However, the synthetic datasets cover a wide range of possible count matrices \((n_{ij})_{ij}\) that one may encounter in applications. Therefore, we recommend using the log-normal proposal, which we also consider for the simulations in Section 5.4 of the main manuscript.

10.2 Computing the coefficients in Proposition 17↩︎

In this section, we show that the coefficients \(c_h\), for \(h = m, \dots, n\), defined in Proposition 17, can be computed through a sequence of discrete convolutions. For each \(\ell = 1, \dots, k\), define the vector \(\boldsymbol{c}(\ell) = (c(\ell,h))_h\) with \[c(\ell; h) = \sum_{\substack{h_1 + \cdots + h_\ell = h \\ m_{\bullet j} \le h_j \le n_{\bullet j}}} \prod_{j=1}^\ell a(j; h_j), \qquad h = \sum_{j=1}^\ell m_{\bullet j}, \dots, \sum_{j=1}^\ell n_{\bullet j},\] where \(a(j; h) = \Gamma(h) \, S(n_{1j},\dots,n_{dj}; h)\), for \(h = m_{\bullet j}, \dots, n_{\bullet j}\) and \(j = 1, \dots, \ell\). Thus, the \(c_h\)’s are such that \(c_h = c(k; h)\), that is, coincide with \(\boldsymbol{c}(k)\). The entries of vector \(\boldsymbol{c}(\ell)\) can be computed from vector \(\boldsymbol{c}(\ell-1)\) as follows: \[\begin{align} c(\ell; h) & = \sum_{\substack{h_1 + \cdots + h_\ell = h \\ m_{\bullet j} \le h_j \le n_{\bullet j}}} \prod_{j=1}^{\ell} a(j; h_j) = \sum_{h_{\ell} = m_{\bullet \ell}}^{n_{\bullet \ell}} \left( \sum_{\substack{h_1 + \cdots + h_{\ell-1} = h - h_{\ell} \\ m_{\bullet j} \le h_j \le n_{\bullet j}}} \prod_{j=1}^{\ell-1} a(j; h_j) \right) a(\ell; h_\ell) \\ & = \sum_{h_{\ell} = m_{\bullet \ell}}^{n_{\bullet \ell}} c(\ell-1; h-h_\ell) \, a(\ell; h_\ell). \end{align}\] In other words, the vector \(\boldsymbol{c}(\ell)\) is obtained by convolution between vector \(\boldsymbol{c}(\ell-1)\) and vector \(\boldsymbol{a}(\ell) = (a(\ell;h),h)\). The computational cost for this operation is \(\mathcal{O}\left( n_{\bullet \ell} \sum_{j=1}^{\ell-1} n_{\bullet j} \right)\). Therefore, the coefficients \((c_h)_h = \boldsymbol{c}(k)\) can be computed through the recursive relation \[\boldsymbol{c}(0) = (1), \qquad \boldsymbol{c}(\ell) = \boldsymbol{c}(\ell-1) * \boldsymbol{a}(\ell) \quad (\ell = 1, \dots, k),\] or, equivalently, \(\boldsymbol{c}(k) = \boldsymbol{a}(1) * \cdots * \boldsymbol{a}(k)\), where \(*\) is the convolution operator. The total computational cost is \(\mathcal{O}\left( \sum_{j<\ell} n_{\bullet j} n_{\bullet \ell} \right)\).

10.3 Optimal choice of parameter \(r\)↩︎

Section 5.2 outlines a rejection sampling algorithm for sampling the latent variable \(\alpha T\) from its density ?? . Specifically, values are proposed from \({\rm Gamma}(\alpha_0 + r, b_0 / \alpha)\) and accepted with probability proportional to \(t^{-r} R(t)\), where \(r\) is a real parameter and \[R(t) = \prod_{i=1}^d \frac{1}{((t))_{n_i}} \left( \sum_{h=m}^n \frac{c_h}{((\alpha_0))_h} \, t^h \right).\] A necessary condition for the rejection sampling scheme is \(t^{-r} R(t)\) to be bounded above for \(t \ge 0\). Since \(R(t)\) is a ratio of polynomials, both having degree \(n\) and non-negative coefficients, it is a continuous function for \(t > 0\). Moreover, when \(q\) is a positive integer, \(((s))_q \sim s\) for \(s \to 0\) while \(((s))_q \sim s^q\) for \(s \to \infty\). Hence, \(R(t) \sim t^{m-d}\) for \(t \to 0\) and \(R(t) \sim c_n / ((\alpha_0))_n\) for \(t \to \infty\), which implies that \(t^{-r} R(t)\) is continuous and bounded for \(t \ge 0\) when \(0 \le r \le m-d\). Our goal is choosing the value of \(r\) within this interval such that the acceptance probability is maximized.

For this purpose, let \(t^*(r)\) be the value of \(t\) that maximizes \(t^{-r} R(t)\) for \(t \ge 0\). The overall acceptance probability is \[\begin{align} & \mathbb{E}_T \left[ \left(\frac{t^*(r)}{t}\right)^r \frac{R(t)}{R(t^*(r))} \right] =\\ & = \frac{t^*(r)^r}{R(t^*(r))} \int_0^\infty t^{-r} R(t) \left( \frac{b_0}{\alpha} \right)^{\alpha_0+r} \Gamma(\alpha_0+r)^{-1} t^{\alpha_0+r-1} e^{-(b_0/\alpha)\,t} \mathrm{d}t \\ & = \frac{t^*(r)^r}{R(t^*(r))} \left( \frac{b_0}{\alpha} \right)^{\alpha_0+r} \Gamma(\alpha_0+r)^{-1} \int_0^\infty t^{\alpha_0-1} e^{-(b_0/\alpha)\,t} R(t) \, \mathrm{d}t. \end{align}\] Finding the value of \(0 \le r \le m-d\) that maximizes the quantity above is equivalent to finding the value maximizing its logarithm, discarding terms not depending on \(r\), that is \[r^* = \arg\max_r \left\{ r \log t^*(r) - \log R(t^*(r)) + (\alpha_0 + r) \log (b_0/\alpha) - \log \Gamma(\alpha_0 + r) \right\}.\] This maximization problem can be further simplified by restricting to a finite set of potentially maximizing values. Indeed, for \(0 < r < m-d\), the function \(R(t)\) is continuous and differentiable for \(t > 0\), and such that \(t^{-r} R(t) \to 0\) for both \(t \to 0\) and \(t \to \infty\). Hence, \(t^*(r)\) is a stationary point for \(t^{-r} R(t)\), which implies \(R'(t^*(r)) \, t^*(r) = r\,R(t^*(r))\). Moreover, by the implicit function theorem, \(t^*(r)\) is a continuous and differentiable function in \(r\). Therefore, the objective function is continuous and differentiable for \(0 < r < m-d\), and the set of potentially maximizing points in \((0,m-d)\) may be restricted to the stationary points (if any), satisfying \[\log t^*(r) + \log (b_0/\alpha) - \Psi(\alpha_0 + r) = 0,\] where \(\Psi\) denotes here the digamma function \(\Psi(x) = \frac{\mathrm{d}}{\mathrm{d}x} \log \Gamma(x)\), and we have used that \(t^*(r)\) is a stationary point for \(t^{-r} R(t)\) and thus \(R'(t^*(r)) \, t^*(r) = r\,R(t^*(r))\). Remarkably, this stationarity condition is also satisfied at the boundary of the maximization set. Indeed, for \(r = 0\), the optimal \(t^*(r)\) may be either a stationary point or \(+\infty\), for which \(R'(t^*(r)) = 0\) necessarily holds. On the other hand, if \(r = m-d\), the optimal \(t^*(r)\) may be either a stationary point or \(0\), for which \(R(t^*(r)) = 0\) holds.

10.4 Sampling from the hierarchy of gamma CRMs↩︎

From Proposition 14(a), the residual component \(\boldsymbol{\tilde{\mu}}^*\) of the posterior distribution of \(\boldsymbol{\tilde{\mu}}\) retains a hierarchical structure, conditionally on latent variables \(\boldsymbol{U}\): \[\begin{align} \tilde{\mu}_1^*,\dots,\tilde{\mu}_d^* \mid \tilde{\mu}_0^*, \boldsymbol{U} & \sim \prod_{i=1}^d \mathrm{CRM}\left( \alpha \, s^{-1}\, e^{-b(1 + U_i/b)\,s} \mathrm{d}s \otimes \tilde{\mu}_0^*\right), \\ \tilde{\mu}_0^* \mid \boldsymbol{U} & \sim \mathrm{CRM}\left( \alpha_0 \, s^{-1}\, e^{-\alpha\, \lambda(\boldsymbol{U})\,s} \mathrm{d}s \otimes P_0 \right). \end{align}\] Therefore, one needs to sample from a hierarchy of gamma CRMs in order to obtain complete samples from the posterior.

At the root of the hierarchy, approximate posterior samples from the rescaled gamma random measure \(\alpha \tilde{\mu}_0^* \mid \boldsymbol{U}\) can be obtained through the Ferguson-Klass representation [43]. This amounts to sampling the largest \(L\) jumps of an infinitely active random measure in decreasing order, and thus provides its best finite-dimensional approximation: for a fixed truncation level, the approximation error is minimized. A straightforward approach uses the inverse Lévy measure algorithm [44], [45]. For \(\ell = 1, \dots, L\), let \(\omega_{0\ell} \ge 0\) be the value solving the equation \[\frac{\xi_\ell}{\alpha_0} = \int_{\omega_{0\ell}}^\infty s^{-1}\, e^{-\alpha\, \lambda(\boldsymbol{U})\,s} \mathrm{d}s = E_1 ( \alpha\, \lambda(\boldsymbol{U}) \,\omega_{0\ell} ) \qquad (\ell = 1, \dots, L),\] where \(E_1\) is the exponential integral function and \(\xi_1 < \dots < \xi_L\) a.s. are the first \(L\) jump times of a unit rate Poisson process, that is, \(\xi_0 = 0\) and the inter-arrival times are \(\xi_{\ell} - \xi_{\ell-1} \stackrel{\scriptsize{\rm iid}}{\sim}{\rm Exp}(1)\), for \(\ell=1,\dots,{L}\). Then the random measure \(\alpha \tilde{\mu}_0^* \mid \boldsymbol{U}\) can be approximated as \[\alpha \tilde{\mu}_0^* \mid \boldsymbol{U} \approx \sum_{\ell=1}^L (\alpha \omega_{0\ell}) \, \delta_{\phi_\ell},\] where \(\phi_1,\dots,\phi_L \stackrel{\scriptsize{\rm iid}}{\sim}P_{0}\) are independent from the Poisson process \((\xi_1, \dots, \xi_L)\). This algorithm requires to sequentially invert the exponential integral function numerically, which is a nontrivial albeit much investigated task. Details about our implementation are provided in the following Section 10.5. Alternative approaches for sampling the largest \(L\) jumps of a random measure in decreasing order, and their specifications for the gamma process, are explored in [47], [48]; see also [46].

The random measures in \(\boldsymbol{\tilde{\mu}}^*\) are conditionally independent given \(\alpha \tilde{\mu}_0^*\) and \(\boldsymbol{U}\), and each component \(\tilde{\mu}_i^* \mid \tilde{\alpha} \mu_0^*,\boldsymbol{U}\) can be approximately sampled from a gamma CRM having Lévy intensity \[\mathrm{d}\nu^*_i(s,x) = \sum_{\ell=1}^L (\alpha \omega_{0\ell}) \, s^{-1} e^{-b(1+U_i/b)\,s}\,\mathrm{d}s \, \mathrm{d}\delta_{\phi_\ell}(x).\] The additive components of the Lévy measure represent independent summands for \(\tilde{\mu}_i^*\) concentrated at different fixed location. Hence, \(\tilde{\mu}_i^* \mid \alpha \tilde{\mu}_0^*,\boldsymbol{U}\) can be approximated as \[\tilde{\mu}_i^* \mid \alpha \tilde{\mu}_0^*,\boldsymbol{U} \approx \sum_{\ell=1}^L \omega_{i\ell}\,\delta_{\phi_\ell},\] where \(\omega_{i1}, \dots \omega_{iL}\) are independent random variables with \(\omega_{i\ell} \sim \text{Gamma}(\alpha \omega_{0\ell}, b(1+U_i/b))\), for \(\ell = 1, \dots, L\). The sampling procedure to obtain an approximation of \(\boldsymbol{\tilde{\mu}}\) by truncation of its infinite sequence of jumps is summarized in Figure 7; the similarities with the sampling algorithms for jumps \(\boldsymbol{J}\), depicted in Figure 1, are evident. Note that the total mass of each random measure \(\tilde{\mu}^*_i\), that is, the mass of the posterior random measure \(\tilde{\mu}_i \mid \boldsymbol{X}_{1:d}\) not assigned to fixed locations, can instead be sampled exactly from a hierarchy of gamma random variables.

This same sampling procedure can be also employed to sample a priori from the gamma-gamma hCRV, with minimal adjustments. Indeed, the prior construction introduced in Example 1 coincides with the hierarchy in Proposition 14(a) displayed above when \(U_i = 0\) for each \(i = 1,\dots,d\), and thus \(\lambda(\boldsymbol{U}) = b_0 / \alpha\). Specifically, approximate samples from the rescaled gamma random measure \(\alpha \tilde{\mu}_0\) can be obtained via the inverse Lévy measure algorithm, while increments of each random measure \(\tilde{\mu}_i \mid \alpha \tilde{\mu}_0\) can be sampled independently from gamma distributions, as discussed in Remark 3.

10.5 Inverting the exponential integral↩︎

The implementation of the inverse Lévy measure algorithm of [45] for the gamma CRM requires to invert the tail integrals of its Lévy density (Section 10.4). This amounts to sequentially invert the exponential integral function \(E_1\), defined for \(x > 0\), \[E_1(x) = \int_x^{+\infty} s^{-1} e^{-s} \mathrm{d}s.\] Note that \(E_1\) is a strictly decreasing function, and thus is invertible. A convenient approach to find its inverse \(E_1^{-1}(y)\), for a given value \(y > 0\), is to determine the unique root of the function \(x \mapsto E_1(x) - y\), exploiting root-finding algorithms such as Newton’s method. This method requires the evaluation of the derivative \(E_1'\) of the exponential integral, which can be computed in closed form, as it equals the opposite of the integrand function, \(E_1'(x) = - x^{-1}e^{-x}\).

A well-known limitation of Newton’s method is the possibility to obtain iteration values that fall outside the domain of the function, where its evaluation is not possible. This situation is particularly relevant to inverse Lévy measure algorithms, as the tail integrals of infinitely divisible Lévy densities diverge to \(+\infty\) as the lower bound of the integration interval goes to \(0\). A simple workaround that typically solves this issue is choosing a starting point \(x_0\) for Newton’s algorithm which falls on the left of the solution, that is \(E_1(x_0) \ge y\). This can be achieved by iteratively halving an initial guess. Remarkably, within the sequential approach required by the algorithm in Section 10.4, the standard choice for the initial guess is the solution at the previous step. However, this falls on the right of the current solution, and thus halving is always necessary. Convergence of this algorithm is guaranteed whenever the Lévy density is a decreasing function, and thus its tail integral is decreasing and convex. Figure 8 shows an illustration of the different behaviours of Newton’s method, depending on the starting point \(x_0\).

Instead, we consider an alternative approach that improves the efficiency of Newton’s method. Specifically, we redefine the exponential integral through the logarithm of its argument; in other words, we invert the function \[f(z) = E_1(e^z) = \int_{e^z}^{+\infty} s^{-1} e^{-s} \mathrm{d}s,\] whose derivative is given by \(f'(z) = E_1'(e^z) \,e^z = - \exp( - e^z)\). This formulation has the advantage of not being restricted to positive values, thus preventing iterations to fall outside the domain. Note that \(f\) is decreasing and convex, which guarantees convergence of Newton’s method for every starting point \(x_0\).

A decisive advantage of this approach is the asymptotic linearity of function \(f\) as \(z\) diverges to \(-\infty\). Specifically, \(f(z) \approx - \gamma - z\) for \(z \to -\infty\), where \(\gamma\) is the Euler-Mascheroni constant. This property is crucial for speeding up convergence of the numerical scheme, as the rate of convergence of Newton’s method is proportional to the second derivative around the solution. The asymptotic expansion is also useful to improve the numerical stability of function evaluations. We argue that redefining the tail integral as a function on the whole real line, removing the constraint to positive values, may be a general technique to enhance the performances of Newton’s method for computing the inverse Lévy measure, beyond the gamma process case.

10.6 Posterior distribution of normalized random measures↩︎

The posterior random probabilities arising from model 4 are the normalization of the posterior random measures characterized in Theorem 10 and specialized in Proposition 14 for the gamma-gamma hCRV. The simplest approach to obtain posterior samples is thus normalizing the posterior samples from the corresponding random measures. In this section, we show that such posterior random probabilities are distributed as conditionally Dirichlet processes with discrete base measures. Therefore, their probability weights can be alternatively sampled from a Dirichlet distribution. This allows for a straightforward comparison with alternative samplers for the HDP, such as the marginal Gibbs samplers of [1]; see the following Section 10.8.

For this purpose, denote by \(\boldsymbol{\pi}\) and \(\boldsymbol{\pi}^*\) the normalized jumps of \(\tilde{\mu}_i \mid \boldsymbol{X}_{1:d}\), for \(i = 1, \dots, d\), \[\pi_{ij} = \frac{J_{ij}}{\sum_{j=1}^k J_{ij} + \tilde{\mu}^*_i(\mathbb{X})}, \qquad \pi^*_{i\ell} = \frac{\omega_{i\ell}}{\sum_{j=1}^k J_{ij} + \tilde{\mu}^*_i(\mathbb{X})},\] where the \(J_{ij}\)’s are jumps at fixed locations and the \(\omega_{i\ell}\)’s are the jumps of \(\tilde{\mu}_i^*\), arranged in decreasing order; see Section 10.4 for details on sampling from \(\tilde{\mu}_i^*\). From Theorem 10 it follows that the random probabilities \(\boldsymbol{\tilde{P}} = \boldsymbol{\tilde{\mu}}/\boldsymbol{\tilde{\mu}}(\mathbb{X})\) are distributed, a posteriori, as \[\label{eq:normalized95weights} \tilde{P}_i = \frac{\tilde{\mu}_i}{\tilde{\mu}_i(\mathbb{X})} \bigm\vert \boldsymbol{X}_{1:d} \, \stackrel{d}{=} \, \sum_{j=1}^k \pi_{ij} \,\delta_{X^*_j} + \sum_{\ell \ge 1} \pi^*_{i\ell} \,\delta_{Y_\ell} \quad (i = 1, \dots, d).\tag{24}\] For the normalized gamma-gamma hCRV, the conditional distribution of each posterior random probability is a Dirichlet process with discrete base measure.

Remarkably, the posterior distribution of \(\boldsymbol{\tilde{P}}\) does not depend on the prior rate parameter \(b\), as already discussed right after Proposition 6. Moreover, it is conditionally independent of latent variables \(\boldsymbol{U}\), given \(\alpha \boldsymbol{J}_0\) and \(\tilde{\mu}_0^*\).

This result is particularly relevant from the algorithmic point of view, as anticipated above. Indeed, after obtaining samples from variables \(\alpha J_{01}, \dots, \alpha J_{0k}\) and an approximation by truncation of \(\alpha \tilde{\mu}^*_0\), as described in Section 10.4, the probability weights in 24 can be sampled from a \((k+L)\)-dimensional Dirichlet distribution. In particular, for each \(i = 1, \dots, d\) and \(L \in \mathbb{N}\), \[\begin{gather} \big(\pi_{i1}, \dots, \pi_{ik}, \pi^*_{i1}, \dots, \pi^*_{iL} \big) \bigm\vert \boldsymbol{X}_{1:d}, \boldsymbol{J}_0, \tilde{\mu}^*_0 \\ \sim \text{Dirichlet}\big(n_{i1} + \alpha J_{01}, \dots, n_{ik} + \alpha J_{0k}, \alpha \omega_{01}, \dots, \alpha \omega_{0L} \big). \end{gather}\] The conditional independence from \(\boldsymbol{U}\) allows further parallelization of the posterior sampling schemes. The resulting structure of conditional dependencies within the complete sampling algorithms for the normalized posterior random measures is summarized in Figure 9.

10.7 Marginal Gibbs samplers for the HDP↩︎

This section provides the details for the derivation of the marginal posterior samplers for the HDP considered in Section 5.3, namely the CRF-based Gibbs sampler and the collapsed Gibbs sampler. The CRF-based Gibbs sampler is based on existing work [1], [49], which we recall here using our notation. Proposition 21 establishes a new result, which determines our collapsed Gibbs sampler; this can be seen as an extension of the sampling schemes proposed in [1] and [55].

The hierarchical Dirichlet process was first described by [1], where the restaurant franchise metaphor is introduced. According to this metaphor, customers in each restaurant eat dishes selected from a countably infinite menu that is shared across restaurants. Moreover, in each restaurant, customers eating the same dish may or may not seat at the same table, while customers eating different dishes are necessarily seated at different tables. In other words, observations are organized within a nested partition structure. Remarkably, this structure is not specific to the hierarchical Dirichlet process, but characterizes every hierarchical structure of normalized random measures considered in [49]; see also [15]. Such partition is typically described through two collections of variables, displaying ties among each collection: for each restaurant \(i = 1,\dots,d\), variables \(X_{i\ell}\) and \(Z_{i\ell}\) denote, respectively, the dish and table of the \(\ell\)-th customer. The distinct values taken by variables \(\boldsymbol{X}_{1:d} = ({\boldsymbol{X}}_1, \dots, {\boldsymbol{X}}_d)\) are denoted by \(X^*_1,\dots,X^*_k\), and represent the distinct dishes from the common menu. Similarly, for each restaurant \(i = 1,\dots,d\), the distinct values taken by the variables \(\boldsymbol{Z}_{i} = (Z_{i1}, \dots, Z_{in_i})\), are denoted by \(Z^*_{ij1},\dots,Z^*_{ijh_{ij}}\) for \(j = 1, \dots, k\), and represent the different tables serving each dish in the \(i\)-th restaurant. The particular choice of hierarchical random measures determines the joint probability distribution of \(\boldsymbol{X}_{1:d}\) and \(\boldsymbol{Z}_{1:d} = ({\boldsymbol{Z}}_1, \dots, {\boldsymbol{Z}}_d)\), namely the dishes and tables allocation, or equivalently the sequence of their predictive distributions. For a hierarchical Dirichlet process with concentration parameters \(\alpha, \alpha_0 > 0\), the joint probability of tables and dishes takes the form \[\label{eq:peppf} \mathbb{P}(\boldsymbol{X}_{1:d}, \boldsymbol{Z}_{1:d} ) = \frac{\alpha_0^k}{((\alpha_0))_h} \prod_{j=1}^k \Gamma(h_{\bullet j}) \, \prod_{i=1}^d \left( \frac{\alpha^{h_{i\bullet}}}{((\alpha))_{n_i}} \prod_{j=1}^k \prod_{r=1}^{h_{ij}} \Gamma(q_{ijr}) \right).\tag{25}\] In the expression above, \(h_{ij}\) denotes the (random) number of tables in which the \(n_{ij}\) customers eating dish \(X^*_j\) in restaurant \(i\) are partitioned, while \(q_{ijr}\), for \(r = 1,\dots,h_{ij}\), denotes the (random) number of customers seated at each of those tables.

Contrary to the approach of [1], which considers kernel mixtures, our model 4 assumes that the dishes \(\boldsymbol{X}_{1:d}\) eaten by each customer are actually observed, and thus fixed. Therefore, one only has to sample the tables \(\boldsymbol{Z}_{1:d}\), in order to obtain a complete description of the nested partition. Remarkably, the full conditional distribution of each \(Z_{i\ell}\), given \(\boldsymbol{X}_{1:d}\) and the other \(\boldsymbol{Z}_{1:d}\), can be easily derived from the joint probability 25 . According to the restaurant franchise metaphor, having observed that the \(\ell\)-th customer from restaurant \(i\) eats dish \(X^*_j\), they may either sit at tables \(Z^*_{ij1},\dots,Z^*_{ijh_{ij}}\) or open a new table. More precisely, given \(X_{i\ell} = X^*_j\), then \[\begin{align} \mathbb{P}(Z_{i\ell} = Z^*_{ijr} \mid \cdots) & \propto \frac{q_{ijr}^{-\ell}}{\alpha + n_i}, \qquad r = 1, \dots, h_{ij}, \\ \mathbb{P}(Z_{i\ell} = {\rm 'new'} \mid \cdots) & \propto \frac{\alpha}{\alpha + n_i} \, \frac{h_{\bullet j}^{-\ell}}{\alpha_0 + h}, \end{align}\] where \(q_{ijr}^{-\ell}\) and \(h_{\bullet j}^{-\ell}\) denote, respectively, the number of customers at table \(Z^*_{ijr}\) and the number of tables serving dish \(X^*_j\) across all restaurants, once customer \(\ell\) is removed from the partition structure. As a result, a straightforward Gibbs sampler consists in sequentially resampling each \(Z_{i\ell}\), for \(i = 1, \dots, d\) and \(\ell = 1, \dots, n_i\), from its conditional distribution. In Section 5.3, we refer to this sampling scheme as CRF-based Gibbs sampler. Note that, since one sequentially updates \(n\) latent variables, the computational cost is at least linear in the number of observations.

The posterior distribution of the hierarchical Dirichlet process 3 , given the allocations of the observations within the nested partition, preserves the hierarchical structure [16], \[\tilde{P}_i \mid {\boldsymbol{X}}_{1:d}, {\boldsymbol{Z}}_{1:d}, \tilde{P}_0 \sim {\rm DP}\left( \alpha \tilde{P}_0 + \sum_{j=1}^k n_{ij} \delta_{X^*_j} \right); \qquad \tilde{P}_0 \sim {\rm DP}\left( \alpha_0 P_0 + \sum_{j=1}^k h_{\bullet j} \delta_{X^*_j} \right).\] Interestingly, these distributions depend only on the numbers of customers \(n_{ij}\)’s for each restaurant and dish, and on the number of tables \(h_{\bullet j}\)’s serving each dish across all restaurants. In our setting, the counts \((n_{ij})_{ij}\) are observed, and thus the posterior distribution is fully characterized by sampling the \(h_{\bullet j}\)’s from their conditional distribution.

The direct sampling of the number of tables serving each dish in each restaurant was first suggested in [1], although conditionally on the random probability \(\tilde{P}_0\). Instead, [55] adopt a marginal approach and rely on the conditional distribution of the \(h_{ij}\)’s for the hierarchical Pitman-Yor process, as derived in [49]. The further marginalization over restaurants in Proposition 21 is particularly suitable for HDP, since its posterior structure depends only on the \(h_{\bullet j}\)’s; moreover, it conveniently reduces the dimension of the sampling space. As a result, one easily devises a Gibbs sampler for the \(h_{\bullet j}\)’s based on their joint distribution. In Section 5.3, we refer to this sampling approach as collapsed Gibbs sampler.

Finally, for a full comparison with the normalized gamma-gamma hCRV, we often consider a gamma prior on the concentration parameter \(\alpha \sim {\rm Gamma}(\alpha_0, \beta)\); see Proposition 6 for the distributional equivalence result. The update of \(\alpha\) within both Gibbs samplers outlined above is performed via a Metropolis-Hastings step, targeting the full conditional distribution \[\alpha \mid {\boldsymbol{X}}_{1:d}, {\boldsymbol{Z}}_{1:d} \sim f(s) \propto s^{\alpha_0 + h - 1}\,e^{-\beta\,s} \prod_{i=1}^d \frac{1}{((s))_{n_i}}.\] An alternative sampling scheme based on an augmentation with auxiliary beta and Bernoulli variables is detailed in the Appendix of [1]. The Julia implementation of the CRF-based and the collapsed Gibbs samplers is available at github.com/claudio
delsole/hCRV.jl, where we also provide an interface to allow their integration within the R environment.

10.8 Numerical illustrations↩︎

This section contains numerical illustrations of the posterior sampling algorithms introduced in Section 5.3. Specifically, we consider the MCMC sampler in Algorithm 2 with symmetric random-walk Metropolis-Hastings steps on the log-scale (mcmc), detailed in Section 10.1, and the exact sampler in Algorithm 3 (exact). As a reference to evaluate their effectiveness, we take the marginal Gibbs sampler of [1] for the HDP, based on the restaurant franchise metaphor, with gamma prior on the concentration parameter; see Section 10.7 above for details. These algorithms target the same posterior distribution for the random probabilities (Proposition 6), hence their posterior estimates should coincide. We consider \(d = 4\) groups of observations, each of size \(n_i = 50\), sampled from independent Poisson distributions with means \(2, 3, 4\) and \(5\). The number of distinct values in the simulated dataset is \(k = 11\). Model parameters are fixed at \(\alpha_0 = \alpha = 1\) and \(b_0 = b = 1\). We draw \(10,000\) posterior samples for each algorithm; for the MCMC schemes, we consider a burn-in of \(1,000\) steps.

Figure 10 provides insights into the posterior sampling of random variable \(\alpha T\). The MCMC algorithm shows good mixing, with an effective sample size of \(1996\). The envelope for the rejection sampling in the exact algorithm is tight, and entails an acceptance rate of \(0.765\); the optimal value for \(r\) is \(8.84\) (see Section 10.3). Instead, Figure 11 displays the MCMC traceplot and posterior distributions for some of the \(J_{0j}\)’s. Specifically, we focus on the latent jumps at values \(X^*_5 = 4\) and \(X^*_8 = 7\), with counts \(n_{\bullet 5} = 36\) and \(n_{\bullet 8} = 7\). Again, the MCMC algorithms show good mixing, with effective sample sizes above \(1,600\) for all \(J_{0j}\)’s. Finally, the posterior distributions for the jumps \(J_{ij}\) at values \(X^*_5 = 4\) and \(X^*_8 = 7\) for the \(d = 4\) groups are displayed in Figure 12. The effective sample sizes are above \(3,300\) for all \(J_{ij}\)’s.

Furthermore, we conduct a visual comparison with the marginal Gibbs sampler for the HDP in terms of the accuracy of posterior distributions. Since this algorithm outputs posterior probability weights, we sample from the posterior random probabilities according to the procedure outlined in 10.6. Considering the same experimental setting described above, we compare the MCMC sampler, the exact sampler, and the marginal Gibbs sampler for HDP, with or without gamma prior on the concentration parameter (HDP and HDPpr, respectively). The posterior distributions for the probability weights \(\pi_{ij}\) at the same values \(X^*_5 = 4\) and \(X^*_8 = 7\) are displayed in Figure 13. The effective sample sizes consistently exceed \(4,000\). As expected, estimates of the posterior distributions for the HDP model without prior visibly differ from those obtained with the other algorithms, which instead target the same posterior distributions. This difference may be likewise observed in the right plot of Figure 5.

10.9 Comparison for increasing number of distinct values↩︎

Section 5.4 of the main manuscript compares the different algorithms for posterior sampling through simulation studies. The comparison is performed in terms of computational time per effective sample, as the number of groups \(d\) or the number of observations per group \(n_i\) increase. In this section, we conduct a further simulation study fixing both the number of groups \(d = 10\) and the number of observations per group \(n_{i} = 25\), while allowing the number of distinct values \(k\) to vary. To enhance the variability of \(k\) across simulated datasets, observations are sampled from a hierarchical Dirichlet process with random concentration parameters \(\alpha \sim \text{Gamma}(5)\) and \(\alpha_0 \sim \text{Gamma}(3)\). Consistently with Section 5.4, results are averaged over \(100\) simulated datasets for each value of \(k\), and execution times are plotted whenever at least \(75\) experiments are completed without errors. Remarkably, numerical overflows in the computation of the multivariate Stirling numbers appear more frequently for smaller values of \(k\); indeed, since the number of observations \(n = 250\) is fixed, a smaller \(k\) entails a larger number of observations for each distinct value. As displayed in Figure 14, the execution time for the MCMC algorithm grows linearly with \(k\): this is in line with the fact that the state of the Markov chain has dimension \(2k+1\). The exact sampler and collapsed Gibbs sampler also display a nearly linear dependence on the number of distinct values, with a smaller slope compared with the MCMC approach. In contrast, the execution time for the CRF-based Gibbs sampler (not shown) is essentially constant.

11 Eliciting the dependence structure↩︎

Normalized hierarchical CRVs induce dependence between the marginal random probabilities, which in turn regulates the borrowing of information across different groups. Hierarchical models are a natural way to induce positive dependence, especially in a Bayesian setting. However, as highlighted in [15], the elicitation of the dependence is more difficult with respect to other proposals because model parameters typically affect both the marginal distribution and the dependence across groups. For this reason, [15] propose two kinds of (weak) flexibility for partially exchangeable models: (i) for every \(\rho \in [0,1]\), there exists a specification of the parameters such that \(\mathrm{Corr}(\tilde{P}_i(A), \tilde{P}_j(A)) = \rho\), where \(A\) is a fixed set s.t. \(E(\tilde{P}_i(A)) \neq 0,1\); (ii) for every \(\rho \in [0,1]\) and for every fixed value of the marginal mean \(E(\tilde{P}_i(A))\) and variance \(\mathrm{Var}(\tilde{P}_i(A))\), there exists a specification of the parameters s.t. \(\mathrm{Corr}(\tilde{P}_i(A), \tilde{P}_j(A)) = \rho\). The expressions in Example 3 show that the normalized gamma-gamma hCRV achieves the flexibility of second kind. In this section, we investigate the role of parameter values in the borrowing of information, for fixed values of the marginal mean and variance.

Let \(\boldsymbol{\tilde{\mu}}\) be a gamma-gamma hCRV as in Example 1 and let \(A\) be a Borel set s.t. \(P_0(A) \neq 0,1\). From Example 3, it is apparent that correlation is not affected by the mean \(E(\tilde{P}_i(A)) = P_0(A)\) but is deeply related to the variance. Indeed, one has \[\begin{align} \mathrm{Corr}( \tilde{P}_i(A), \tilde{P}_j(A)) & = \left(1 + \frac{\alpha_0}{\alpha} e^{1/\alpha} E_{\alpha_0}\big(1/\alpha\big)\right)^{-1}, \\ \mathrm{Var}( \tilde{P}_i(A) ) & = \left(1 + \frac{\alpha_0}{\alpha} e^{1/\alpha} E_{\alpha_0}\big(1/\alpha \big)\right) \frac{ P_0(A)(1-P_0(A))}{1+\alpha_0}. \end{align}\] We refer to Table ¿tbl:table:limits? for some limiting behaviours, which are compared with those of the HDP. Since \(\alpha\) and \(\alpha_0\) impact both the variance and the dependence structure, it may be difficult to elicit them in practice. From the point of view of the practitioner, it is certainly more intuitive to elicit mean, variance and correlation independently. Hence, one can choose \(\sigma^2, \rho \in (0,1)\) such that \(\mathrm{Var}(\tilde{P}_i(A)) = \sigma^2 \, P_0(A)(1-P_0(A))\) and \(\mathrm{Corr}(\tilde{P}_i(A), \tilde{P}_j(A)) = \rho\), and then find the corresponding values of \(\alpha\) and \(\alpha_0\) by solving the system of non-linear equations in Section 6, namely \[\rho(1+\alpha_0/\alpha \, e^{1/\alpha} E_{\alpha_0}(1/\alpha)) - 1 = 0, \qquad \sigma^2 (1 + \alpha_0) - 1/\rho = 0,\] Figure 15 displays the values of \(\alpha,\alpha_0\) that correspond to a range of values for \(\sigma^2\) and \(\rho\), along with the corresponding analysis for the HDP.

To visualize the effect of the borrowing, we consider the same simulated dataset as in Section 6, namely \(d = 3\) groups of independent Poisson observations, each of size \(n_{i}=10\), with means equal to \(2\), \(3\) and \(4\). Figure 16 displays the expected values of the posterior random means \(E \big( \int x \, \mathrm{d}\tilde{P}_i(x) \bigm\vert \boldsymbol{X}\big)\), which coincide with the means of the predictive distributions for the three groups, as we vary the prior pairwise correlation coefficient \(\rho \in (0,1)\) and the prior variance through \(\sigma^2\in (0,1)\), keeping the base measure \(P_0 = N(10,1)\) fixed. Posterior samples are obtained using the exact algorithm. For a fixed variance \(\sigma^2\), a higher correlation \(\rho\) induces more borrowing, and thus the posterior means are closer to each other. Depending on \(\sigma^2\), the estimates are closer to the prior (which is pushing them towards higher values, being centered at \(10\)), or to their empirical means. On the other hand, for a fixed correlation \(\rho\), a higher variance \(\sigma^2\) reduces the weight of the prior, and thus pushes the estimates towards their empirical means. Interestingly, lower values of \(\rho\) can induce estimates that are closer to each other than higher values of \(\rho\), depending on the value of \(\sigma^2\). Indeed, low values of \(\sigma^2\) will force estimates to be closer to the prior, and thus also closer to each other; e.g., compare \(\rho = \sigma^2 = 0.1\) vs. \(\rho = 0.5\) and \(\sigma^2 = 0.9\). This is the effect of the shrinkage, which can sometimes be difficult to distinguish from the borrowing of information, especially if the prior mean is close to the grand mean of the observations.

References↩︎