Determinantal point process mixtures via spectral density approach

TL;DR

This paper introduces a spectral density approach to DPP mixture models that encourages well-separated location parameters, extends to covariate incorporation, and develops Bayesian inference with empirical validation.

Contribution

It proposes a general spectral representation for DPP mixture models, allowing flexible approximation and extension to covariate-dependent settings, with comprehensive Bayesian inference.

Findings

Effective separation of mixture components demonstrated

Model extensions incorporate covariate information

Bayesian inference performs well in simulations and data examples

Abstract

We consider mixture models where location parameters are a priori encouraged to be well separated. We explore a class of determinantal point process (DPP) mixture models, which provide the desired notion of separation or repulsion. Instead of using the rather restrictive case where analytical results are available, we adopt a spectral representation from which approximations to the DPP intensity functions can be readily computed. For the sake of concreteness the presentation focuses on a power exponential spectral density, but the proposed approach is in fact quite general. We later extend our model to incorporate covariate information in the likelihood and also in the assignment to mixture components, yielding a trade-off between repulsiveness of locations in the mixtures and attraction among subjects with similar covariates. We develop full Bayesian inference, and explore model properties and posterior behavior using several simulation scenarios and data illustrations.

Tables

Table 1. Table 1: Prior specification for ( ρ , ν ) 𝜌 𝜈 (\rho,\nu) and K 𝐾 K and posterior summaries for the Galaxy dataset; ( a ρ , b ρ ) subscript 𝑎 𝜌 subscript 𝑏 𝜌 (a_{\rho},b_{\rho}) appear in ( 21 ); here 𝒱 1 subscript 𝒱 1 \mathcal{V}_{1} is { 1 , 2 , 5 , 10 , 20 } 1 2 5 10 20 \{1,2,5,10,20\} and 𝒱 2 = { 0.5 , 1 , 2 , 3 , 5 , 10 , 15 , 20 , 30 , 50 } subscript 𝒱 2 0.5 1 2 3 5 10 15 20 30 50 \mathcal{V}_{2}=\{0.5,1,2,3,5,10,15,20,30,50\} .

Test	$\rho$	$\nu$	$𝔼\left(K\right)$	$Var\left(K\right)$	$𝔼\left(K|data\right)$	$Var\left(K|data\right)$	MSE	LPML
1	2	2	2	1.67	6.09	1.10	78.95	-171.72
2	5	10	5.00	7.12	6.07	1.09	78.33	-167.96
3	$a_{\rho }=1,b_{\rho }=1$	2	2.18	1.978	6.10	1.10	73.89	-164.47
4	$a_{\rho }=1,b_{\rho }=1$	10	2.73	2.15	6.11	1.12	74.93	-162.71
5	$a_{\rho }=1,b_{\rho }=1$	discr(${𝒱}_1$)	2.47	2.21	6.06	1.08	74.02	-172.54
6	$a_{\rho }=1,b_{\rho }=1$	discr(${𝒱}_2)$	2.51	2.27	6.10	1.13	76.64	-170.94

Table 2. Table 2: Prior specification for α 𝛼 \alpha and posterior summaries for the Galaxy dataset using the function DPdensity in DPpackage .

Test	$\alpha$	$𝔼\left(K|data\right)$	$Var\left(K|data\right)$	MSE	LPML
7	$gamma\left(a_0=0.5,b_0=1\right)$	6.166	1.549	62.703	-151.797
8	$0.8$	5.936	1.25	61.255	-151.146
9	0.45	4.371	1.142	139.659	-169.978
10	$gamma(4,2)$	7.271	1.594	36.708	-149.258

Table 3. Table 3: Root MSE for estimating 𝔼 ​ ( y ∣ x 1 , x 2 , x 3 ) 𝔼 conditional 𝑦 subscript 𝑥 1 subscript 𝑥 2 subscript 𝑥 3 \mathbb{E}(y\mid x_{1},x_{2},x_{3}) for 12 combinations of covariates ( x 1 , x 2 , x 3 ) subscript 𝑥 1 subscript 𝑥 2 subscript 𝑥 3 \left(x_{1},x_{2},x_{3}\right) and P ​ P ​ M x 𝑃 𝑃 subscript 𝑀 𝑥 PPM_{x} as competing model of reference (compare also the results in Table 1 of Müller et al., 2011 ).

$x_1$	$x_2$	$x_3$	DPP	$PP{M}_X$
-1	0	0	6.1	7.9
0	0	0	6.7	3.9
1	0	0	7.2	2.8
-1	1	0	6.5	5.4
0	1	0	6.5	4.6
1	1	0	6.8	4.0
-1	0	1	6.8	6.1
0	0	1	6.1	4.2
1	0	1	5.7	4.5
-1	1	1	5.9	9.5
0	1	1	6.6	8.3
1	1	1	5.8	6.2
avg			6.4	5.6

Table 4. Table 4: Comparison with competitors for the simulated dataset with covariates: best values according to each index are in bold. DPPx denotes our model, while LDDP is the linear dependent Dirichlet process mixture, P ​ P ​ M x 𝑃 𝑃 subscript 𝑀 𝑥 PPM_{x} is the product partition model with covariates, and LDTFP is the linear dependent tailfree process model.

	DPPx	LDDP	$PP{M}_X$	LDTFP
Root $MS{E}_{train}$	324.531	304.742	278.395	304.374
Root $MS{E}_{test}$	216.675	215.1694	217.2459	212.761
$LPM{L}_{train}$	-871.8	-902.2295	-873.1671	-901.465

Table 5. Table 5: Prior specification for β k subscript 𝛽 𝑘 \beta_{k} ’s and σ k 2 superscript subscript 𝜎 𝑘 2 \sigma_{k}^{2} ’s parameters and posterior summaries for the Biopics dataset; m 𝑚 m and v 𝑣 v are prior mean and variance, respectively, of σ k 2 superscript subscript 𝜎 𝑘 2 \sigma_{k}^{2} . Posterior mean and variance of the number K 𝐾 K of mixture components are in the fifth and sixth columns, respectively, while the last two columns report MSE and LPML, respectively.

Test	$\varphi$	$m$	$v$	$𝔼\left(K\mid data\right)$	sd($K\mid$data)	MSE	LPML
$A$	50	5	1	4.49	1.10	1126.32	-960.89
$B$	200	5	10	4.45	1.19	983.55	-954.55
$C$	50	3	$+\infty$	5.66	1.27	501.22	-918.74
$D$	200	10	5	4.21	1.33	1805.83	-980.61
$E$	100	2	1	5.31	1.21	564.26	-935.56
$F$	200	2	10	5.51	1.26	557.44	-925.22

Table 6. Table 6: Posterior summaries for the tests on the Biopics dataset under a linear dependent Dirichlet process mixture.

Case	$𝔼\left(K\mid data\right)$	sd$\left(K\mid data\right)$	MSE	LPML
$G$	2.95	1.03	1282.49	-937.51
$H$	3.56	2.36	682.98	-914.00

Table 7. Table S.7: Prior specification for ( ρ , ν ) 𝜌 𝜈 (\rho,\nu) and the corresponding mean and variance induced on K 𝐾 K (top). Hyperparameters ( a ρ , b ρ ) subscript 𝑎 𝜌 subscript 𝑏 𝜌 (a_{\rho},b_{\rho}) appear in ( 21 ), while 𝒱 1 = { 1 , 2 , 5 , 10 , 20 } subscript 𝒱 1 1 2 5 10 20 \mathcal{V}_{1}=\{1,2,5,10,20\} and 𝒱 2 = 𝒱 1 ∪ { 0.5 , 3 , 15 , 30 , 50 } subscript 𝒱 2 subscript 𝒱 1 0.5 3 15 30 50 \mathcal{V}_{2}=\mathcal{V}_{1}\cup\{0.5,3,15,30,50\} . Posterior summaries for the simulated dataset from a mixture with 8 components are in the bottom subtable.

	Prior specification
Test	$\rho$	$\nu$	$𝔼\left(K\right)$	$V\left(K\right)$
$S_0$	9.00	1	8.98	45.12
$S_1$	9	10	9	23.05
$S_2$	$a_{\rho }=1,b_{\rho }=1$	1	1.94	1.99
$S_3$	$a_{\rho }=1,b_{\rho }=1$	2	2.18	1.99
$S_4$	$a_{\rho }=1,b_{\rho }=1$	10	2.74	2.17
$S_5$	$a_{\rho }=1,b_{\rho }=1$	discr(2,5,20)	2.52	2.11
$S_6$	$a_{\rho }=1,b_{\rho }=1$	discr(${𝒱}_1$)	2.45	2.18
$S_7$	$a_{\rho }=1,b_{\rho }=1$	discr(${𝒱}_2$)	2.5	2.25

Table 8. Table S.8: Prior mean and variance of K 𝐾 K and posterior summaries for the Galaxy dataset with Whittle-Matérn (top) and generalized Cauchy (bottom) spectral densities.

	Whittle-Matérn
Test	$𝔼\left(K\right)$	Var$\left(K\right)$	$𝔼\left(K\mid data\right)$	Var$\left(K\mid data\right)$	MSE	LPML
$\left(i\right)$	10.21	17.29	6.07	1.09	73.67	-167.22
$\left(ii\right)$	2.09	2.15	6.08	1.09	73.89	-167.68
$\left(iii\right)$	3.53	9.87	6.07	1.10	75.80	-167.33

Table 9. Table S.9: Prior specification for the Air quality index dataset. The scale parameter ϕ italic-ϕ \phi appears in the g 𝑔 g -prior specification of ( 18 ), while m 𝑚 m and v 𝑣 v denote prior mean and variance, respectively, of σ k 2 subscript superscript 𝜎 2 𝑘 \sigma^{2}_{k} as in ( 16 ).

Test	$\varphi$	$m$	$v$	mean($K$)	sd($K$)	MSE	LPML
$A{Q}_1$	1000	2	1	6.999	1.469	861.048	-1101.013
$A{Q}_2$	500	10	$+\infty$	5.192	1.143	870.689	-1235.988
$A{Q}_3$	1000	0.1	1	9.045	2.243	840.0685	-1071.931
$A{Q}_4$	500	5	$+\infty$	7.811	2.665	835.596	-1160.631

Table 10. Table S.10: Posterior summaries for the Air Quality Index dataset under the linear dependent Dirichlet process mixture for two different prior specifications.

Test	$𝔼\left(K\mid data\right)$	Var$\left(K\mid data\right)$	MSE	LPML
$A{Q}_5$	5.14	0.38	827.72	-1100.73
$A{Q}_6$	5.03	0.16	827.04	-1100.06

Equations

$$y_i\mid k,\mathbf{\theta },\mathbf{\pi }\stackrel{\mathrm{ind}}{\sim }\sum _{j=1}^k{\pi }_{j}f(y_i\mid {\theta }_j),i=1,\dots ,n,$$

$$\mathbb{E}\left(\sum _{x_1,\dots ,x_n\in X}^{\ne }h(x_1,\dots ,x_n)\right)={\int }_{B^n}{\rho }^{(n)}(x_1,\dots ,x_n)h(x_1,\dots ,x_n)d{x}_1\cdots d{x}_n,$$

$${\rho }^{(n)}(x_1,\dots ,x_n)=det[C](x_1,\dots ,x_n),(x_1,\dots ,x_n)\in B^n,n=1,2,\dots ,$$

$$C(x,y)=\sum _{k=1}^{+\infty }{\lambda }_k^S{\varphi }_k(x){\varphi }_k(y),(x,y)\in S\times S,S\text{ compact subset},$$

$$N(S)\stackrel{d}{=}\sum _{k=1}^{+\infty }{B}_k,\mathbb{E}(N(S))=\sum _{k=1}^{+\infty }{\lambda }_k^S,\mathrm{Var}(N(S))=\sum _{k=1}^{+\infty }{\lambda }_k^S(1-{\lambda }_k^S),$$

$$f(\{x_1,\dots ,x_n\})={\mathrm{e}}^{|S|-D_S}det[\tilde{C}](x_1,\dots ,x_n),n=1,2,\dots ,$$

$$\tilde{C}(x,y)=\sum _1^{+\infty }\frac{{\lambda }_k^S}{1-{\lambda }_k^S}{\varphi }_k(x){\varphi }_k(y),x,y\in S.$$

$$C_0(x)={\int }_{\mathbb{R}{}^{{}^d}}\phi (y)\cos (2\pi x\cdot y)\mathrm{d}y,x\in \mathbb{R}{}^{{}^d}$$

$$f^{app}(\{x_1,\dots ,x_n\})={\mathrm{e}}^{|S|-D_{app}}det[{\tilde{C}}_{app}](x_1,\dots ,x_n),\{x_1,\dots ,x_n\}\subset S,$$

$${\tilde{C}}_{app}(x,y)=C_{app,0}(x-y)=\sum _{k\in {\mathbb{Z}}^d}\left[\frac{\phi (k)}{1-\phi (k)}\right]\cos (2\pi k\cdot (x-y)),x,y\in S,$$

$$D_{app}=\sum _{k\in {\mathbb{Z}}^d}\log \left(1+\frac{\phi (k)}{1-\phi (k)}\right).$$

$$Re\left({\int }_S{C}_0(y){\mathrm{e}}^{-2\pi ik\cdot y}\mathrm{d}y\right)\approx \phi (k):=Re\left({\int }_{\mathbb{R}{}^{{}^d}}{C}_0(y){\mathrm{e}}^{-2\pi ik\cdot y}\mathrm{d}y\right).$$

$$f^{app}(\{x_1,\dots ,x_n\})=|R{|}^{-n}{\mathrm{e}}^{|R|-|S|}{f}_Y^{app}(T(\{x_1,\dots ,x_n\})),\{x_1,\dots ,x_n\}\subset R.$$

$$\phi (x;\rho ,\nu )=s^d\exp \left\{-{\left(\frac{s}{\sqrt{\pi }}\right)}^{\nu }{\left(\frac{\Gamma (\frac{d}{2}+1)}{\Gamma (\frac{d}{\nu }+1)}\right)}^{\nu d}{\rho }^{\nu d}\Vert x{\Vert }^{\nu }\right\},\rho ,\nu >0,$$

$$\phi (x;\rho ,\nu )=s^d\exp \left\{-\frac{s^2{\rho }^{2d}}{\sqrt{\pi }}\Vert x{\Vert }^2\right\},\rho >0,$$

$$y_i\mid s_i=k,\{{\mu }_k\},\{{\sigma }_k^2\},K$$

$$X=\{{\mu }_1,{\mu }_2,\dots ,{\mu }_K,K\}$$

$$(\rho ,\nu )$$

$$p(s_i=k)$$

$$w_1,\dots ,w_K\mid K$$

$${\sigma }_k^2\mid K$$

$$p(s_i=k)=w_k(x_i)$$

$${\beta }_2,\dots ,{\beta }_K\mid K$$

$$y_i\mid s_i=k,x_i,\{{\mu }_k\},\{{\sigma }_k^2\},K$$

$$({\gamma }_1,{\sigma }_1^2),\dots ,({\gamma }_K,{\sigma }_K^2)|K$$

$$\rho \mid \nu \stackrel{d}{=}M(s,\epsilon ,\nu )+{\rho }_0,{\rho }_0\sim gamma(a_{\rho },b_{\rho }).$$

$$M(s,\epsilon ,\nu )=\frac{2s\Gamma (1/\nu +1){\pi }^{1/2}}{\Gamma (3/2){\left(\log \left(\frac{s}{\epsilon }\right)\right)}^{1/\nu }}.$$

$$\rho \stackrel{d}{=}\sqrt{\frac{\pi }{\log (\frac{1}{2\epsilon })}}+{\rho }_0,{\rho }_0\sim gamma(a_{\rho },b_{\rho }).$$

$$p(s_i=k\mid rest)\propto w_k\mathcal{N}\left(y_i;{\mu }_k,{\sigma }_k^2\right).$$

$${\sigma }_k^2\mid rest\sim inv-gamma\left(a_0+\frac{n_k}{2},b_0+\frac{1}{2}\sum _{i:s_i=k}{\left(y_i-{\mu }_k\right)}^2\right),k=1,\dots ,K.$$

$$p({\mu }_1,\dots ,{\mu }_K\mid rest)\propto det[\tilde{C}](\{{\tilde{\mu }}_1,\dots ,{\tilde{\mu }}_K\};\rho ,\nu )\prod _{k=1}^K\prod _{i:s_i=k}\mathcal{N}\left(y_i;{\mu }_k,{\sigma }_k^2\right)$$

$$\propto \prod _{k=1}^K\left\{\left(\tilde{C}({\tilde{\mu }}_k,{\tilde{\mu }}_k)-b{\tilde{C}}_{-k}^{-1}{b}^T\right)\prod _{i:s_i=k}\mathcal{N}\left(y_i;{\mu }_k,{\sigma }_k^2\right)\right\},$$

$$p(\rho ,\nu \mid rest)\propto det\left(\tilde{C}\right)[{\tilde{\mu }}_1,\dots ,{\tilde{\mu }}_K,\rho ,\nu ]\exp \left(-\sum _{k=-N}^N\log (1+\tilde{\phi }(k;\rho ,\nu ))\right)\pi (\rho ,\nu ).$$

$$w_{j_1}^{new}=w_{j_1}+w_{j_2}$$

Full text

Keywords: Density Estimation; Nonparametric Regression; Repulsive Mixtures; Reversible Jumps

1 Introduction

Mixture models are an extremely popular class of models, that have been successfully used in many applications. For a review, see, e.g. Frühwirth-Schnatter (2006). Such models are typically stated as

[TABLE]

where $\bm{\pi}=(\pi_{1},\ldots,\pi_{k})$ are constrained to be nonnegative and sum up to 1, $\bm{\theta}=(\theta_{1},\ldots,\theta_{k})$, and $1\leq k\leq\infty$, with $k=\infty$ corresponding to a nonparametric model. A common prior assumption is that $\bm{\pi}\sim\operatorname{Dirichlet}(\delta_{1},\ldots,\delta_{k})$ and that the components of $\bm{\theta}$ are drawn i.i.d. from some suitable prior $p_{0}$. However, the weights $\bm{\pi}$ may be constructed differently, e.g. using a stick-breaking representation (finite or infinite), which poses a well-known connection with more general models, including nonparametric ones. See, e.g., Ishwaran and James (2001) and Miller and Harrison (2017). A popular class of Bayesian nonparametric models is the Dirichlet process mixture (DPM) model, introduced in Ferguson (1983) and Lo (1984). It is well-known that this class of mixtures usually overestimates the number of clusters, mainly because of the “rich gets richer” property of the Dirichlet process. By this we mean that both prior and posterior distributions are concentrated on a relatively large number of clusters, but a few are very large, and the rest of them have very small sample sizes. Mixture models may even be inconsistent; see Rousseau and Mengersen (2011), where concerns about over-fitted mixtures are illustrated, and Miller and Harrison (2013), for inconsistency features of DPMs.

Despite their success, mixture models like (1) tend to use excessively many mixture components. As pointed out in Xu et al. (2016), this is due to the fact that the component-specific parameters are a priori i.i.d., and therefore, free to move. This motivated Petralia et al. (2012), Fúquene et al. (2016) and Quinlan et al. (2017) to explicitly define joint distributions for $\bm{\theta}$ having the property of repulsion among its components, i.e. that $p(\theta_{1},\ldots,\theta_{k})$ puts higher mass on configurations such that components are well separated. For a different approach, via sparsity in the prior, see Malsiner-Walli et al. (2016).

Xu et al. (2016) explored a similar way to accomplish separation of mixture components, by means of a Determinantal Point Process (DPP) acting on the parameter space. DPPs have recently received increased attention in the statistical literature (Lavancier et al., 2015). DPPs are point processes having a product density function expressed as the determinant of a certain matrix constructed using a covariance function evaluated at the pairwise distances among points, in such a way that higher mass is assigned to configurations of well-separated points. We give details below. DPPs have been used in a number of different modeling efforts. Bardenet and Titsias (2015) and Affandi et al. (2014) applied DPPs to model spatial patterns of nerve fibers in diabetic patients, a basic motivation being that such fibers become more clustered as diabetes progresses. The latter discussed also applications to image search, showing how such processes could be used to study human perception of diversity in different image categories. Similarly, Kulesza and Taskar (2012) show how DPPs can be applied to various problems such as finding diverse sets of high-quality search results, building informative summaries by selecting diverse sentences from documents, modeling non-overlapping human poses in images or video, and automatically building timelines of important news stories.

We discuss full Bayesian inference for a class of mixture densities where the locations follow stationary DPPs. The approach is based on the spectral representation of the covariance function defining the determinant as the joint distribution of component-specific parameters. While our methods can be used with any such valid spectral representation, we focus for the sake of concreteness on the power exponential spectral case; see examples with different spectral densities in the on-line Supplementary Material. This particular specification allows for flexible repulsion patterns, and we discuss how to set up different types of prior behavior, shedding light on the practical use of our approach in that particular scenario. We discuss applications in the context of synthetic and real data applications. We later extend our model to incorporate covariate information in the likelihood and also in the assignment to mixture components. In particular, subjects with similar covariates are a priori more likely to co-cluster, just as in mixtures of experts models (see, e.g., McLachlan and Peel, 2005), where weights are defined as normalized exponential functions.

We implement a reversible jump (RJ) MCMC posterior simulation scheme for each of the models we propose. In all cases, we estimate clusters for the subjects in the sample, by considering the partition that minimizes the posterior expectation of Binder’s loss function (Binder, 1978) under equal misclassification costs. This is a common choice in the applied Bayesian nonparametric literature (Lau and Green, 2007).

The first mixture model we propose is as in Xu et al. (2016). They consider, as a prior for the location points, a particular case of DPPs, called L-ensambles. These are defined as a direct generalization from the same definition but for the finite state space case. Our proposal does not use L-ensembles. Instead, we stick to the more general case, following the spectral approach discussed in Lavancier et al. (2015). Although we limit ourselves to the case of isotropic DPPs, inhomogeneous DPPs can be obtained by transforming or thinning a stationary process; see Lavancier et al. (2015), Section 3 and Appendix A. A crucial point in our models and algorithms is the DPP density expression, which is only defined for DPPs restricted to compact subsets $S$ of the state space, with respect to the unit rate Poisson process. When this density exists, it explicitly depends on $S$. A sufficient condition for the existence of the density is that all the eigenvalues of the covariance function, restricted to $S$, are smaller than 1. We follow the spectral approach and assume that the covariance function defining the DPP has a spectral representation. A basic motivation for our choice is that conditions for the existence of a density become easier to check. We review here the basic theory on DPPs, making an effort to be as clear and concise as possible in the presentation of our subsequent models.

We acknowledge that the RJ scheme for our first mixture model is as in Xu et al. (2016). In both cases the algorithm requires computing the DPP density with respect to the unit rate Poisson process. We explain how to adapt the calculations to our case, and discuss the need to restrict the process to (any) compact subset. When extending the model to incorporate covariate information in both, likelihood and prior assignment to mixture components, the RJ MCMC algorithm requires modifications, as discussed below.

The rest of this article is organized as follows. Section 2 presents notation and theoretical background necessary to understand how DPP mixture models are constructed. Section 3 illustrates the covariate-dependent extension. Here we build on regular mixture models, incorporating covariate dependence in the mixture weights and optionally in the likelihood, which still allows for repulsion among components after correcting for the regression effect. Section 4 presents results from a simulation study and for the Galaxy dataset, while data applications are discussed in Section 5. We conclude in Section 6 with final comments and discussion. The on-line Supplementary Material contains a description of the two RJ MCMC algorithms, additional illustrative examples on simulated and real data, and supplemental figures.

2 Using DPPs to induce repulsion

We review here the basic theory on DPPs to the extent required to explain our mixture model. We use the same notation as in Lavancier et al. (2015), where further details on this theory may be found.

2.1 Basic theory on DPPs

Let $B\subseteq\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$; we mainly consider the cases $B=\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$ and $B=S$, a compact subset in $\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$. By $X$ we denote a simple locally finite spatial point process defined on $B$, i.e. the number of points of the process in any bounded region is a finite random variable, and there is at most one point at any location. See Daley and Vere-Jones (2003; 2007) for a general presentation on point processes. The class of determinantal point processes we consider is defined in terms of their moments, expressed by their product density functions $\rho^{(n)}:B^{n}\rightarrow[0,+\infty)$, $n=1,2,\ldots$. Intuitively, for any pairwise distinct points $x_{1},\ldots,x_{n}\in B$, $\rho^{(n)}(x_{1},\ldots,x_{n})dx_{1}\cdots dx_{n}$ is the probability that $X$ has a point in an infinitesimal small region around $x_{i}$ of volume $dx_{i}$, for each $i=1,\ldots,n$. More formally, $X$ has $n$-th order product density function $\rho^{(n)}:B^{n}\rightarrow[0,+\infty)$ if this function is locally integrable (i.e. $\int_{S}|\rho^{(n)}(x)|\mathrm{d}x<+\infty$ for any compact $S$) and, for any Borel-measurable function $h:B^{n}\rightarrow[0,+\infty)$,

[TABLE]

where the $\neq$ sign over the summation means that $x_{1},\ldots,x_{n}$ are pairwise distinct. See also Møller and Waagepetersen (2007). Let us consider a covariance function $C:B\times B\rightarrow{\mathbb{R}}$. The theory recalled here is valid also for complex-valued covariance functions, but such extensions are not needed in what follows.

Definition**.**

A simple locally finite spatial point process $X$ on $B$ is called a determinantal point process with kernel $C$ if its product density functions are

[TABLE]

where $[C](x_{1},\ldots,x_{n})$ is the $n\times n$ matrix with entries $C(x_{i},x_{j})$. We write $X\sim DPP_{B}(C)$; when $B=\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$ we write $X\sim DPP(C)$.

Remark**.**

If $A$ is a Borel subset of $B$, then the restriction $X_{A}:=X\cap A$ of $X$ to $A$ is a DPP with kernel given by the restriction of $C$ to $A\times A$.

By Theorem 2.3 in Lavancier et al. (2015), first proved by Macchi (1975), such DPP’s exist under the two following conditions:

• •

$C$ is a continuous covariance function; hence, by Mercer’s Theorem,

[TABLE]

where $\{\lambda_{k}^{S}\}$ and $\{\phi_{k}(x)\}$ are the eigenvalues and eigenfunctions of $C$ restricted to $S\times S$, respectively.

• •

$\lambda_{k}^{S}\leq 1$ for all compact $S$ in $\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$ and all $k$.

Formula (2.9) in Lavancier et al. (2015) reports the distribution of the number $N(S)$ of points of $X$ in $S$, for any compact $S$:

[TABLE]

where $B_{k}\stackrel{{\scriptstyle\mathrm{ind}}}{{\sim}}Be(\lambda_{k}^{S})$, i.e. the Bernoulli random variable with mean $\lambda_{k}^{S}$. When restricted to any compact subset $S$, the DPP has a density with respect to the unit rate Poisson process which, when $\lambda_{k}^{S}<1$ for all $k=1,2,\ldots$, has the following expression:

[TABLE]

where $|S|=\int_{S}\mathrm{d}x$, $D_{S}=-\sum_{1}^{+\infty}\log(1-\lambda_{k}^{S})$ and

[TABLE]

When $n=0$ the density (as well as the determinant) is defined to be equal to 0. See Møller and Waagepetersen (2007) for a thorough definition of absolute continuity of a spatial process with respect to the unit rate Poisson process. However, note that from the first part of (3) we have $\mathbb{P}(N(S)=0)=\prod_{k=1}^{+\infty}(1-\lambda_{k}^{S})$; this probability could be positive due to the assumption $\lambda_{k}^{S}<1$ for all $k=1,2,\ldots$.

From now on we restrict our attention to stationary DPP’s, that is, when $C(x,y)=C_{0}(x-y)$, where $C_{0}\in L^{2}(\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}})$ is such that its spectral density $\varphi$ exists, i.e.

[TABLE]

and $x\cdot y$ is the scalar product in $\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$. If $\varphi\in L^{1}(\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}})$ and $0\leq\varphi\leq 1$, then the $DPP(C)$ process exists. Summing up, the distribution of a stationary DPP can be assigned by its spectral density; see Corollary 3.3 in Lavancier et al. (2015).

To explicitly evaluate (4) over ${\displaystyle S=\left[-\dfrac{1}{2},\dfrac{1}{2}\right]^{d}}$, we approximate $\widetilde{C}$ as suggested in Lavancier et al. (2015). In other words, we approximate the density of $X$ on $S$ by

[TABLE]

where

[TABLE]

[TABLE]

The approximation $C(x,y)\approx C_{app,0}(x-y)$ ($x-y\in S$) follows because the exact Fourier expansion of $C_{0}(x-y)$ in $S$ is as in (7) with the real part of $\int_{S}C_{0}(y)\mathrm{e}^{-2\pi ik\cdot y}\mathrm{d}y$ instead of $\varphi(k)$; if we assume $C_{0}$ such that $C_{0}(t)\approx 0$ for $t\not\in S$, then

[TABLE]

See also Lavancier et al. (2015), Section 4.1. See Figure S.9 in the Supplementary Material, where we display an example of the function $C_{0}$.

When $S=R$ is a rectangle in $\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$, we can always find an affine transformation $T$ such that $T(R)=S=\left[-\dfrac{1}{2},\dfrac{1}{2}\right]^{d}$. Define $Y=T(X)$. If $f^{app}_{Y}$ is the approximate density of $Y$ as in (6), we can then approximate the density of $X_{R}$ by

[TABLE]

In practice, the summation over ${\mathbb{Z}}^{d}$ in (7) above is truncated to ${{\mathbb{Z}}_{N}}^{d}$, where ${\mathbb{Z}}_{N}:=\{-N,-N+1,\ldots,0,\ldots,N-1,N\}$; see further details in Lavancier et al. (2015), Section 4.3.

One particular example of spectral density that we found useful is

[TABLE]

for fixed $s\in(0,1)$ (e.g. $s=\frac{1}{2}$) and $\|x\|$ is the Euclidean norm of $x\in\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$. This function is the spectral density of a power exponential spectral model (see (3.22) in Lavancier et al. (2015) when $\alpha=s\ \alpha_{max}(\rho,\nu)$). In this case, we write $X\sim PES-DPP(\rho,\nu)$. The corresponding spatial process is isotropic. When $\nu=2$, the spectral density is

[TABLE]

corresponding to the Gaussian spectral density.

Remark**.**

Note that $\varphi(x;\rho,\nu)<1$ when $0<s<1$ for any $x\in\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$, $\rho,\nu>0$, so that $X_{S}$ has a density as described in (4). Figure S.8 in the Supplementary material shows a plot of the power exponential spectral density (10) for different values of parameters $\rho$, $\nu$. Note that $\nu$ controls the shape of $\varphi(x;\rho,\nu)$, which ranges from a slowly decreasing function of $x$ to an indicator function. On the other hand $\rho$ plays the role of a centering parameter, with higher values retarding the decay speed of $\varphi(x;\rho,\nu)$. As discussed above, knowledge about the spectral density is all that is needed for the approximations to work. Moreover, even if the analytic expression of $C(x,y),(x,y)\in\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}\times\mathbb{R}{\scriptscriptstyle{{}^{{}^{d}}}}$ is known, as in, e.g. (10), we still need to compute the eigenvalues and eigenfunctions of $C$ restricted to $S\times S$ for any compact $S$, and this may be analytically impossible. A potential disadvantage derived from this is that parameter interpretation in the spectral domain becomes unclear. A possible way to alleviate this difficulty consists of conducting extensive simulations that help understanding the role of such parameters. We do that for the case of (10); see Section 4.

Typically, DPPs have been used to make inference mostly on spatial data (including images, videos and point patterns related to towns, trees and cells); see, for instance, Shirota and Gelfand (2016) who describe an approximate Bayesian computation method to fit DPPs to spatial point pattern data. Historically, the first paper where DPPs were adopted as a prior for statistical inference in mixture models is Affandi et al. (2013).

The statistical literature includes a number of papers illustrating theoretical properties for estimators of DPPs from a non-Bayesian viewpoint. Biscio and Lavancier (2017) study asymptotic properties of minimum contrast estimators for parametric stationary determinantal point processes. Biscio and Lavancier (2016) quantify the repulsiveness of a DPP, based on its second-order properties, and address the question of how repulsive a stationary DPP can be. Moreover, Bardenet and Titsias (2015) derive bounds on the likelihood of a DPP, and Kulesza and Taskar (2012) provide an introduction to DPPs, focusing on the intuitions, algorithms, and extensions that are most relevant to the machine learning community.

2.2 The mixture model with repulsive means

To deal with limitations of model (1) or DPMs, we consider repulsive mixtures. Our aim is to estimate a random partition of the available subjects, and we want to do so using “few” groups. By repulsion we mean that cluster locations are a priori encouraged to be well separated, thus inducing fewer clusters than if they were allowed to be independently selected. We start from parametric densities $f(\cdot;\theta)$, which we take to be Gaussian, and assume that the collection of location parameters follows a DPP. We specify a hierarchical model that achieves the goals previously described. Concretely, we propose:

[TABLE]

where the PES-DPP$(\rho,\nu)$ assumption (12) is regarded as a default choice that could be replaced by any other valid DPP alternative. We note that, as stated, the prior model may assign a positive probability to the case $K=0$. This case of course makes no sense from the viewpoint of the model described above. Nevertheless, we adopt the working convention of redefining the prior to condition on $K\geq 1$, i.e., truncating the DPP to having at least one point. In practice, the posterior simulation scheme later described simply ignores the case $K=0$, which produces the desired result. Note also that we have assumed prior independence among blocks of parameters not involving the locations $\mu_{k}$.

Model (11)-(16) is a DPP mixture model along the lines proposed in Xu et al. (2016). Indeed, we both use DPPs as priors for location points in the mixture of parametric densities. However, the specific DPP priors are different, as they restrict to a particular case of DPPs (L-ensambles), and choose a Gaussian covariance function for which eigenvalues and eigenfunctions are analytically available. We adopt instead a spectral approach for assigning the prior (12). Similar to Xu et al. (2016), we carry out posterior simulation using a reversible jump step as part of the Gibbs sampler. However, when updating the location points $\mu_{1},\ldots,\mu_{K}$ we refer to formulas (6)-(9). Xu et al. (2016) take advantage of the analytical expressions that we do not have for our case, and that are also unavailable in other possible specific choices of the spectral density. The posterior simulation algorithm we propose for our model is described in Section S.1 in the Supplementary material.

3 Generalization to covariate-dependent models

The methods discussed in Section 2 were devised for density estimation-like problems. We now extend the previous modeling to the case where $p$-dimensional covariates $x_{1},\ldots,x_{n}$ are recorded as well. We do so by allowing the mixture weights to depend on such covariates. In this case, there is a trade-off between repulsiveness of locations in the mixtures and attraction among subjects with similar covariates. We also entertain the case where covariate dependence is added to the likelihood part of the model. Our modeling choice here is akin to mixtures of experts models (see, e.g., McLachlan and Peel, 2005), i.e., the weights are defined by means of normalized exponential function.

Building on the model from Section 2.2, we assume the same likelihood (11) and the DPP prior for $X=\{\mu_{1},\mu_{2},\dots,\mu_{K},K\}$ in (12)-(13), but change (14) and (15) to

[TABLE]

where the $\beta_{1}=0$ assumption is to ensure identifiability. To complete the model, we assume (16) as the conditional marginal for $\sigma^{2}_{k}$; the prior for $(\rho,\nu)$ in (13) is later specified. Here $\beta_{0}\in{\mathbb{R}}^{p}$, and to choose $\Sigma_{0}$, we use a g-prior approach, namely $\Sigma_{0}=\phi\times\left(X^{T}X\right)^{-1}$, where $\phi$ is fixed, typically of the same order of magnitude of the sample size (see Zellner, 1986).

Assuming now (11) on top of (17)-(18) rules out the case of a likelihood explicitly depending on covariates, which instead would generally achieve a better fit than otherwise. Of course, there are many ways in which such dependence may be added. For the sake of concreteness, we assume here a Gaussian regression likelihood, where only the intercept parameters arise from the DPP prior. More precisely, we assume

[TABLE]

where the $\gamma_{k}$’s are $p$-dimensional regression coefficients. The notation in (20) means that $\gamma_{k}\mid\sigma_{k}^{2}\sim\mathcal{N}_{p}(\gamma_{0},\sigma_{k}^{2}\Lambda_{0})$, and $\sigma_{k}^{2}\sim inv-gamma(a_{0},b_{0})$, where $\gamma_{0}\in{\mathbb{R}}^{p}$ and $\Lambda_{0}$ is a covariance matrix. The prior for $\{s_{i}\}$ and $\beta_{j}$’s is given in (17)-(18) as in the previous model. Note that (19) implies that only the intercept term is distributed according to the repulsive prior. Thus, we allow the response mean to be corrected by a linear combination of the covariates with cluster-specific coefficients, with the repulsion acting only on the residual of this regression. The result is a more flexible model than the repulsive mixture (11)-(16). Observe that there is no need to assume the same covariate vector in (19) and (17), but we do so for illustration purposes only.

The Gibbs sampler algorithm employed to carry out posterior inference for this model is detailed in Section S.4 of the Supplementary material. However, it is worth noting that the reversible jump step related to updating the number of mixture components $K$ and the update of the coefficients $\left\{\beta_{2},\beta_{3},\dots,\beta_{K}\right\}$ are complicated by the presence of the covariates. For the $\beta$ coefficients, we resort to a Metropolis-Hastings step, with a multivariate Gaussian proposal centered in the current value. For $K$, we employ an ad hoc Reversible Jump move.

4 Simulated data and reference datasets

Before illustrating the application of our models to specific datasets, we discuss some general choices that apply to all examples. Every run of the Gibbs sampler (implemented in R) produced a final sample size of 5,000 or 10,000 iterations (unless otherwise specified), after a thinning of 10 and initial burn-in of 5,000 iterations. In all cases, convergence was checked using both visual inspection of the chains and standard diagnostics available in the CODA package. Elicitation of the prior for $(\rho,\nu)$ requires some care, as the role of these parameters is difficult to understand. Therefore, an extensive robustness analysis with respect to $\pi(\rho,\nu)$ for those datasets was carried out. See Sections 4.1 and S.2 of the Supplementary material. We point out that an initial prior independence assumption $\pi(\rho,\nu)=\pi(\rho)\ \pi(\nu)$ produced bad mixing of the chain. In particular, when $\rho$ is small with respect to $\nu$, the spectral function $\varphi(\cdot)$ has a very narrow support, concentrated near the origin, forcing the covariance function $\widetilde{C}_{app}(x,y)$ to become nearly constant for $x,y\in S$ and thus producing nearly singular matrices. We next investigated the case $\pi(\rho,\nu)=\pi(\rho\mid\nu)\times\pi_{\nu}(\nu)$, where

[TABLE]

Here, $M(s,\varepsilon,\nu)$ is a constant that is the minimum value of $\rho$ such that $\varphi(2)>\varepsilon$ (here $\varphi(2)$ is a reference value chosen to avoid a small support), and $\varepsilon$ is a threshold value, assumed to be small (0.05, for instance). From Figure S.8 in the Supplementary Material, it is clear that $\varphi(\cdot;\rho,\nu)$ goes to [math] too fast when $\nu$ is small relative to $\rho$. It follows that

[TABLE]

On the other hand, two different choices for $\pi_{\nu}$ were considered: a gamma distribution, which gave a bad chain mixing, and a discrete distribution on $\mathcal{V}_{2}=\{0.5,1,2,3,5,10,15,20,30,50\}$ (or on one of its subsets). In this case, the mixing of the chain was better, but the posterior for $\nu$ did not discriminate among the values in the support. For this reason, in Sections 5.3, 6 and 7, we assume $\nu=2$, $s=1/2$ and

[TABLE]

4.1 Galaxy data

This popular dataset contains $n=82$ measured velocities of different galaxies from six well-separated conic sections of space. Values are expressed in Km/s, scaled by a factor of $10^{-3}$. We set the hyperparameters in this way: for the variance $\sigma^{2}_{k}$ of the components, $\left(a_{0},b_{0}\right)=\left(3,3\right)$ (such that the mean is 1.5 and the variance is $9/4$) and for the weights $\{w_{k}\}$ the Dirichlet has parameter $(1,1,\dots,1)$. The other hyperparameters are modified in the tests, as in Table 1, where we report summaries of interest, such as the prior and posterior mean and variance for the number of components $K$. In addition, we also display the mean squared error (MSE) and the log-pseudo marginal likelihood (LPML) as indices of goodness of fit, defined as $MSE=\sum_{i=1}^{n}(y_{i}-\hat{y_{i}})^{2}$ and $LPML=\sum_{i=1}^{n}\log\left(f(y_{i}\mid y^{(-i)})\right)$, where $\hat{y_{i}}$ is the posterior predictive mean and $f(y_{i}\mid y^{(-i)})$ is the $i-$th conditional predictive ordinate, that is the predictive distribution obtained using the dataset without the $i-$th observation. Figure 1 (left) shows density estimates and the estimated partition of the data, obtained as the partition that minimizes the posterior expectation of Binder’s loss function under equal misclassification costs (see Lau and Green, 2007). The points at the bottom of the plots represent observations, while colors refer to the corresponding cluster. Figure 1 (right) displays the posterior distribution of $K$ for Test 4 and 6 in Table 1.

As a comparison, the same posterior quantities than in Table 1 were computed using the DPM, fitted via the function DPdensity available from DPpackage (Jara et al., 2011); see Table 2. The same prior for $\sigma_{k}^{2}$ as in our model was assumed. Note that $\alpha$ is fixed in Test 8 in such a way that a priori $\mathbb{E}(K)=2$, while for Test 3, $a_{0}$ and $b_{0}$ are chosen so that $\mathbb{E}(\alpha)=2$ and $\operatorname{Var}(\alpha)=1$. Clearly, Test 10 (DPM) shows the best indexes of goodness of fit but at the cost of overestimating the number of clusters. It is well-known that, in general, clustering in the context of nonparametric mixture models as DPMs is strongly affected by the base measure (see, e.g. Miller and Harrison, 2017). Our model, on the other hand, avoids the delicate choice of the base measure leading to more robust estimates of $K$. See also Figure 1 (right) which displays the posterior distribution of $K$ under the DPM mixture and our models.

Finally, we recall that in Section S.3 in the Supplementary material we report some further tests on the Galaxy dataset to show the influence of various choices of spectral density on the inference. We conclude that there is evidence to robustness with respect to the choice of spectral density.

4.2 Simulated data with covariates

We consider the same simulated dataset as in Müller et al. (2011), Section 5.2; the simulation “truth” consists of 12 different distributions, corresponding to different covariate settings (see Figure 1 of that paper). Model (11)-(13), (16)-(18) was fitted to the dataset, assuming $\beta_{0}=0$, $\Sigma_{0}=400\times\left(X^{T}X\right)^{-1}$, $a_{\rho}=1,b_{\rho}=1.2$, and $a_{0},b_{0}$ such that the prior mean of $\sigma^{2}_{k}$ is 50 and variance is 300. Recall also that here we assume $\nu=2$.

As an initial step, inference for the complete dataset (1000 observations) was carried out, yielding a posterior of $K$, not reported here, mostly concentrated over the set $\{8,9,\dots,16\}$, with a mode at 11. Figure S.10 in the Supplementary Material shows posterior predictive distributions for the 12 different reference covariate values, along with 90% credibility intervals. These are in good accordance with the simulation truth (compare Figure 1 in Müller et al., 2011).

To replicate the tests in Müller et al. (2011), a total of $M=100$ datasets of size 200 were generated by randomly subsampling 200 out of the 1000 available observations. Computational burden over multiple repetitions was controlled by limiting the posterior sample sizes to 2,000. Table 3 displays the root MSE for estimating $\mathbb{E}(y\mid x_{1},x_{2},x_{3})$ for each of the 12 covariate combinations defining the different clusters for our model and for the PPMx, as in Table 1 of Müller et al. (2011). The computations also include evaluation of the root MSE and LPML for all the 100 datasets for estimating the data used to train the model, with $MSE_{train}=\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2},$ where $\hat{y}_{i}$ is the expected value of the estimated predictive distribution, and for a test dataset of 100 new data, $MSE_{test}=\sum_{i=1}^{n}\left(y^{test}_{i}-\hat{y}_{i}\right)^{2}$. In addition, we report $LPML_{train}$, value of the Log Pseudo Marginal Likelihood for the training dataset. Table 4 shows the values compared to other competitor models, i.e the linear dependent Dirichlet process mixture (LDDP) defined in De Iorio et al. (2004), the product partition model (PPMx) in Müller et al. (2011) and the linear dependent tailfree process model (LDTFP) in Jara and Hanson (2011). The best values are in bold: our model performs well according to the LPML, while the MSE suggests to use PPMx or LDTFP. In general, our model is competitive with respect to other popular models in the literature. Moreover, in the LDDP case, we have that the average number of clusters is 20.6 with variance 2.266, thus indicating a less parsimonious model compared to ours.

In summary, our extensive simulations suggest that the proposed approach tends to require less mixture components than other well-known alternative models, while still providing a reasonably good fit to the data.

5 Biopic movies dataset

For this illustrative example we consider the Biopics data available in the R package fivethirtyeight (Ismay and Chunn, 2017). This dataset is based on the IMDB database, related to biographical films released from 1915 through 2014. An interesting explorative analysis of the data can be found in https://fivethirtyeight.com/features/straight-outta-compton-is-the-rare-biopic-not-about-white-dudes/.

We consider the logarithm of the gross earnings at US box office as a response variable, with the following covariates: (i) year of release of the movie (in a suitable scale, continuous); (ii) a binary variable that indicates whether the main character is a person of color; and (iii) a categorical variable that considers if the country of the movie is US, UK or other. After removing the missing data from the dataset, we were left with $n=437$ observations and $p=4$. We note that 76 biopics have a person of color as a subject and the frequencies of the category “origin” are $\left(256,79,64\right)$ for US, UK and “other”, respectively; “other” means mixed productions (e.g. US and Canada, or US and UK). In what follows, the hyperparameters in model (19)-(20), (12)-(13), (17)-(18) are chosen as $\beta_{0}=0$, $(a_{\rho},b_{\rho})=(1,1)$. The prior mean and variance of $K$ induced by these hyperparameters are 2.162 and 1.978, respectively. The scale hyperparameter $\phi$ in the $g$-prior for $\beta$ and $(a_{0},b_{0})$ vary as determined in Table 5, where $m$ and $v$ denote the prior mean ${\displaystyle b_{0}/(a_{0}-1)}$ and variance ${\displaystyle b_{0}^{2}/((a_{0}-1)^{2}(a_{0}-2))}$, respectively, of the inverse gamma distribution for $\sigma_{k}^{2}$ as in (20). We also assume $\gamma_{0}$ equal to the vector of all 0’s, while $\Lambda_{0}$ is such that the marginal a priori variance of $\gamma_{k}$ is equal to $diag(0.01,0.1,0.1,0.1)$, in accordance to the variances of the corresponding frequentist estimators.

The posterior of $K$ is robust with respect to the choice of prior hyperparameters; on the other hand, our results show that by not including covariates in the likelihood, i.e. setting all $\gamma_{k}$’s are equal to 0, inference on $K$ is much more sensitive to the choice of $(a_{0},b_{0})$ (results not shown here).

Predictive inference was also considered, by evaluating the posterior predictive distribution at the following combinations of covariate values: $(i)$ (mean value for covariate year, US, white); $(ii)$ (mean value for covariate year, US, color); $(iii)$ (mean value for covariate year, UK, white); $(iv)$ (mean value for covariate year, UK, color); $(v)$ (mean value for covariate year, “other”, white); and $(vi)$ (mean value for covariate year, “other”, color). Corresponding plots are shown in Figure S.11 in the Supplementary material. These distributions appear to be quite different in the six cases: in particular, we can observe that in cases $(i)$ and $(ii)$, the posterior is shifted towards higher values. This is quite easy to interpret, since the measurements are given by the earnings in the US box offices; therefore, we expect that in general US movies will be more profitable in that market. The difference due to the race is, on the other hand, less evident. However, the predictive densities show slightly higher earnings for movies where the subject is a person of color, if the origin is “other” ($(v)$ and $(vi)$). Movies from the UK, on the other hand, exhibit the opposite behavior ($(iii)$ and $(iv)$).

We report here the posterior cluster estimate for Test B in Table 5. We found three groups, with sizes 10, 193, 234, respectively; see Figure 2 for the estimated clusters and boxplots of the response. As a comparison, it can be useful to report the total average values for the response, 15.36, and for the covariates: 7.89 (year), 0.18 (UK), 0.15 (“other”), 0.83 (white). These 3 groups have a nice interpretation in terms of covariates: group 1 is the smallest, with a high average response (17.18), and it is characterized by a high percentage of movies from “other” countries, with a person of color as its subject. Group 2 corresponds also to a high average response (16.42), but the average values of UK, “other” and person of color are similar to the total averages (0.14, 0.09, 0.84, respectively). The average response in group 3 is smaller (14.40) than the total sample mean, while the average values of UK, “other” and person of color are 0.22, 0.17, 0.84, respectively.

To assess effectiveness of the proposed model, we compare the results with the linear dependent Dirichlet process mixture model introduced in De Iorio et al. (2004) and implemented in the LDDPdensity function of DPpackage (Jara et al., 2011). Prior information has been fixed as follows: for Test G the mass parameter of the Dirichlet process $\alpha$ is set equal to 0.3 such that $\mathbb{E}(K)=2.87$ and $Var(K)=1.81$, that approximately match the prior information we gave on the parameter $K$. Similarly, under Test H, $\alpha$ is distributed according to the $gamma(1/4,1/2)$, such that the prior mean on $K$ is 3.6 and variance 22.18. The normal baseline distribution is a multivariate Gaussian with mean vector 0 and a random covariance matrix which is given a non-informative prior and the inverse-gamma distribution for the variances of the mixture components has parameters such that mean and variance are equal to 5, 1, respectively similarly as in Table 5. Posterior summaries can be found in Table 6.

As a comparison between the estimated partitions under our model (Figure 2) and the LDDP mixture model, Figure S.12 in the Supplementary Material displays the estimated partition obtained under the LDDP model under Test $G$, that has 3 groups with sizes $\{300,127,10\}$.

An additional illustrative application, concerning data on air quality in different locations in North and South America is reported in Section S.5 of the Supplementary material.

6 Conclusion

This work deals with mixture models where the prior has the property of repulsion across location parameters. Specifically, the discussion is centered on mixtures built on determinantal point processes (DPPs), that can be constructed using a general spectral representation. The methods work with any valid spectral density, but for the sake of concreteness, illustrations were discussed in the context of the power exponential case.

Though we limit ourselves to the case of isotropic DPPs, inhomogeneous DPPs can be obtained by transforming or thinning a stationary process. However, we believe that this case is not very interesting, unless there is a strong reason to assume non-homogeneous locations a priori.

Our computational experiments and data illustrations show that the repulsion induced by the DPP priors indeed tends to eliminate the annoying case of very small clusters that commonly arises when using models that do not constrain location/centering parameters. This happens with very small sacrifice of model fit compared to the usual mixture models.

Another advantage of our model over DPMs is that we avoid the delicate choice of the base measure of the Dirichlet process, leading to more robust estimates on the number $K$ of components in the mixture.

Supplementary Material

S.1 Gibbs sampler for model (11)-(16)

Posterior inference for our DPP mixture model as in (11)-(16) is carried out using a Gibbs sampler algorithm. The full-conditionals are outlined below: we provide the details of the computation only when the conditional posterior distribution is not straightforward. In what follows, $rest$ refers to the data and all parameters except for the one to the left of “$\mid$”.

• •

The labels $\{s_{1},\dots,s_{n}\}$ are independently distributed according to a discrete distribution whose support is $\{1,2,\dots,K\}$:

[TABLE]

• •

The distribution of the weights $\{w_{1},\dots,w_{K}\}$ is conjugate: the conditional distribution is still a Dirichlet distribution, where the parameters are $\delta+n_{k}$, $k=1,\dots,K$.

• •

The variances in each component of the mixture $\{\sigma^{2}_{1},\dots,\sigma^{2}_{K}\}$ are generated independently according to the following distribution:

[TABLE]

• •

Sampling the means $\{\mu_{1},\dots,\mu_{K}\}$ needs more care: following the reasoning in Xu et al. (2016), this full-conditional can be written as

[TABLE]

[TABLE]

thanks to the Schur determinant identity. Note that $det[\tilde{C}](\{\tilde{\mu}_{1},\dots,\tilde{\mu}_{K}\};\rho,\nu)$ in the above expression follows from the expression of the density of a DPP on a compact set; see (9). Then, $b$ is a vector defined as $b=\tilde{C}(\tilde{\mu}_{k},\tilde{\mu}_{-k})$, $\tilde{\mu}_{-k}=\{\tilde{\mu}_{j}\}_{j\neq k}$ and $\tilde{C}_{-k}^{-1}$ is a matrix of dimension $(K-1)\times(K-1)$ defined as $\tilde{C}\left(\tilde{\mu}_{-k},\rho,\nu\right)$. Moreover, $\tilde{\mu}_{k}=T(\mu_{k})$ is the transformed variable that takes values on the set $S=[-1/2,1/2]^{d}$. Typically, the rectangle $R$ such that $T(R)=S$ is fixed in such a way that it is large and contains abundantly all the datapoints.

We update each mean $\mu_{k}$ separately for $k=1,\dots,K$ using a Metropolis-Hastings update.

• •

The full-conditional for the parameters $(\rho,\nu)$ is as follows

[TABLE]

The adaptive Metropolis-Hastings algorithm of Roberts and Rosenthal (2009) is employed in this case, in order obtain a better mixing of the chains and to avoid the annoying choice of the parameters for the proposal distribution.

• •

In order to sample $K$ we need a Reversible Jump step: standard proposals to estimate mixtures of densities with a variable number of components are based on moment matching (Richardson and Green, 1997) and have been relatively often used in the literature. The idea is to build a proposal that preserves the first two moments before and after the move, as in Xu et al. (2016). In particular, the only possible moves are the splitting move, passing from $K$ to $K+1$, and the combine move, from $K$ to $K-1$.

• (i)

Choose move type: uniformly choose among split and combine move (however, if $K=1$ the only possibility is to split)

• (ii.a)

Combine: randomly select a pair $(j_{1},j_{2})$ to merge into a new parameter indexed with $j_{1}$. The following relations must hold:

[TABLE]

[TABLE]

[TABLE]

• (ii.b)

Split: randomly select a component $j$ to split into two new components. In this case, we need to impose the following relationships:

[TABLE]

[TABLE]

[TABLE]

where $\alpha\sim Beta(1,1)$, $\beta\sim Beta(1,1)$ and $r\sim Beta(2,2)$.

• (iii)

Probability of acceptance: the proposed move is accepted with probability $\alpha=\min\left(1,\dfrac{1}{q(\mbox{proposed, old})}\right)$ if we selected a combine step, $\min\left(1,q(\mbox{old, proposed)}\right)$ in the split case. In particular,

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

Moreover, $p_{K+1}^{split}=0.5$ if $K>1$, 1 otherwise; $p_{K}^{comb}=0.5$ if $K>1$, 0 otherwise.

S.2 Tests on data from a mixture with 8 components

We simulated a dataset with $n=100$ observations from a mixture of 8 components. Each component is the Gaussian density with mean $\theta_{k}$ and $\sigma^{2}_{k}=\sigma^{2}=0.05$: the means $\{\theta_{k}\}$ are evenly spaced in the interval $(-10,10)$. In the model (11)-(16), we set $a_{0}=2.0025$, $b_{0}=0.050125$ so that $\mathbb{E}(\rho_{0})=0.05$ and $Var(\rho_{0})=1$; again, $s=0.5$ and $\delta=1$.

Table S.7 reports hyperparameters values for different tests and posterior summaries of interest, as well as prior mean and variance of $K$. In particular, we show the posterior mean and variance for the number of components $K$ (with which we assess the effectiveness of the model for clustering), the mean squared error (MSE) and the log-pseudo marginal likelihood (LPML) (that helps quantifying the goodness-of-fit). In all cases we obtained a pretty satisfactory estimate of the exact number of components, which is 8: the posterior is concentrated around the true value with a very small variance. See also Figure S.3.

From the density estimation viewpoint, we have from Table S.7 that both MSE and LPML are similar for all the tests, thus indicating robustness with respect to the prior choice of parameters $\rho$ and $\nu$. However, preferable tests seem to be $S_{2}$ and $S_{7}$; see Figure S.4, where density estimates and estimated partitions for these two cases are displayed. The posterior density of $\rho$ under Tests $S_{2}$ and $S_{7}$ is shown in Figure S.5.

S.3 Different spectral densities: application to the Galaxy dataset

We consider the proposed model with different spectral densities, to check its robustness on the inference. All the models presented in this paper are, as a matter of fact, general and in principle any spectral density $\varphi(\cdot)$ satisfying the conditions for the existence of the DPP process can be employed. The choice of the spectral density in (10) is motivated by its strong repulsiveness (see Lavancier et al., 2015). However, in this section we show inference on the Galaxy dataset obtained when spectral representations other than the power spectral density, drive the DPP.

We choose isotropic covariance functions that are well-known in the spatial statistics literature: the Whittle-Matérn and the generalized Cauchy. Both densities depend on three parameters: intensity $\rho>0$, scale $\alpha>0$ and shape $\nu>0$. In order to assure $\varphi(x)<1$ for all $x$, $\rho$ must be smaller than $\rho_{max}=\alpha^{-d}M$, where $M$ needs to be specified for each of the two cases. For the Whittle-Matérn we have

[TABLE]

and for the generalized Cauchy

[TABLE]

where $d$ is the dimension of the space where $x$ lives ($d=1$ in what follows) and $\mathcal{K}_{\nu}(\cdot)$ is the modified Bessel function of the second kind.

We fix $\rho=\dfrac{1}{2}\rho_{max}$ and $(\alpha,\nu)$ equal to: $(i)$ (0.1,0.1), $(ii)$ (0.1,2), $(iii)$ (1,0.1) in the tests below. To fit the Galaxy data to the model in Section 2.2, the selected hyperparameter values are $\delta=1$ and $(a_{0},b_{0})=(3,3)$ (see (15) and (16)).

Table S.8 displays posterior summaries for the two families of spectral densities under hyperparameter settings $(i)$, $(ii)$ and $(iii)$. Posterior summaries of the number of components $K$ and goodness-of-fit values are close to those of Table 1. This gives evidence to robustness with respect to the choice of the spectral density.

S.4 Gibbs sampler in presence of covariates

The Gibbs sampler algorithm employed to carry out posterior inference for model (19)-(20), (12)-(13), (17)-(18) is different from the one in Section S.1 except for the full conditionals of $(\rho,\nu)$. The sampling of labels $\{s_{i}\}_{i=1}^{n}$ differs from (S.22), since now

[TABLE]

The sampling of the $\{\mu_{k}\}_{k=1}^{K}$ is similar as the same step in Section S.1, but now

[TABLE]

However the substantial change from the model without covariates to the model with covariates is due to the update of $K$, the number of components in the mixture, and of $\{\beta_{2},\dots,\beta_{K}\}$ (recall that $\beta_{1}=0$ for identifiability reasons); these are indeed complicated by the presence of the covariates. Moreover, the update of $\{\sigma_{k}^{2}\}$ is now replaced by

[TABLE]

where $n_{k}=\#\{i:\,s_{i}=k\}$; here we assume the vector of the prior mean of $\gamma_{k}$, $\gamma_{0}$, to be equal to the 0-vector. This full-conditional is the posterior of the standard conjugate normal likelihood, normal - inverse gamma regression model. In particular, we have that

[TABLE]

with $\Lambda^{*}=\left(\Lambda_{0}^{-1}+\sum_{i:\,s_{i}=k}x_{i}x_{i}^{T}\right)^{-1}$ and $m^{*}=\Lambda^{*}\left(\sum_{i:\,s_{i}=k}y_{i}x_{i}\right)$. Moreover

[TABLE]

The full-conditional for the coefficients $\beta_{k}$, $k=2,\dots,K$ is as follows:

[TABLE]

which has no known form. Therefore we resort to a Metropolis Hastings step with a multivariate Gaussian proposal, centered in the current value of the vector and with a diagonal covariance matrix, i.e. $\zeta\mathbb{I}_{p\times p}$, where $\zeta$ is a tuning parameter chosen to guarantee convergence of the chains.

On the other hand, the update of $K$ requires a Reversible Jump-type move. However, the approach used in Section S.1 above is difficult to implement when mixing weights depend on covariates, as in this case, so that we need to find another way to define a transition probability. Our approach is similar to that of Norets (2015), with some differences that will be highlighted in the next lines.

As before, there are two available moves: split or combine. The probability of proposing one of them is 0.5, except if $K=1$, when only the split move can be proposed.

Split: if this move is picked, $K^{prop}=K+1$, so that we need to create a new group and its corresponding parameters (the other parameters are kept fixed):

• (i)

randomly pick one cluster, say $j$, containing at least two items

• (ii)

randomly divide data associated to this group, $\tilde{y}_{j}$, into two subgroups, $\tilde{y}_{j_{1}}$ and $\tilde{y}_{j_{2}}$;

• (iii)

set $\gamma_{j_{1}}=\gamma_{j}$, $\sigma^{2}_{j_{1}}=\sigma^{2}_{j}$, $\beta_{j_{1}}=\beta_{j}$, $\mu_{j_{1}}=\mu_{j}$. Now we need to choose a value for $\gamma_{j_{2}}$, $\sigma^{2}_{j_{2}}$, $\beta_{j_{2}}$ and $\mu_{j_{2}}$. In Norets (2015), this is done by sampling the new values from the posterior, conditioning also on the other parameters (even if, for practical purposes, Gaussian approximations for conditional posteriors are used in the implementation of the algorithm). Instead, we sample $\left(\mu_{j_{2}},\gamma_{j_{2}},\sigma^{2}_{j_{2}}\right)$ from the posterior of the following auxiliary model

[TABLE]

where $\tilde{x}_{j_{2}}$ and $\tilde{y}_{j_{2}}$ represent covariates and responses in the new group with label $j_{2}$, respectively. Parameter $\beta_{j_{2}}$ is sampled from a $p$-dimensional Gaussian distribution with mean $\beta^{mode}$ and variance covariance matrix $\Sigma^{mode}$. In particular, $\beta^{mode}$ is the argmax of the following expression

[TABLE]

which corresponds to an approximation of the full-conditional of the $\beta_{k}$ (we dropped the dependence on the other $\beta_{j}$’s by considering the expected value in the denominator). Note that $\mathbb{E}\left(\exp\left(\beta_{j}^{T}x_{i}\right)\right)$ is not other than the moment generating function, thus it is equal to $\exp\left(\beta_{0}^{T}x_{i}+x_{i}^{T}\Sigma_{0}x_{i}/2\right)$.

Combine: here $K^{prop}=K-1$, so that it suffices to collapse two groups into one. Specifically, we randomly choose one group to delete, say $j_{1}$, and remove the corresponding parameters $\beta_{j_{1}}$, $\mu_{j_{1}}$ and $\sigma^{2}_{j_{1}}$. Then, we choose another group, $j_{2}$, and assign all the data $\tilde{y}_{j_{1}}$ to it.

Acceptance rate: this is simply given by

[TABLE]

[TABLE]

where $\theta_{K}=(\sigma^{2}_{1:K},\gamma_{1:K},\mu_{1:K},\beta_{1:K})$. Moreover, $p_{s}(j)$ is the probability of splitting component $j$ and similarly for the other terms.

S.5 Air quality index dataset

The Air Quality Index (AQI) is an index for reporting air quality, see for instance https://airnow.gov/index.cfm?action=aqibasics.aqi. It describes how clean or polluted the air is, and what associated health effects might be a concern for the population. The Environmental Protection Agency calculates the AQI for five major air pollutants regulated by the Clean Air Act: ground-level ozone, particle pollution, carbon monoxide, sulfur dioxide, and nitrogen dioxide. Data can be obtained from several sources, for instance, from http://aqicn.org/. For a real-time map, see https://clarity.io/airmap/.

For the purpose of this illustration, we investigate the spatial relations in measurements of the AQI made on September 13th, 2015, at 16pm. We consider 1147 locations scattered around North and South America (the values of AQI have been standardized).

We ran the MCMC algorithm to fit model (11)-(13), (16)-(18), with a burn-in of 10,000, a thinning of 10 and a final sample size of 5,000. As before, $\beta_{0}=0$ and $\nu=2$. Table S.9 displays different settings of the hyperparameters for which the prior mean on the number of groups is 1.996 and the prior standard deviation is 1.290 (computed using a Monte Carlo approach). The different hyperparameter settings differ for the specification of $\phi$, the scale hyperparameter in the $g$-prior in (18), and the prior mean $m$ and variance $v$ of $\sigma^{2}_{k}$; see (16).

Figure S.6 shows the estimated clusters obtained under Test $AQ_{1}$. The north - east coast seems to be associated with better environmental conditions, and it is clear that important urban sprawls are generally grouped together. More in detail, the Binder loss function method estimated 6 groups characterized by the following mean and standard deviations of the AQI: (0.95, 0.44) in the red group, (-0.27, 0.45) in the yellow, (-0.70, 0.21) in the green, (1.7,1.64) in the light blue, (0.28, 0.54) in the blue, (-0.51, 0.29) in the pink group; yellow, pink and green points are associated to lower values of AQI, while red and light blue to higher values. The boxplots of the AQI by cluster in Figure S.6 are clearly interpretable: the cluster depicted in light blue gathers the polluted cities in south America and big cities in the West coast of the U.S. (Las Vegas, Los Angeles, Seattle, for instance). On the other hand, yellow and green points indicate less dangerous environmental conditions that characterize the North-East coast: however, the small red cluster contains the big cities that are present in this area (Chicago, New York, Philadelphia, Boston).

Figure S.7 displays three different predictive laws that correspond to different locations: Sacramento, which shows the lowest predicted values of AQI, New York, where the environmental conditions are worse, and Monterrey, that presents an intermediate situation.

Similarly as for the Biopics dataset, we compare the inference under our model with the linear dependent Dirichlet process mixture model introduced in De Iorio et al. (2004). Prior information is fixed as follows: $\alpha$ is distributed according to the $gamma(1,1)$ distribution for Test $AQ_{5}$, so that the prior mean and variance of $K$ are 7.15 and 36, respectively, i.e. the prior of $K$ is vague. On the other hand, in Test $AQ_{6}$ the mass parameter $\alpha$ of the Dirichlet process is set equal to 0.15 such that $\mathbb{E}(K)=2.09$ and $Var(K)=1.02$, which approximately matches the prior information given on $K$ (mean 1.996 and variance 1.29). The baseline distribution is a multivariate Gaussian with mean vector 0 and a random covariance matrix which is given a non-informative prior and the hyperparameters of the inverse-gamma distribution for the variances of the mixture components are such that prior mean and variance are equal to 5 and 1, respectively. Posterior summaries can be found in Table S.10.

S.6 Additional plots

• •

Figure S.8 displays the graph of the power exponential spectral density $\varphi(x;\rho,\nu)$ in (10) when $\rho$ is 2 (left) and 100 (right) and $\nu$ varies.

• •

Figure S.9 displays the value of $C_{0}(t)$ corresponding to the Gaussian spectral density where $s=0.5$ and $\rho$ varies as in the legend. The vertical dashed line represents the upper limit of the set $S=\left(-\dfrac{1}{2},\dfrac{1}{2}\right)$. The approximation $C_{0}(t)\simeq 0$ for $t\notin S$ is perfectly adhered to when $\rho$ is small. The higher the $\rho$ is, the slower is the decay rate of the function $C_{0}(t)$.

• •

Figure S.10 refers to Section 4.2 of the paper, where we considered a simulated dataset with three covariates. Predictive distributions corresponding to the 12 different reference values of the covariates are shown. The simulation truth can be found in Figure 1 of Müller et al. (2011).

• •

Figure S.11 shows the predictive distribution for the response “box-office earnings” in the Biopics dataset application for cases $(i)-(vi)$ with parameter setting E in Table 5.

• •

Figure S.12 diplays the cluster estimate for the Biopics dataset obtained under a linear dependent Dirichlet process model with prior specification $G$ in Table 6 of the paper.