Robust Bayesian nonparametric variable selection for linear regression

Spike‐and‐slab and horseshoe regressions are arguably the most popular Bayesian variable selection approaches for linear regression models. However, their performance can deteriorate if outliers and heteroskedasticity are present in the data, which are common features in many real‐world statistics and machine learning applications. This work proposes a Bayesian nonparametric approach to linear regression that performs variable selection while accounting for outliers and heteroskedasticity. Our proposed model is an instance of a Dirichlet process scale mixture model with the advantage that we can derive the full conditional distributions of all parameters in closed‐form, hence producing an efficient Gibbs sampler for posterior inference. Moreover, we present how to extend the model to account for heavy‐tailed response variables. The model's performance is tested against competing algorithms on synthetic and real‐world datasets.


Introduction
Bayesian variable selection is a popular tool in statistics and machine learning that can be used for feature selection in linear regression models.The two most popular models are arguably the spike-and-slab and horseshoe regression models.Both models rely on choosing a sparsity inducing prior over the regression coefficients β β β ∈ R p .In the spike-and-slab model, the prior is a mixture between a point mass at zero and a diffuse prior, while in the horseshoe model, a continuous prior balances the local and global shrinkage (see Section 2).Compared to the estimator for the regression coefficients in standard linear regression, the posterior distribution under these priors produces a shrinkage effect toward zero in the posterior point estimates.This is especially useful in the p >> n regime, in which the number of attributes p can be much larger than the number of observations n.
In this paper we propose extensions of the spike-and-slab and horseshoe regression models to deal with heteroskedasticity and outliers in the data.Specifically, we propose a Bayesian nonparametric model in which each observation is endowed with its own specific variance σ 2 i , which is sampled from an unknown distribution P .The distribution P is a Dirichlet process prior and the resulting model is a Dirichlet process scale mixture model.Due to the discreteness of the Dirichlet process, the resulting vector of variances (σ 2 1 , . . ., σ 2 n ) will be partitioned into groups, hence also producing a corresponding clustering of the observations, in which observations belonging to the same group will have the same conditional variance.Moreover, some observations may be allocated to a single group having much larger variance than the others, hence allowing for outliers in the data.
The key features of our proposed model are: • Parsimonious regression construction: Our nonparametric method performs feature selection with the posterior concentrating on the most relevant prediction coefficients.However, unlike variable coefficient models, there is no increase in the degrees-of-freedom associated with the number of regression parameters (or smoothness of parameters to control, c.f. functional regression).
• Interpretable model structure: The proposed linear regression model changes the structure of the marginal variance components while retaining the highly interpretable properties of a standard sparse regression formulation.• Posterior credible intervals: A fully Bayesian approach gives credible intervals for the regression coefficients.We propose several ways to perform posterior inference selection using these intervals, which provide uncertainty quantification for decision makers working with high-dimensional data.• Efficient inference: Under our model, full conditional distributions of all parameters can be derived in closed form, hence producing an efficient Gibbs algorithm that gives a tuning-free and rejection-free Markov chain Monte Carlo (MCMC) sampler.
The rest of the paper is organized as follows.Section 2 briefly reviews the spike-andslab and horseshoe regression models.Section 3 introduces the Bayesian nonparametric framework and our Dirichlet process mixture model for linear regression along with details of our MCMC sampler.In Subsection 3.3, we provide an extension to our model to account for heavy-tailed response variables.Finally, Section 4 presents an empirical comparison of our proposed model against popular alternative models on synthetic data and real-world data examples.

Bayesian Variable Selection for Linear Regression
The Bayesian approach to variable selection for linear regression models is to introduce sparsity-inducing priors on the regression coefficients.The two most popular approaches in the literature, which we shall review here, are the discrete mixture priors known as the spike-and-slab (Mitchell and Beauchamp, 1988;George and McCulloch, 1993), and the continuous shrinkage priors, most notably the horseshoe prior (Carvalho et al., 2010).

Spike-and-slab priors
The original spike-and-slab model was initially proposed by Mitchell & Beauchamp (1988) and significantly developed by Madigan & Raftery (1994) and George & McCulloch (1997).The final adjustments to the model were completed by Ishwaran & Rao (2005) in what they refer to as the stochastic variable selection model.The spike-andslab prior is intuitively simple and consists of two components.The spike is a delta function centered at zero indicating β j ≈ 0, and the slab gives probability mass to non-zero coefficients.
In our regression framework, we have data y ∈ R n that is explained by a matrix of attributes X ∈ R n×p and coefficients β ∈ R p .Assuming a spike-and-slab prior for β, we have the following model, The likelihood for the data assumes a standard Gaussian regression model, and the prior for β is a scale mixture of Gaussians.The variable η j is a latent indicator and δ v0 (•) denotes a point mass at v 0 , where v 0 is a value chosen close to 0, and thus the variance of the prior on β is either almost zero or a broad Gaussian distribution.

Horseshoe priors
Spike-and-slab priors are intuitively appealing, but in practice their discrete nature makes posterior inference computationally difficult.A popular alternative perspective on sparse Bayesian regression is given by the horseshoe prior construction (Carvalho et al., 2009(Carvalho et al., , 2010, see), see).The horseshoe prior is a continuous shrinkage prior which makes posterior computation more efficient when using gradient-based MCMC sampling tools such as STAN (Carpenter et al., 2017).
Under the horseshoe regression model, the parameters λ j and τ are the local and global shrinkage parameters, respectively.Following Carvalho et al. (2010), we choose half-Cauchy priors for λ j and τ .The intuition behind the horseshoe prior is that the global parameter τ will force the regression coefficients towards zero, while the heavy tails of the half-Cauchy prior for the local shrinkage parameters λ j will allow for non-zero β coefficients.
The horseshoe model used in our experiments follows the form originally proposed by Carvalho et al. (2010).A standard Gibbs sampling approach on this parametrization of the model is difficult to implement due to the non-conjugate posterior form of the shrinkage parameters in the model.In our implementation, we introduce auxiliary variables that lead to conjugate full conditionals for all parameters as suggested by Makalic and Schmidt (2015a), and allow for a straightforward implementation of Gibbs sampling.
The parameterization given by Makalic and Schmidt (2015a) makes use of an inverse gamma scale mixture representation of a random variable with a half-Cauchy distribution.It can be shown that Using this decomposition leads to the reparametrized horseshoe model for which the sampling schemes presented in the next section easily follow.
3 Nonparametric Stochastic Variable Selection

Background to Dirichlet processes
In Bayesian nonparametric statistics, unknown parameters are infinite dimensional and cannot be parametrized by a subset of Euclidean space.The most common examples of nonparametric problems are the estimation of density functions, distribution functions and nonparametric regression (Hjort et al., 2010).
The most popular nonparametric distribution function is the Dirichlet process, introduced in Ferguson (1973).In this setting, P has a Dirichlet process distribution with base measure P 0 and concentration parameter α, denoted P ∼ DP (α, P 0 ).For every measurable partition (A 1 , . . ., A k ), the random vector (P (A 1 ), . . ., P (A k )) has a Dirichlet distribution on the k-dimensional simplex with parameters (αP 0 (A 1 ), . . ., αP 0 (A k )).The base measure P 0 is the mean of the prior (i.e.E(P ) = P 0 ), while the concentration parameter α regulates prior uncertainty around P 0 .A large α value implies a strong belief in the prior.

Dirichlet process mixture model
The Dirichlet process mixture model (DPM) (Lo, 1984) assumes each observation y i is sampled from a countable mixture model.The likelihood of each mixture component, K(y i ; θ * k ), is parametrized by an unknown parameter vector θ * k .In the following, K will be a Gaussian kernel and θ * k is the group specific variance.Both the mixture parameters θ * and the mixture weights w k can be encoded into a random measure P = k≥1 w k δ θ * k which is endowed with a Dirichlet process prior, hence producing the following mixture model for i = 1, . . ., n

By introducing latent class variables, {θ
The DPM model has been used in a variety of applications in statistics and machine learning for density estimation and clustering.In density estimation, the goal is to estimate the unknown density of the observations through a mixture model, while in clustering the goal is to cluster observations into groups with a similar distribution.For the latter task, the discreteness property of the Dirichlet process is very convenient, since with positive probability, we will observe ties among the latent variables θ 1:n .Two observations y i and y j having the same value of the latent variables θ i and θ j will be assigned to the same cluster and have the same distribution.
In Section 3.3 we utilize the properties of the Dirichlet process mixture model to propose a Bayesian nonparametric extension to the spike-and-slab and horseshoe regression models.We show that our new model is able to capture heteroskedasticity and outliers in the observations, as well as capturing clustering behaviour in the variance structure.Under both DP variable selection models the mixture component has a likelihood with a cluster-specific variance term which is not fixed apriori, but is learned from the data.

Dirichlet process variable selection
In both the spike-and-slab (2.1) and horsehoe (2.3) models, one assumes that each y i has the same conditional variance σ 2 .Under our nonparametric Dirichlet process model, we introduce an observation dependent variance σ 2 i for each data point.The vector σ 2 := σ 2 1:n is assumed to be sampled from an unknown discrete distribution P sampled from a Dirichlet process.Specifically, we consider the following general hierarchical model, where the prior π for β β β either follows the spike-and-slab construction in lines 2-5 of eq.(2.1), or the horseshoe model in lines 2-4 of eq.(2.3).We refer to the general model in eq.(3.2) as the Dirichlet process variable selection model and recognize that the spike-and-slab and horseshoe models are special cases, which we denote as Dirichlet process spike-and-slab (DPSS) and Dirichlet process horseshoe (DPHS), respectively.
In this model, we are assuming a variance σ i for each observation y i , where the vector of variances σ 1 , . . ., σ n is assumed to be conditionally i.i.d.from an unknown distribution P .Specifically, we use a Dirichlet process as a nonparametric prior for P , centered at IG(b 1 , b 2 ) with concentration parameter α.From the properties of the DP, P will almost surely be a discrete distribution.This implies that, with positive probability, we will observe ties among σ 1 , . . ., σ n .In other words, some σ i will take the same value.We denote by σ * 1 , . . ., σ * Kn the K n distinct values assumed by σ 1 , . . ., σ n .We will then have K n clusters among our observations {y i } n i=1 , where observations within a cluster have the same variance, but different clusters have different variances.Specifically, if σ i = σ j , then y i and y j will be in the same cluster and have the same conditional variance, but with possibly different means, depending on their attributes.

Posterior inference
Under the DPSS and DPHS models, posterior inference can be carried out efficiently using a Markov chain Monte Carlo sampler.We propose a general Gibbs sampler to sample σ 2 1:n (see Algorithm 1) based on the algorithm by Escobar and West (1995), Algorithm 1 DP variable selection Gibbs sampler for t in 1: number of iterations do for i in 1:n do Sample classification value c i with probability (3.3).end for for k in 1: where Sample β, τ j , η j and ω using Algorithm 2 else if Horseshoe then Sample β, λ j and τ using Algorithm 3 end if Sample α, by sampling the latent variables where end for where at each iteration, we first resample a classification vector c 1:n assigning each σ 2 i to a block in the partition, and then resample the distinct values σ 2 * 1:k .The partition generated by the variance value assigned to each observation is distributed, a priori, as a Chinese Restaurant Process (Aldous, 1985).This distribution can be understood intuitively as a Chinese restaurant serving an infinite amount of customers arriving in succession.Each customer needs to be seated on a unbounded number of tables, where the probability of seating a customer at an occupied table is proportional to the number of people already seated at that table.Alternatively, the customer could be seated at a new table with probability proportional to the concentration parameter α.Assuming i.i.d.variances, we can assume, a posteriori, that each observation is the last customer to arrive and each table represents a cluster of variances.The customer/observation will be seated at an already occupied table, or a new table with probability proportional to where K − are the number of distinct c j for j = i, labeled {1, . . ., K − }, and n −i,c is the number of c j = c for j = i.The likelihood of observation i being assigned a new variance, i.e. seated at a new table, is defined as .
Finally, the concentration parameter α is sampled using a data augmentation trick (Ghosal and Van der Vaart, 2017, see pg.89).
We sample the regression coefficients β for the spike-and-slab and horseshoe models using Algorithms 2 and 3, respectively.The Gibbs sampler for the spike-and-slab model follows from Ishwaran and Rao (2005) while the Gibbs sampler for the horseshoe model is based on the data augmentation construction of Makalic and Schmidt (2015b).The full conditional distributions of all parameters in both models can be easily derived and have closed form expressions.Alternatively, particularly for the case of p n, we can get a random sample from the Multivariate Normal distribution of β using the linear solver of Bhattacharya et al. (2016) adapted to our Dirichlet process mixture model.This extension is detailed in Algorithm 4, this algorithm would replace the first step of Algorithms 2 and 3 once the parameters of matrices Σ and Λ are fixed.Since both Σ and Λ are diagonal matrices, the extension changes the cost of sampling β from O(p 3 ) when sampling directly from the multivariate normal to O(n 2 p) when using the linear solver.
An advantage of our Dirichlet process model construction is that the number of clusters K n that partition the vector of variance parameters σ 2 1:n are learned by the DP model and do not need to be fixed apriori.
Algorithm 2 Sample β with the spike-and-slab model Sample where where where 1 , . . ., σ 2 p ) and Λ = diag(τ 2 λ 2 1 , . . ., τ 2 λ 2 p ), for j in 1:p do Sample the local shrinkage parameter λ 2 j and auxiliary variable ν j , ; end for Sample the global shrinkage parameter τ 2 the auxiliary variable ξ, Algorithm 4 Linear solver sampling of β given Σ and Λ. Sample u ∼ N (0, Λ) and δ ∼ N (0, Clustering the variance parameters means that we can account for outlier observations if one of the variances σ * 1 , . . ., σ * Kn is very large.Outlier observations belonging to that specific cluster will have much higher variance than the others and can therefore take more extreme values.In this way, the model can both simultaneously account for heteroskedasticity and outliers which are characterised by clusters with extreme variances.

Heavy Tails: Student-t extension
In many real-world datasets, the response variable y may contain some larger than expected values which are commonly referred to as outliers.A common approach to model such data is to replace the normal distribution assumption for the conditional distribution of y i given x x x i with a Student-t distribution.Compared to a normal distribution, the Student-t distribution has heavier tails and therefore permits outlier observations with a higher probability than under the normal model.We can account for outliers in our Dirichlet process variable selection model (3.2) by replacing the conditional distribution of y i with y i |x i , β, σ 2 ∼ Student-t(x i β, σ 2 i , ν) i = 1, . . ., n, for some degrees of freedom ν.The parameter ν regulates the thickness of the tails of the distribution.A small value for ν leads to thicker tails with the expectation that large or extreme values for y i will be observed.As ν → ∞, the Student-t distribution recovers the normal distribution.
A Gibbs sampler to perform posterior inference for the Student-t model can be derived from Algorithm 1 by including an additional data augmentation step.The required conditional distributions are derived by representing the Student-t distribution as a normal distribution with an inverse gamma variance.The remainder of this section describes these procedures.
We assume our response variables are distributed as Y ∼ Student-t(µ, σ 2 , ν) with density, with mean E(Y ) = µ and variance Var(Y ) = σ 2 ν ν−2 .A well-known reparameterization of the model follows that if, Algorithm 5 DP variable selection Gibbs sampler, Student-t extension for t in 1: number of iterations do for i in 1:n do Sample classification vector c i with probability (3.6).end for for k in 1: end for Sample α, by sampling the latent variables where The extension of the DPSS and DPHS models to a Student-t model is as follows Using the reparameterized Student-t model, and integrating out P from (3.5), the joint where as before, π(σ 1:n |α) denotes the marginal likelihood of σ 1:n when P is integrated out, which can be represented using the Chinese restaurant process.The proposed Gibbs sampler is specified in Algorithm 5, where at each iteration, we first resample a class vector c 1:n assigning each σ 2 i and G i to a block in the partition, and then resample the distinct values σ 2 * 1:k followed by each latent variable G i .In the heavy tailed case, the observation is related to a variance through its latent variable and hence the posterior probability of the classification vector is where K − be the number of distinct c j for j = i, labeled {1, . . ., K − }, and n −i,c is the number of c j = c for j = i.Also, the likelihood of observation i being assigned a new variance is modified to Finally, regression coefficients β are sampled using Algorithms 2 and 3 with small modifications to include the augmented variable G i .

Experiments
In this section, we compare the DPSS and DPHS models against the standard spikeand-slab and horseshoe models of Ishwaran and Rao (2005) and Carvalho et al. (2010).We begin by studying their predictive performance, their ability to model heterogeneous data and their variable selection capabilities over a range of scenarios utilizing synthetic data, then we test the models' performance on a real-world dataset.Specifically, on the reconstruction of a network of genomic data through variable selection on a linear model.
Code for synthetic and real experiments can be found in https://github.com/albcab/RobustVariableSelection.
For all experiments, each algorithm is run for J = 10, 000 iterations with a burn-in period of J/2.Convergence and mixing of the Markov chain is confirmed through visual diagnostics.Hyper-parameters for the models are set to a 1 = b 1 = 2.01, a 2 = b 2 = d 1 = 1 and d 2 = 1/2, which leads to weakly informative priors.This claim is verified with simple cross validation for all the datasets.It is worth noting that one could place additional hyper-priors on these parameters and add extra steps to the MCMC algorithm.However, in doing so we would run the risk of rejecting steps potentially leading to longer mixing times.It is worth noting that in comparable studies (e.g.Ishwaran and Rao, 2005), it is common to simply fix these parameters when assessing performance.Finally, both the dependent and independent variables are normalized to have zero mean for the testing of each model.
In each scenario, and for each combination of attributes and data lengths, 100 synthetic datasets are created and parameters are estimated for each of these datasets.The results presented in the remainder of this section consider all 100 estimates for robustness.

Homoskedastic and Heteroskedastic Errors
We compare the Dirichlet process variable selection models against the standard spikeand-slab and horseshoe models in the presence of homoskedastic (Scenario 1 ) and heteroskedastic (Scenario 2 ) noise. Figure 1 presents boxplots of the posterior error β − β 0 2 / β 0 2 for every combination of number of observations n and parameters p, with Scenario 1 plotted on the left and Scenario 2 given on the right, while Figure 2 presents the same error boxplots when testing the heavy tailed extension of the models.
Robustness -it is clear from the upper right plots on Figure 1 that for Scenario 2 the robustness provided by the DPSS and DPHS models improves the predictive estimates of parameters β when data is sufficient to make the model overdetermined, i.e. when n is at least as large as p.As expected, under Scenario 1 the results are similar for the models allowing heterogeneous and homogeneous variances.Interesting behaviour is observed when the model is underdetermined at varying degrees.When the model is highly underdetermined n p, models with the horseshoe priors provide more robust results in comparison to spike-and-slab priors.However, as n approaches p the spike-and-slab prior provides lower errors than the horseshoe (albeit at a smaller scale), until the model becomes overdetermined and the horseshoe prior models are preferred.Similar results are observed for the heavy tailed extension of the models.In contrast to the Gaussian models, when the models are at least determined n ≥ p the positive effect of robust variance modelling is less noticeable.
Coefficients -looking at Figure 1 and 2 from top to bottom the effect of adding coefficients to the model is illustrated on the predictive capabilities of the models.The Dirichlet process model performs well in the presence of heterogeneity, particularly as the number of model coefficients is less than the number of observations.The robustness added by the heterogeneous variances benefits the horseshoe prior in the overdetermined case and the spike-and-slab prior in the underdetermined case.
Clusters -one of the key advantages of the nonparametric approach we adopt is the availability of posterior densities for the number of variance components K.The Dirichlet process' assumption of a number of clusters which scale logarithmically towards infinity with n can be observed by looking at the boxplots for samples of K for an increasing n in Figures 3 and 4 for the heavy tailed extension.A potentially useful caveat of the heavy tailed extension of the robust models is that the number of clusters K increases at a slower rate than the case of Gaussian errors.In the case where a finite number of clusters in the data is known, or the number of clusters K should scale faster or slower than logarithmically, priors that generalize the Dirichlet process can be used, e.g. the Pitman-Yor process (Pitman and Yor, 1997).It is interesting to notice how, as the number of coefficients increases but the number of observations remain constant, the number of clusters K retract themselves on their bias when the number of coefficients becomes much larger than the number of observations, as shown in Figures 5 and 6.A comparison of these two Figures also shows how the heavy tailed extension of the robust model provides a more gradual increase on the number of clusters K.

Support Recovery
In addition to estimating the regression coefficients, a further task of interest, especially in high-dimensional regression, is to select a subset of attributes which are deemed significant for the predictive model.In a frequentist context, this is usually achieved via forward selection with sequential F-tests, or with 1 or 0 shrinkage and/or selection, for instance through the use of the AIC or BIC criteria.
One advantage of the Bayesian approach is that we can sample from the joint posterior of the coefficients, thus we can construct credible intervals with relative ease.In this section, we compare two methods for estimating the support of the regression model: • The first approach is to simply look at the posterior of the inclusion parameter η j for selecting either the spike, or slab.Specifically, if the mode of the posterior ηj is above level 1 − ζ we add the index j to the estimated support set M. This method works only in the DPSS and SS models.• The second approach is to look at the empirical posterior credible interval of βi at the percentiles ζ/2 × 100 and (1 − ζ/2) × 100.If the constructed interval excludes zero, then we add the index to the estimated support set M. This method is used for the DPHS and HS models.
We note that the second approach is somewhat similar to the z-cut method discussed in Ishwaran and Rao (2005) To summarise the performance of the two approaches, we take the number of true and false positive (TP, FP) coefficients included in the model and consider how this changes as a function of n.Clearly, these values will change as a function of ζ, and in practice this may be tuned to favour either TP or FP results.For simplicity, we report results with ζ = 0.05.Table 1 summarises the performance of the posterior inclusion thresholding method for different values of n in the heteroskedastic scenario (Scenario 2 ) and the standard Gaussian formulation of the models, similar results where observed on their heavy tailed extension.In the case of an overdetermined model, both the standard and robust implementations learn correctly the underlying data generating process.However, when the model is underdetermined the standard spike-and-slab tends to be more conservative in the inclusion of variables while the robust method favour a better estimate TP with the drawback of more FP.Similar results can be observed in the case of the horseshoe prior.The support recovery, or model selection, of the robust models will be detailed further when applied to real world data, where the benefits of each prior are illustrated.

Reconstruction of transcription regulatory networks
We consider the problem of reconstructing genetic regulatory networks from gene expression data (Marbach et al., 2010).This problem can be modelled as a directed network, where each node corresponds to a different gene and each connection represents a directed interaction between two genes at the transcription level.The sparse spike-and-slab and horseshoe priors provide an attractive approach to reconstructing these networks using support recovery methods.
The models are tested on data from the challenge posed at the The Dialogue for Reverse Engineering Assessments and Methods (DREAM) 2009 conference, specifically the multifactorial subchallenge1 .The challenge consists of reverse engineering five independent networks using 100 steady-state measurements for each network with 100 genes.The levels of expression of all the genes are measured under different perturbed conditions.Each gene X i for i = 1, ..., 100 is treated independently and modelled through a linear relationship with the rest of the genes in the network X, i.e.X i = X T β i + .Every β i has an independent spike-and-slab or horseshoe prior and is a normal random variable with mean zero and homogeneous variances in the case of SS and HS, or heterogeneous variances, in the case of DPSS and DPHS.For each gene, the support is recovered using the method described in Section 4.2 for the case of a spike-and-slab prior with β i .In the case of a horseshoe prior, we follow Steinke et al. (2007) and approximate the posterior probability of a connection from gene j to gene i by the probability of the event |(β i ) j | > 0.1 under the posterior for β i .
The performance of the different approaches is evaluated using the mean of the logarithmic loss of the probabilities of connections p ij from gene j to each gene i, defined as 1 100 where y ij = 1 if there is a directed connection from j to i and y ij = 0 otherwise.Table 2 presents these results for each tested model.The improvement in prediction given by heterogeneous variances in the model is substantial.This supports the results from Section 4.2 on synthetic data, illustrating the improvement offered by robust models to correctly classify variables included in the model, specifically in their capability of reducing type II error in the predictions while correctly identifying the significant coefficients in the linear relationship.It is interesting to note that, while the DPSS model provided better support recovery results in the synthetic data example, the DPHS provides better results in the real data example of gene network reconstruction.

Conclusion
This paper outlines a nonparametric extension of popular Bayesian variable selection models to account for heterogeneity, outliers and clustering effects.We also provide an extension to the model to allow for heavy-tailed data distributions.process variable selection models is computationally efficient using a Gibbs sampling construction.The results presented for both synthetic and real data examples show a robust improvement in both predictive accuracy on test data and improved efficiency in identifying key model attributes.

Figure 5 :
Figure 5: Histogram of the derived number of clusters K from the posterior of σ 1 , . . ., σ n for n = 200 and p = 50, 200, and 2000, where the two number of clusters is K = 9.

Figure 6 :
Figure 6: Histogram of the derived number of clusters K from the posterior of σ 1 , . . ., σ n for n = 200 and p = 50, 200, and 2000 on their heavy tailed Student-t extension, where the two number of clusters is K = 9.