GENERALIZED DOUBLE PARETO SHRINKAGE
Artin Armagan, David B. Dunson and Jaeyong Lee
SAS Institute Inc., Duke University and Seoul National University
Abstract: We propose a generalized double Pareto prior for Bayesian shrinkage estimation and inferences in linear models. The prior can be obtained via a scale mixture of Laplace or normal distributions, forming a bridge between the Laplace and NormalJeffreys’ priors. While it has a spike at zero like the Laplace density, it also has a Student’s like tail behavior. Bayesian computation is straightforward via a simple Gibbs sampling algorithm. We investigate the properties of the maximum a posteriori estimator, as sparse estimation plays an important role in many problems, reveal connections with some wellestablished regularization procedures, and show some asymptotic results. The performance of the prior is tested through simulations and an application.
Key words and phrases: Heavy tails, highdimensional data, LASSO, maximum a posteriori estimation, relevance vector machine, robust prior, shrinkage estimation.
1. Introduction
There has been a great deal of work in shrinkage estimation and simultaneous variable selection in the frequentist framework. The LASSO of Tibshirani (1996) has drawn much attention to the area, particularly after the introduction of LARS (Efron et al. (2004)) due to its superb computational performance. There is a rich literature analyzing the LASSO and related approaches (Fu (1998), Knight and Fu (2000), Fan and Li (2001), Yuan and Lin (2005), Zhao and Yu (2006), Zou (2006), Zou and Li (2008)), with a number of articles considering asymptotic properties.
Bayesian approaches to the same problem became popular with the works of Tipping (2001) and Figueiredo (2003). By expressing Student’s priors for basis coefficients as scale mixtures of normals (West (1987)), and relying on type II maximum likelihood estimation (Berger (1985)), Tipping (2001) developed the relevance vector machine for sparse estimation in kernel regression. In this setting, however, exact sparsity comes with the price of forfeiting propriety of the posterior by driving the scale parameter of the Student’s distribution toward zero. In fact, driving both the scale parameter and the degrees of freedom to zero yields the socalled NormalJeffreys’ prior, . The name emerges due to the fact that the hierarchy follows as , , where the latter is the Jeffreys’ prior on the prior variance of . Figueiredo (2003) proposed an expectationmaximization algorithm for maximum a posteriori estimation under Laplace and NormalJeffreys’ priors, with estimates under the Laplace corresponding to the LASSO. The NormalJeffreys’ prior leads to substantially improved performance with finite samples due to the property of strongly shrinking small coefficients to zero while minimally shrinking large coefficients due to the heavy tails; however, it has no meaning from an inferential aspect as it leads to an improper posterior.
A Bayesian LASSO was proposed by Park and Casella (2008) and Hans (2009). However, these procedures inherit the problem of overshrinking large coefficients due to the relatively light tails of the Laplace prior. StrawdermanBerger priors (Strawderman (1971), Berger (1980)) have some desirable properties yet lack a simple analytic form. Recently proposed priors have been designed to have high density near zero and heavy tails without the impropriety problem of NormalJeffreys. The horseshoe prior of Carvalho, Polson, and Scott (2009, 2010) is induced through a carefullyspecified mixture of normals, leading to such desirable properties as an infinite spike at zero and very heavy tails. They studied sparse shrinkage estimation properties of the horseshoe in a normal means problem. Griffin and Brown (2007, 2010) proposed an alternative class of hierarchical priors for shrinkage with some similarities to the prior we propose, but it lacks a simple analytic form that facilitates the study of some properties.
There is a need for alternative shrinkage priors that lead to sparse point estimates if desired, do not overshrink coefficients that are not close to zero, facilitate straightforward computation even in large cases, and result in a joint posterior distribution that does a good job of quantifying uncertainty. We propose the generalized double Pareto prior which independently finds mention in Cevher (2009). It has a simple analytic form, yields a proper posterior, and possesses such appealing properties as a spike at zero, Student’s like tails, and a simple characterization as a scale mixture of normals that leads to a straightforward Gibbs sampler for posterior inferences. We consider both fully Bayesian and frequentist penalized likelihood approaches based on this prior. We show that the induced penalty in the regularization framework yields a consistent thresholding rule having the continuity property in the orthogonal case, with a simple expectationmaximization algorithm described for sparse estimation in nonorthogonal cases. In another independent work motivated by applications to genome wide associations studies, Lee et al. (2011) consider a generalized prior (McDonald and Newey (1988)) that includes the generalized double Pareto as a special case. Similarities to previous work are limited and our contributions beyond them are (i) the formal introduction of a generalized Pareto density, thresholded and folded at zero, as a shrinkage prior in Bayesian analysis, (ii) the scale mixture representation of the generalized double Pareto in Proposition 1 which is central to our work, (iii) its connection to the Laplace and NormalJeffreys’ priors as limiting cases in Proposition 2, (iv) the resulting fully conditional posteriors in a linear regression setting along with a simple Gibbs sampling procedure, (v) a detailed discussion on the hyperparameters and and their treatment, along with the incorporation of a griddy sampling scheme into the Gibbs sampler, (vi) a detailed analysis of the induced penalty by the generalized double Pareto prior and the properties of the resulting thresholding rule, (vii) an explicit analytic form for the maximum a posteriori estimator in orthogonal cases, (viii) an expectationmaximization procedure to obtain the maximum a posteriori estimate in nonorthogonal cases using the normal mixture representation, (ix) the onestep estimator (Zou and Li (2008)) resulting from the Laplace mixture representation, revealing the connection of the resulting procedure to the adaptive LASSO of Zou (2006), and (x) the oracle properties of the resulting estimators.
2. Generalized Double Pareto Prior
The generalized double Pareto density is
(2.1) 
where is a scale parameter and is a shape parameter. In contrast to (2.1), the generalized Pareto density of Pickands (1975) is parametrized in terms of a location parameter , a scale parameter , and a shape parameter as
(2.2) 
with for and for . The mean and variance for the generalized Pareto distribution are for and for . If we let , (2.2) becomes an exponential density as with mean and variance .
To modify the generalized Pareto density to be a shrinkage prior, we let and reflect the positive part about the origin, assuming , for a density that is symmetric about zero. The mean and variance for the generalized double Pareto distribution are for and for . The dispersion is controlled by and , with controlling the tail heaviness and corresponding to Cauchylike tails and no finite moments.
Figure 2.1 compares the density in (2.1) to Cauchy and Laplace densities for the special case , so that . We refer to this form as the standard double Pareto. Near zero, the standard double Pareto resembles the Laplace density, suggesting similar sparse shrinkage properties of small coefficients in maximum a posteriori estimation. It also has Cauchylike tails, which is appealing in avoiding overshrinkage away from the origin. This is illustrated in Figure 2.1(a). Figure 2.1(b) illustrates how the density in (2.1) changes for different values of and .
Prior (2.1) can be represented as a scale mixture of normal distributions leading to computational simplifications. As shorthand notation, let denote that has density (2.1).
Proposition 1.
Let , , and , where and . The resulting marginal density for is .
Proposition 1 reveals a relationship between the prior in (2.1) and the prior of Griffin and Brown (2007), with the difference being that Griffin and Brown (2007) place a mixing distribution on leading to a marginal density on with no simple analytic form.
In Proposition 2 we show that the prior in (2.1) forms a bridge between two limiting cases – Laplace and NormalJeffreys’ priors.
Proposition 2.
Given the representation in Proposition 1, implies

for and ,

for , and .
Proof.
For the first item, setting implies placing a Jeffreys’ prior on , . Integration over yields , which implies the NormalJeffreys’ prior on . For the second item, notice that , where denotes the Dirac delta function, since and . Thus, . ∎
As noted in Polson and Scott (2010), if has exponential or lighter tails, observations are shrunk towards zero by some nondiminishing amount, regardless of size. This phenomenon is wellunderstood and commonly observed in estimation under the Laplace prior, where an exponential density mixes a normal density. The higherlevel mixing (over ) in Proposition 1 allows to have heavier tails, remedying the unwanted bias.
As grows, the density becomes lighter tailed, more peaked and the variance becomes smaller, while as grows, the density becomes flatter and the variance increases. Hence if we increase , we may cause unwanted bias for large signals, though causing stronger shrinkage for noiselike signals; if we increase we may lose the ability to shrink noiselike signals, as the density is not as pronounced around zero; and finally, if we increase and at the same rate, the variance remains constant but the tails become lighter, converging to a Laplace density in the limit. This leads to overshrinking of coefficients that are away from zero. As a typical default specification for the hyperparameters, one can take . This choice leads to Cauchylike tail behavior, which is wellknown to have desirable Bayesian robustness properties.
To motivate this default choice, we assess the behavior of the prior shrinkage factor , where is the parameter of interest (Carvalho et al. (2010)). As , the prior imposes no shrinkage, while as it has a strong pull towards zero. The generalized double Pareto distribution implies a prior on upon integration over in Proposition 1. For the standard double Pareto, this is
where denotes the complementary error function. In Figure 2.2, we compare under the standard double Pareto, StrawdermanBerger, horseshoe, and Cauchy priors, which may all be considered default choices. The priors behave similarly for , implying similar tail behavior. The behavior of for governs the strength of shrinkage of small signals. As , tends towards zero for the Cauchy, implying weak shrinkage, while is unbounded for the horseshoe, suggesting a strong pull towards zero for small signals. The StrawdermanBerger and standard double Pareto priors are a compromise between these extremes, with bounded for in both cases. The standard double Pareto assigns higher density close to one than the StrawdermanBerger prior, and has the advantage of a simple analytic form over the StrawdermanBerger and horseshoe priors.
Of course it is best to adjust and according to any available prior information pertaining to the sparsity structure of the estimated vector. For general and values, the prior on is
(2.3)  
where denotes the confluent hypergeometric function. Note that takes a “horseshoe” shape when . Carvalho, Polson, and Scott (2010) show that implies a NormalJeffreys’ prior on , which can also be observed by setting in (2.3) in conjunction with Proposition 1. Hence is unbounded at forcing to be unbounded at only if . The effects of and are now observed with better clarity from Figure 2.3. As increases, less and less density is assigned to the neighborhood of , repressing shrinkage. On the other hand, increasing values place more and more density in the neighborhood of promoting further shrinkage. This notion is later reinforced by Proposition 3, such that the prior induces a thresholding rule under maximum a posteriori estimation if . Hence, we need to carefully pick these hyperparameters, in particular , as there is a tradeoff between the magnitude of shrinkage and tail robustness.
3. Bayesian Inference in Linear Models
Consider the linear regression model , where is an dimensional vector of responses, is the design matrix and . Letting independently for ,
(3.1) 
From Proposition 1, this prior is equivalent to , with and . We place the Jeffreys’ prior on the error variance, .
Using the scale mixture of normals representation, we obtain a simple data augmentation Gibbs sampler having the conditional posteriors , , , , where and InvGauss denotes the inverse Gaussian distribution with location and scale parameters and . In our experience, this Gibbs sampler is efficient with fast rates of convergence and mixing.
In the absence of any prior information on and , one may either set them to their default values or, as an alternative, choose hyperpriors to allow the data to inform about the values of and . We use and to correspond to generalized Pareto hyperpriors with location parameter , scale parameter and shape parameter . The median value of the resulting distribution for and is , centered at the default choices suggested earlier, while the mean and variance do not exist.
For sampling purposes, let and . These transformations suggest a uniform prior on and in given the generalized Pareto priors on and . Consequently, the conditional posteriors for and are
We propose the embedded griddy Gibbs (Ritter and Tanner (1992)) sampling scheme:

Form a grid of points in the interval .

Calculate .

Normalize the weights, .

Draw a sample from the set with probabilities , and set to be used at the current iteration of the Gibbs sampler.
Repeat the same procedure for and obtain a random draw for . We also experiment with fixing as while treating as unknown. In this case, the prior variance of is determined by .
In what follows we establish the ties between the Bayesian approach we have taken and some frequentist regularization approaches. The simple analytic structure of the generalized double Pareto prior facilitates analyses while its hierarchical formulation leads to straightforward computation.
4. Sparse Maximum a Posteriori Estimation
The generalized double Pareto distribution can be used not only as a prior in a Bayesian analysis, but also to induce a sparsityfavoring penalty in regularized least squares:
(4.1) 
where is initially assumed to have orthonormal columns and denotes the penalty function implied by the prior on the regression coefficients. Following Fan and Li (2001), let , and denote the minimization problem in (4.1) for a component of as
(4.2) 
with the penalty function that simply retains the term in that depends on .
From Fan and Li (2001), a good penalty function should result in an estimator that is (i) nearly unbiased when the true unknown parameter is large, (ii) a thresholding rule that automatically sets small estimated coefficients to zero to reduce model complexity, and (iii) continuous in data () to avoid instability in model prediction. In the following, we show that the penalty function induced by prior (3.1) may achieve these properties.
4.1. Nearunbiasedness
The first order derivative of (4.2) with respect to is , where is the term causing bias in estimation. Although it is appealing to introduce bias in small coefficients to reduce the mean squared error and model complexity, it is also desirable to limit the shrinkage of large coefficients with as . In addition, it is desirable for to approach zero rapidly, implying shrinkage, and the associated introduction of bias rapidly decreases as coefficients get further away from zero. In fact, the rate of convergence of to zero is of the same order under the generalized double Pareto and NormalJeffreys’ priors, with . As controls the tail heaviness in the generalized double Pareto prior, with lighter tails for larger values of , convergence of the ratio to is intuitive. In the case of LASSO, the bias, , remains constant regardless of , which can also be observed in Figure 4.4(b).
4.2. Sparsity
As noted in Fan and Li (2001), a sufficient condition for the resulting estimator to be a thresholding rule is that the minimum of the function is positive.
Proposition 3.
Under the formulation in Proposition 1, prior (3.1) implies a penalty yielding an estimator that is a thresholding rule if .
This result is obtained by finding the minimum of and taking it greater than zero. The thresholding is a direct consequence of the fact that when , which requires that , the derivative of (4.2) is positive for all positive and negative for all negative . In this case, the penalized least squares estimator is zero. When , two roots may exist. The larger one (in absolute value) or zero is the penalized least squares estimator. To elaborate more on this, the root(s) may exist for only when . A helpful illustration is Figure 3 of Fan and Li (2001).
4.3. Continuity
Continuity in data is important if an estimator is to avoid instabilities in prediction. As in Breiman (1996), “a regularization procedure is unstable if a small change in data can make large changes in the regularized estimator”. Discontinuities in the thresholding rule may result in inclusion or dismissal of a signal with minor changes in the data used (see Figure 4.4(b)). Hardthresholding, the “usual” variable selection, is an unstable procedure, while ridge and LASSO estimates are considered stable.
A necessary and sufficient condition for continuity is that the minimum of the function is at zero (Fan and Li (2001)). For our prior, the minimum of this function is obtained at . Therefore yields an estimator with this property.
Proposition 4.
Under the formulation in Proposition 1, a subfamily of prior (2.1) with yields an estimator with the continuity property.
In this particular case, the penalized likelihood estimator is set to zero if . When ,
(4.3) 
As can be observed in Figure 4.4(a), ensuring continuity by letting creates a tradeoff between sparsity and tailrobustness. As the thresholding region becomes wider, the larger values are penalized further, yet not nearly at the level of LASSO.
4.4. Maximum a Posteriori Estimation via ExpectationMaximization
We assume a normal likelihood to formulate the procedure for nonorthogonal linear regression. Estimation is carried out via the expectationmaximization (EM) algorithm.
4.4.1. Exploiting the Normal Mixture Representation
We take the expectation of the logposterior with respect to the conditional posterior distributions of and at the th step, then maximize with respect to and to get the values for the step.

Estep:

Mstep: Letting , we have
We refer to this estimator as GDP(MAP).
4.4.2. Exploiting the Laplace Mixture Representation and the Onestep Estimator
In the proof of Proposition 1, the integration over leads to a Laplace mixture representation of the prior. Since the mixing distribution of the Laplace is a known distribution the required expectation is obtained with ease, resulting in the maximization step,
(4.4)  
where , , and . The componentspecific multiplier on is obtained from the expectation of with respect to its conditional posterior distribution, . Similar results to (Ch4.Ex11) are in Candes, Wakin and Boyd (2008), Cevher (2009), and Garrigues (2009).
An intuitive relationship to the adaptive LASSO of Zou (2006) and the onestep sparse estimator of Zou and Li (2008) can be seen via the Laplace mixture representation. As a computationally fast alternative to estimating the exact mode via the above EM algorithm, we can obtain a “onestep estimator” and exploit the LARS algorithm as in Zou and Li (2008). The onestep estimator is
(4.5) 
with and . This estimator resembles the adaptive LASSO. The LARS algorithm can be used to obtain very quickly. We refer to this estimator as GDP(OS).
Remark 1. For , the GDP(OS) solution path for varying is identical to the adaptive LASSO solution path with (see (4) in Zou (2006)) using identical .
Remark 2. GDP(OS) forms a bridge between the LASSO and the adaptive LASSO: as and , GDP(OS) gives the LASSO solution with penalty parameter .
We derive the GDP(OS) estimator only to reveal a close connection with the adaptive LASSO of Zou (2006) and do not use it in our experiments.
4.4.3. Normal vs. Laplace Representations in Computation
As pointed out by an anonymous referee, it is appropriate to compare the convergence behavior of the EM algorithms that exploit different mixture representations. We generated observations from , where the were independent standard normals for , , and . We set the first components of to be and the rest to . For each combination we simulated data sets and ran the EM algorithms obtained from normal and Laplace scale mixture representations. Figure 4.5 illustrates the number of iterations taken by the two algorithms until . As expected, the convergence under the Laplace mixture representation was much faster with the intermediary mixing parameter integrated out rather than using the expectation step in the EM algorithm.
4.5. Oracle Properties
Following Zou (2006) and Zou and Li (2008), we show that the GDP(MAP) and GDP(OS) estimators possess oracle properties. Relaxing the normality assumption on the error term leads to two conditions for Theorem 2 and Theorem 3.

where are independent and identically distributed with mean and variance .

, where is a positive definite matrix.
In what follows, , retains the entries of indexed by , and retains the rows and columns of indexed by .
Theorem 1.
Let
denote the GDP(MAP) estimator, where and . Let . Suppose that , and, . Then is

consistent in variable selection in that ;

asymptotically normal with .
Remark 3. More generally, the above results hold if and .
Theorem 2.
Let denote the GDP(OS) estimator in (4.5) and . Suppose that , , and . Then is

consistent in variable selection in that ;

asymptotically normal with .
The proofs are deferred to Section 8.
5. Experiments
5.1. Simulation
In this section, we compare the proposed estimators to the posterior means obtained under the normal, Laplace, and horseshoe priors, to the Bayesian model averaged (BMA) estimator, as well as to the sparse estimates resulting from LASSO (Tibshirani (1996)) and SCAD (Fan and Li (2001)). GDP(PM) and GDP(MAP) denote the posterior mean and the MAP estimates, respectively, under the generalized double Pareto prior. Hyperparameter values are provided in footnotes of Tables 5.1 and 5.2 when fixed in advance and are otherwise treated as random with the priors specified in Section 3. When not fixed, we first obtain the posterior means of the hyperparameters from an initial Bayesian analysis, then use them in the calculation of the MAP estimates.
We generated observations from , where the were standard normals with , , and . We used the following configurations:
Model 1: randomly chosen components of set to and the rest to .
Model 2: randomly chosen components of set to and the rest to .
Model 3: randomly chosen components of set to and the rest to .
Model 4: randomly chosen components of set to and the rest to .
Model 5: .
In our experiments and the columns of were centered and the columns of scaled to have unit length. For the calculation of competing estimators we used lars (Hastie and Efron (2011)), SIS (Fan et al. (2010)), monomvn (Gramacy (2010)) and BAS (Clyde and Littman (2005), Clyde, Ghosh, and Littman (2010)) packages in R. We mainly followed the default settings provided by the packages. Under the normal prior, the socalled “ridge” parameter was given an inverse gamma prior with shape and scale parameters . Under the Laplace prior, as a default choice, a gamma prior was placed on the “LASSO parameter” , as given in (6) of Park and Casella (2008), with shape and rate parameters and , respectively. Under the horseshoe prior, the monomvn package uses the hierarchy given in Section 1.1 of Carvalho, Polson, and Scott (2010). For BMA, we used the default settings of the BAS package that employs the ZellnerSiow prior given in Section 31 of Liang et al. (2008). The tuning for LASSO and SCAD were carried out by the criteria given in Yuan and Lin (2005) and Wang, Li, and Tsai (2007), respectively, avoiding crossvalidation.
data sets were generated for each case. In Table 5.1, we report the median model error. Model error was calculated as , where is the variancecovariance matrix that generated and denotes the estimator in use. The values in the subscripts give the bootstrap standard error of the median model error values obtained. The bootstrap standard error was calculated by generating 500 bootstrap samples from 100 model error values, finding the median model error for each case, and then calculating the standard error for it. Under each model, the best three performances are boldfaced in the tables.
GDP(PM) estimates showed a similar performance to that of horseshoe under sparse setups. GDP(PM) (with and unknown) also showed great flexibility in adapting to dense models with small signals. GDP(MAP) estimates performed similarly to SCAD and much better than LASSO, particularly so with increasing sparsity, signal and/or sample size. The GDP(PM) and GDP(MAP) calculations are straightforward and computationally inexpensive due to the normal (and Laplace) scale mixture representation used. Being able to use a simple Gibbs sampler (especially when ) makes the procedure attractive for the average user.
Letting may be somewhat restrictive if the underlying model is very dense or very sparse, but in the cases we considered, it performed comparably to others and we believe that it constitutes a good default prior similar to standard Cauchy with the added advantage of thresholding ability. Although we do not take up cases in this paper, in such situations much larger values of would need to be chosen to adjust for multiplicity.
5.2. Inferences on Hyperparameters
Here we take a closer look at the inferences on the hyperparameters obtained from an individual data set for Models 2 and 5 from Section 5.1. This gives us some insight into how and are inferred with changing sample size and sparsity structure. Note that is more restrictive than GDP(PM) as is fixed, treating only as unknown. Figure 5.6 gives the marginal posteriors of and in cases of and GDP(PM) as described in Section 5.1, while Table 5.2 reports the posterior means for and , as well as model error (ME) performance (as calculated in Section 5.1) on the particular data set used. We clearly observe the adaptive nature and higher flexibility of GDP(PM) moving from a sparse to a dense model with a big increase, particularly in , flattening the prior on . There is not quite as much wiggle room in the case of . All it can do is to drive smaller to allow heavier tails to accommodate a dense structure. As observed in Table 5.1, however, performs comparably in sparse cases.
6. Data Example
We consider the ozone data analyzed by Breiman and Friedman (1985) and by Casella and Moreno (2006). The original data set contains 13 variables and 366 observations. The modeled response is the daily maximum onehour averaged ozone reading in Los Angeles over 330 days in 1976. There are predictors considered and deleting incomplete observations leaves observations. For validation, the data were split into a training set containing 180 observations and a test set containing 23 observations. We considered models including main effects, quadratic, and twoway interaction terms resulting in possible subsets. The complex correlation structure of the data is illustrated in Figure 6.7.
Figure 6.8 summarizes the performance of the proposed estimators and their competitors. Median values for and the standard error intervals were obtained by running the methods on different random trainingtest splits. Standard errors were computed via bootstrapping the medians times.
The median number of predictors retained in the model by all three GDP(MAP) estimates was only while it was and for LASSO and SCAD. Hence GDP(MAP) promoted much sparser models. In terms of prediction, yielded the second best results after BMA, with , GDP(PM), and the horseshoe estimator all having somewhat worse performance. These shrinkage priors are designed to mimic model averaging behavior, so we expected to obtain results that were competitive with, but not better than, BMA. The improved performance for may be attributed to the use of default hyperparameter values that were fixed in advance at values thought to produce good performance in sparse settings. Treating the hyperparameters as unknown is appealing from the standpoint of flexibility, but in practice the data may not inform sufficiently about their values to outperform a good default choice. and SCAD both performed within the standard error range of LASSO, while retaining a smaller number of variables in the model. As it is important to account for model uncertainty in prediction, the posterior mean estimator under the GDP prior is appealing in mimicking BMA. In addition, obtaining a simple model containing a relatively small number of predictors is often important, since such models are more likely to be used in fields in which predictive black boxes are not acceptable and practitioners desire interpretable predictive models.
7. Discussion
We have proposed a hierarchical prior obtained through a particular scale mixture of normals where the resulting marginal prior has a folded generalized Pareto density thresholded at zero. This prior combines the best of both worlds in that fully Bayes inferences are feasible through its hierarchical representation, providing a measure of uncertainty in estimation, while the resulting marginal prior on the regression coefficients induces a penalty function that allows for the analysis of frequentist properties under maximum a posteriori estimation. The resulting posterior mean estimator can be argued to be mimicking a Bayesian model averaging behavior through mixing over higher level hyperparameters. Although Bayesian model averaging is appealing, it can be argued that allowing parameters to be arbitrarily close to zero instead of exactly equal to zero may be more natural in some problems. Hence we have a procedure that not only bridges two paradigms – Bayesian shrinkage estimation and regularization – but also yields three useful tools: a sparse estimator with good frequentist properties through maximum a posteriori estimation, a posterior mean estimator that mimics a model averaging behavior, and a useful measure of uncertainty around the observed estimates. In addition, the proposed methods have substantial computational advantages in relying on simple blockupdated Gibbs sampling, while BMA requires sampling from a model space with models. Given the simple and fast computation and the excellent performance in small sample simulation studies, the generalized double Pareto should be useful as a shrinkage prior in a broad variety of Bayesian hierarchical models, while also suggesting close relationships with frequentist penalized likelihood approaches. The proposed prior can be used in generalized linear models, shrinkage of basis coefficients in nonparametric regression, and in such settings as factor analysis and nonparametric Bayes modeling.
8. Technical Details
Proof of Theorem 1.
The proof follows along similar lines as does the proof of Theorem 2 in Zou (2006). We first prove asymptotic normality. Let and
Let , suggesting . Now
and we know that and . If , then . Consider the limiting behavior of the third term, noting that