Asymptotically minimax empirical Bayes estimation of a sparse normal mean vector

Introduction

This normal means model is by now a classic one which has been widely studied from both a mathematical and applied point of view. Despite the extent to which the many-normal-means model has been studied, it is still a practically important model in a variety of problems. For example, the sparse normal mean model is the cornerstone for many modern Bayes and empirical Bayes multiple testing procedures, e.g., Scott and Berger (2006), Jin and Cai (2007), Bogdan et al. (2008), Efron (2008), and Martin and Tokdar (2012). More recently, Scott et al. (2013) have presented a novel use of the same classical model considered here but in the regression setting. Clearly, research on this classical model is alive and well, and the results provided by our unique approach, namely, asymptotically minimax concentration rates and superior finite-sample performance compared to many existing methods, are useful contributions.

Recently, Castillo and van der Vaart (2012) have considered the performance of several Bayesian methods for this problem. They focus on frequentist properties of a Bayesian posterior distribution, and the corresponding Bayes estimators, for priors with a two-groups structure. In sparse estimation problems, a two-groups prior puts positive probability on $\theta$ vectors with some exact zero entries, so the marginal prior for each component is a mixture of a continuous distribution and a point-mass at zero. Castillo and van der Vaart (2012) show that, for a suitably chosen two-groups prior, the posterior concentrates around the true signal at the asymptotically optimal minimax rate. From this, concentration properties of posterior quantities, such as the posterior mean, can be derived. An important message in their paper is that care is needed in choosing the prior for the non-zero $\theta$ entries. In particular, they show that priors with too light tails, e.g., Gaussian, give sub-optimal concentration properties. The results presented herein provide similar guidance, though our perspective is quite different.

Here we take a novel empirical Bayes approach. In particular, we present a hierarchical two-groups prior where, given a weight $\omega$, the $\theta_{i}$’s are modeled as independent, with $\theta_{i}=0$ with probability $g_{i}(\omega)$, and $\theta_{i}\sim h_{i}(\theta\mid\omega)$ with probability $1-g_{i}(\omega)$, where the functions $g_{i}$ and $h_{i}$ depend on data $X_{i}$. These functions are defined explicitly in Section 2. To complete the hierarchy, $\omega$ is assigned a prior concentrated near 1. We argue that the effect of the data-dependent prior is mitigated by preventing the posterior from tracking the data too closely. This approach provides some new insights, which we compare with those coming from the fully Bayesian framework of Castillo and van der Vaart (2012).

In Section 3 we present our theoretical framework. First, we show that our empirical Bayes posterior concentrates, with probability 1, around the true mean vector at the optimal minimax rate (with respect to square error loss) for the assumed sparsity class. Concentration rate theorems for empirical Bayes posteriors are relatively scarce in the literature, and our technique for handling the challenges that arise from data appearing in both the likelihood and prior might be useful in other problems; one possible extension is discussed briefly in Section 5. We then show that our empirical Bayes posterior mean is an asymptotically minimax estimator of $\theta$. Finally, we show that, asymptotically, the support of our empirical Bayes posterior has, up to a logarithmic factor, the same effective dimension as the true sparse $\theta$. An interesting observation is that the particular form of the prior on $\omega$ is the main catalyst for concentration of our empirical Bayes posterior.

Section 4 describes computation of our empirical Bayes posterior mean via a straightforward Markov chain Monte Carlo. Simulation results are presented to show that our empirical Bayes posterior mean generally outperforms those Bayesian and non-Bayesian competitors with comparable large-sample properties. In particular, we compare our method with a two hard thresholding estimators (Donoho and Johnstone 1994), Bayes and empirical Bayes estimators based on priors with a two-groups structure (Castillo and van der Vaart 2012; Johnstone and Silverman 2004), and a new estimator based on the one-group Dirichlet–Laplace prior (Bhattacharya et al. 2014). Our proposed empirical Bayes estimator is competitive in all cases considered here, and, in many cases, is strikingly better than the others. Some concluding remarks are given in Section 5.

An empirical Bayes model

For the independent normal mean model, $X_{i}\sim\mathsf{N}(\theta_{i},1)$, $i=1,\ldots,n$, let $p_{\theta_{i}}(x_{i})$ denote the density of $X_{i}$, and, for $x=(x_{1},\ldots,x_{n})$, let $p_{\theta}^{n}(x)=\prod_{i=1}^{n}p_{\theta_{i}}(x_{i})$ denote the corresponding joint density of $X=(X_{1},\ldots,X_{n})$. Define a data-dependent hierarchical prior $\Pi_{X}$ for $\theta=(\theta_{1},\ldots,\theta_{n})$ as follows. Introduce a weight parameter $\omega\in(0,1)$, and take the joint prior distribution for $(\theta_{1},\ldots,\theta_{n},\omega)$, under $\Pi_{X}$, to have density proportional to

where $\alpha>0$, $\kappa\in(0,1)$, and $\sigma^{2}>0$ are parameters to be discussed further in Sections 3–4. A representation of this as a genuine empirical Bayes plug-in prior is given in Section 3.1. The dependence of the prior on $(\alpha,\kappa,\sigma^{2})$ will not be reflected in our notation.

Observe that if $\sigma^{2}<(1-\kappa)^{-1}$, then the prior for $\theta_{i}$ is proper, a mixture of a point mass and a Gaussian centered at $X_{i}$. When $\sigma^{2}>(1-\kappa)^{-1}$, the prior is improper. In any case, the posterior is proper, so this possible impropriety of the prior is not a concern. In fact, $\sigma^{2}=(1-\kappa)^{-1}$ is a critical boundary, corresponding to an improper uniform prior for the non-zero $\theta_{i}$’s; see Section 3.3. The term $\omega^{\alpha n-1}$ in the joint density, which resembles a beta density, turns out to be critical to the success of our proposed method, both in theory and in implementation.

Given data $X=(X_{1},\ldots,X_{n})$ from the normal mean model and the empirical Bayes prior distribution $\Pi_{X}$ for $\theta$, we could combine these to form an empirical Bayes posterior distribution via Bayes theorem. That is, for a suitable set $A$ in the $\theta$-space, define the probability measure

We will investigate concentration properties of the empirical Bayes posterior in Section 3. In particular, we show that the empirical Bayes posterior mean derived from $Q_{n}$ is an asymptotically minimax estimator of $\theta$.

It might seem that our apparent double-use of the data—in the prior and in the likelihood—could lead to a posterior $Q_{n}$ that tracks the data too closely. To see that this is not the case, note that if $|X_{i}|$ is large, then the prior probability for $\theta_{i}=0$, under $\Pi_{X}$, would be rather large. Thus, the prior has an unexpected shrinkage effect, pushing $\theta_{i}$ corresponding to $X_{i}$ with large magnitude towards zero. On the other hand, an $X_{i}$ with large magnitude shifts the prior on the non-zero part further from zero, effectively making the tails heavier, to accommodate large signals. These two phenomena suggest that using data in both the prior and the likelihood will not result in a posterior that tracks data too closely. In fact, our theoretical and numerical results demonstrate that the posterior is doing the right thing, namely, concentrating on the true $\theta$.

Empirical Bayes posterior asymptotics

To start, it will help to look at the proposed model from a different perspective. For mathematical convenience, we shift our focus and rewrite the empirical Bayes posterior $Q_{n}$ using a fractional likelihood. That is, we write $p_{\theta}^{n}(X)=p_{\theta}^{n}(X)^{\kappa}p_{\theta}^{n}(X)^{1-\kappa}$ and move the $1-\kappa$ fraction into the prior $\Pi_{X}$ defined above. The effect of this is an alternative prior for $(\theta,\omega)$ of a very simple form:

To provide some further intuition for the prior (1) presented in Section 2, we may consider a data-free version of the prior in (2), where the $X_{i}$’s are replaced by hyperparameters $\mu_{i}$. The marginal likelihood for $\mu=(\mu_{1},\ldots,\mu_{n})$, given $\omega$, is

and $X_{i}$ is clearly the maximum marginal likelihood estimate of $\mu_{i}$. The use of plug-in estimates for mean hyperparameters was considered in Babenko and Belitser (2010) though in a slightly different context. We get the empirical Bayes prior (1) by plugging in $X_{i}$ for $\mu_{i}$ and undoing the fractional likelihood.

Within this alternative setup, we introduce independent binary latent variables $I_{1},\ldots,I_{n}$, where $I_{i}=1$ if and only if $\theta_{i}=0$. Then, given $\omega$, the indicators $I_{1},\ldots,I_{n}$ are independent $\mathsf{Ber}(\omega)$ variables. These indicators characterize the support of the vector $\theta$; in particular, $\sum_{i=1}^{n}(1-I_{i})$ is the number of non-zero $\theta_{i}$ and is distributed as $\mathsf{Bin}(n,1-\omega)$. The beta prior for $\omega$ is concentrated near 1 for $n$ large, so the support size will tend to be small, consistent with the assumption of sparsity. Castillo and van der Vaart (2012), on the other hand, focus primarily on priors directly on the support size, though this kind of beta–binomial prior is considered in their Example 2.2. We find that direct use of the weight $\omega$ is both theoretically and computationally convenient; see Remark 1.

This version of the empirical Bayes posterior is particular amenable for our asymptotic analysis; see, also Walker and Hjort (2001). The use of pseudo-posteriors, where an inverse temperature parameter plays the role of $\kappa$, has been considered in the statistics and machine learning literature (e.g., Zhang 2006; Jiang and Tanner 2008; Dalalyan and Tsybakov 2008), but our context is different.

2 Lower bound on the denominator

In the normal mean model, let $\theta^{\star}$ denote the true mean vector. Assume that $\theta^{\star}$ is sparse in the sense that most of its entries are zero. To make this more precise, let $\mathcal{S}^{\star}\subset\{1,2,\ldots,n\}$ denote the support of $\theta^{\star}$, i.e., $\theta_{i}^{\star}\neq 0$ if and only if $i\in\mathcal{S}^{\star}$. Let $s_{n}=\#\mathcal{S}^{\star}$ be the cardinality of $\mathcal{S}^{\star}$, and say that $\theta^{\star}$ is $s_{n}$-sparse. Then by sparse we mean that $s_{n}\to\infty$ but $s_{n}=o(n)$ as $n\to\infty$. That is, although $\theta^{\star}$ is $n$-dimensional, its effective dimension is actually much smaller.

Start by rewriting the empirical Bayes posterior $Q_{n}$ once more as

Our overall goal is to show that $Q_{n}$ concentrates its mass near $\theta^{\star}$ with $\mathsf{P}_{\theta^{\star}}$-probability 1. The strategy is to show that the denominator of $Q_{n}$ is not too small, and the numerator, for sets $A_{n}$ away from $\theta^{\star}$, is not too large.

Our first result gives a bound on the denominator of $Q_{n}$, like that which obtains from the familiar Kullback–Leibler property (e.g., Schwartz 1965; Ghosal et al. 1999; Barron et al. 1999; Ghosal et al. 2000; Shen and Wasserman 2001). This lower bound will be used in Section 3.3 to derive vanishing upper bounds on the $Q_{n}$-probability assigned to complements of balls around $\theta^{\star}$. But besides as a tool for proving other things, the following lemma suggests that our empirical Bayes-style prior is sufficiently concentrated around $\theta^{\star}$. As Castillo and van der Vaart (2012) show, without suitable prior concentration, the desired posterior concentration is not possible. Therefore, if we associate lower bounds on the denominator of $Q_{n}$ in (3) with adequate prior concentration, then Lemma 1 says that our prior is sufficiently concentrated around $\theta^{\star}$.

For given $\omega$, the inner expectation involves an average over all configurations of the indicators $(I_{1},\ldots,I_{n})$ defined in Section 3.1. This average is clearly larger than just the case where the indicators exactly match up with the support $\mathcal{S}^{\star}$ of $\theta^{\star}$, times the probability of that configuration. That is,

The term $\omega^{s_{n}-n}(1-\omega)^{s_{n}}$ corresponds to the probability for the configuration of $(I_{1},\ldots,I_{n})$ matching the support $\mathcal{S}^{\star}$. The integral for $i\in\mathcal{S}^{\star}$ is the expectation of the normal density ratio for non-zero $\theta_{i}$ with respect to the $\mathsf{N}(X_{i},\sigma^{2})$ prior. Finally, the product over $i\not\in\mathcal{S}^{\star}$ disappears because $p_{0}(X_{i})=p_{\theta_{i}^{\star}}(X_{i})$ for $i\not\in\mathcal{S}^{\star}$. To further bound this quantity, first pull out the terms $\exp\{\frac{\kappa}{2}(X_{i}-\theta_{i}^{\star})^{2}\}$ in the latter integrand that do not depend on $\theta_{i}$. Since, by the law of large numbers, $s_{n}^{-1}\sum_{i\in\mathcal{S}^{\star}}(X_{i}-\theta_{i}^{\star})^{2}\to 1$, as $n\to\infty$, with $\mathsf{P}_{\theta^{\star}}$-probability 1, this part contributes a factor $\exp\{\frac{\kappa}{2}s_{n}+o(s_{n})\}$ to the lower bound for $D_{n}$. Next,

So, the remaining product over $i\in\mathcal{S}^{\star}$ equals $(1+\kappa\sigma^{2})^{-s_{n}/2}$, and we can conclude that the entire product over $i\in\mathcal{S}^{\star}$ in the lower bound for $D_{n}$ is itself lower bounded by

It remains to bound the first integral over $\omega$. Since $\pi(d\omega)=\alpha n\omega^{\alpha n-1}\,d\omega$, we have

The last inequality follows since $(1-b)^{1-b}>b^{b}$ for small $b>0$. Next, if we write

and use the approximation $-\log(1-x)=x+o(x)$, for $x\approx 0$, then we get a lower bound on the $\omega$-integral of the form:

for $c=\alpha/(1+\alpha)>0$. Putting these pieces together, gives the lower bound

3 Concentration

then we will demonstrate that $Q_{n}(A_{M\varepsilon_{n}})\to 0$ with $\mathsf{P}_{\theta^{\star}}$-probability 1.

The theorem below requires a restriction on $(\kappa,\sigma^{2})$. In particular, we require that, for some $\beta>1$, $(\kappa,\sigma^{2})$ reside in the feasible region

We are particularly interested in large $\beta$, so that $\kappa$ arbitrarily close to 1 can be included. Figure 1 displays a portion of the region $R_{\beta}$, for $\beta=200$. The condition $\sigma^{2}=(1-\kappa)^{-1}$ discussed in Section 2 defines the boundary of $R_{\beta}$, for large $\beta$ and $\kappa\approx 1$.

For any fixed $\beta>1$, take $(\kappa,\sigma^{2})$ in the feasible set $R_{\beta}$. If $\theta^{\star}$ is $s_{n}$-sparse, then there exists $M>0$ such that $Q_{n}(A_{M\varepsilon_{n}})\to 0$ with $\mathsf{P}_{\theta^{\star}}$-probability 1.

Let $N_{n}$ be the numerator for $Q_{n}(A_{M\varepsilon_{n}})$ in (3), i.e.,

Taking expectation of $N_{n}$, with respect to $\mathsf{P}_{\theta^{\star}}$, we get

Write $J_{\omega}(d\theta_{i})$ for the measure defined in the $i$-th product term. Split this into discrete and continuous pieces:

For clarity, we shall work with the discrete and continuous parts separately.

Discrete part. Using the Renyi divergence formula for normal distributions, the discrete term simplifies to $\omega\exp\{-\frac{\kappa(1-\kappa)}{2}(\theta_{i}-\theta_{i}^{\star})^{2}\}\,\delta_{0}(d\theta_{i})$.

Continuous part. An application of Hölder’s inequality, with coefficients $\frac{\beta}{\beta-1}$ and $\beta$, whose reciprocals sum to one, gives

For $(\kappa,\sigma^{2})\in R_{\beta}$, we have $\frac{\kappa\beta}{\beta-1}<1$. Then the same Renyi divergence formula used above gives $\exp\{-\frac{\kappa}{2}\frac{\beta(1-\kappa)-1}{\beta-1}(\theta_{i}-\theta_{i}^{\star})^{2}\}$. The second term in the upper bound equals

After some tedious algebra, this can be rewritten as

Combining the two terms in the upper bound, ignoring the normal density, gives

For $(\kappa,\sigma^{2})$ in the feasible region $R_{\beta}$ in (5), the coefficient on $(\theta_{i}-\theta_{i}^{\star})^{2}$ in the exponential term above is negative.

We can now find a constant $c>0$, depending on $(\kappa,\sigma^{2},\beta)$, such that

Next, take $M$ such that $cM>2$, and then take $K\in(2,cM)$. Then Markov’s inequality gives the upper bound

This upper bound has a finite sum over $n\geq 1$, so the Borel–Cantelli lemma gives that $N_{n}\leq e^{-K\varepsilon_{n}}$, with $\mathsf{P}_{\theta^{\star}}$-probability 1 for all large $n$. Putting together this bound on $N_{n}$ and the one on $D_{n}$ from Lemma 1, we get

Since $s_{n}=o(\varepsilon_{n})$, the exponent diverges to $-\infty$ regardless of the sign on $\eta$. Therefore, $Q_{n}(A_{M\varepsilon_{n}})\to 0$ as $n\to\infty$ with $\mathsf{P}_{\theta^{\star}}$-probability 1.∎

The $\varepsilon_{n}$ concentration rate is driven primarily by the beta prior on the weight $\omega$. In particular, it comes from the term $(s_{n}/n)^{2s_{n}}$ in the lower bound (4) in Lemma 1. This means that the prior for $\theta$, given $\omega$, should be selected so that it does not interfere with the correct rate coming from the lower bound on the denominator of $Q_{n}$.

Castillo and van der Vaart (2012) show that the minimax concentration rate will not hold if the prior on non-zero $\theta$ has too light of tails, e.g., Gaussian. A way to understand this point, from our perspective, is that the Gaussian conditional prior interferes with what the beta prior for the weight $\omega$ is doing. As we have demonstrated, this does not necessarily mean that Gaussian is wrong, but that some adjustments should be made to prevent this interference.

4 Asymptotic minimaxity of the posterior mean

Since the empirical Bayes posterior concentrates around the right place and the right rate, it ought to produce an estimator of $\theta$ with good properties. For this problem, perhaps the most natural choice of estimator is the empirical Bayes posterior mean,

Next we show that $\hat{\theta}_{n}$ is a minimax estimator if $\theta^{\star}$ is $s_{n}$-sparse.

Take $(\kappa,\sigma^{2})$ as in Theorem 1. If $\theta^{\star}$ is $s_{n}$-sparse, then there exists a universal constant $M^{\prime}>0$ such that $\mathsf{E}_{\theta^{\star}}\|\hat{\theta}_{n}-\theta^{\star}\|^{2}\leq M^{\prime}\varepsilon_{n}$ for all large $n$.

For the integration over $A_{M\varepsilon_{n}}$, we again look at the numerator and denominator of $Q_{n}$ separately, as in the previous subsection. The denominator has the same lower bound as in Lemma 1. Take $n$ large enough that, with $\mathsf{P}_{\theta^{\star}}$-probability 1, the lower bound in the lemma holds; then the expectation of the ratio can be bounded by upper bounding the expectation of the numerator, together with the lower bound on the denominator. Expectation of the numerator, with respect to $\mathsf{P}_{\theta^{\star}}^{n}$, proceeds just like in the proof of Theorem 1. This time, we get

But $\|\hat{\theta}_{n}-\theta^{\star}\|^{2}\leq\int\|\theta-\theta^{\star}\|^{2}\,Q_{n}(d\theta)$ by Jensen’s inequality, so $\mathsf{E}_{\theta^{\star}}\|\hat{\theta}-\theta^{\star}\|^{2}\leq M\varepsilon_{n}(1+e^{-\nu\varepsilon_{n}})$. Take $M^{\prime}=2M$ to complete the proof. ∎

5 Effective posterior dimension

where $D_{\theta}=\#\{i:\theta_{i}=0\}$; this fact derives from the full conditionals in Section 4.1 below. So, if $\alpha$ is not too large, and $\omega$ concentrates around $1-s_{n}n^{-1}$, then $D_{\theta}$ concentrates around $n-s_{n}$. Therefore, the posterior distribution for $\theta$ must reside on a space with effective dimension proportional to $s_{n}$.

Let $\delta_{n}=K\varepsilon_{n}n^{-1}$, where $\varepsilon_{n}=s_{n}\log(n/s_{n})$ as before, and $K>0$ is a suitably large constant. Then, under the conditions of Theorem 1,

Write the numerator of $\mathsf{P}(1-\omega>\delta_{n}\mid X)$ as

This is similar to the first display in the proof of Theorem 1. Just as in that proof, we get the following bound on the expectation:

where $c$ is a positive constant and $v=v(\sigma^{2},\beta)$ is a variance term that depends on the particular $\sigma^{2}$ and $\beta$ values. Each integral in the inside product is bounded above by 1, so we get

From Lemma 1, we have that the denominator of $\mathsf{P}(1-\omega>\delta_{n}\mid X)$ is lower bounded by $\exp\{-2\varepsilon_{n}+O(s_{n})\}$ with probability 1 for large $n$. So, for large $n$, we get

If we pick $K$ such that $K\alpha>2$, then the fact that $s_{n}=o(\varepsilon_{n})$ implies that this upper bound approaches zero as $n\to\infty$, proving the claim. ∎

Since the logarithmic term $\log(n/s_{n})$ is small, the practical implication of this result is that the posterior distribution of $\omega$ concentrates around $1-s_{n}n^{-1}$. The simulation results displayed in Figure 3 below confirm this.

Numerical results

Computation of the empirical Bayes posterior mean can be carried out via a simple Gibbs sampler for $\omega$ and $\theta=(\theta_{1},\ldots,\theta_{n})$ based on the full conditionals:

where $D_{\theta}=\#\{i:\theta_{i}=0\}$. That is, first sample from the $\theta\mid\omega$ conditional posterior in (8a), then from the $\omega\mid\theta$ conditional posterior in (8b). Repeat this process to obtain a sample from the full posterior. R code for this Gibbs sampling procedure is available at www.math.uic.edu/~rgmartin. Once the posterior sample is available, the empirical Bayes estimator $\hat{\theta}$, the posterior mean, is obtained by computing a coordinate-wise average of the posterior $\theta$ samples. Besides the posterior mean, many other quantities of interest can be calculated. For example, inclusion probabilities, $\mathsf{P}(\theta_{i}\neq 0\mid X)$, $i=1,\ldots,n$, can be easily calculated. Also, in a function estimation problem, where $\theta_{1},\ldots,\theta_{n}$ are coefficients attached to the fixed basis functions, the posterior samples of the unknown functions are readily available.

Theory and experience suggest that good numerical results are obtained for large $\kappa$ and large $\sigma^{2}$. Throughout, we use $\kappa=0.99$ and $\sigma^{2}=(1-0.99)^{-1}=100$, on the boundary of the feasible region. For $\alpha$, (7) suggests that relatively small values of $\alpha$ are appropriate, so that the $\omega$ posterior can learn from $X$ through $D_{\theta}$. We have found that choosing $\alpha$ to be decreasing with $n$ is a reasonable choice. (This has no consequence on the results in Theorems 1–3.) In particular, in the three examples below, with $n=200,500,1000$ we take $\alpha=0.25,0.10,0.05$, respectively. Alternatively, one could use the data to choose $\alpha$. For example, a method-of-moments estimator of $\alpha$ can be obtained as follows. First, estimate $D=D_{\theta}$ via universal hard thresholding, i.e., $\hat{D}$ equals the number of $X_{i}$ such that $|X_{i}|\leq(2\log n)^{1/2}$. Under the assumed prior, $D$ has a beta–binomial distribution, with expectation $n^{2}\alpha/(n\alpha+1)$. If we set this expectation equal to $\hat{D}$, then solving for $\alpha$ gives a method-of-moments estimator, in particular, $\hat{\alpha}=\hat{D}\{n(n-\hat{D})\}^{-1}$. In our examples below, we use the $n$-dependent but data-free choices of $\alpha$ mentioned above.

2 Simulation studies

For illustration, we first reproduce a simulation study presented in Bhattacharya et al. (2014). In particular, we take samples $X=(X_{1},\ldots,X_{n})$ of dimension $n=200$ from the normal mean model $X_{i}\sim\mathsf{N}(\theta_{i}^{\star},1)$. Recall the sparsity level $s_{n}$ is the number of non-zero $\theta_{i}^{\star}$’s. In this case, we consider $s_{n}=10,20,40$, and the signals are fixed at values $A=7,8$. Table 1 displays estimates of the mean squared error obtained from 100 replications of $X$. In addition to the proposed empirical Bayes posterior mean estimator (EBM), based on $\kappa=0.99$, $\sigma^{2}=100$, and $\alpha=0.25$, the methods being compared are a Dirichlet–Laplace estimator (DL) of Bhattacharya et al. (2014), an empirical Bayes median estimator (EBMed) of Johnstone and Silverman (2004), and a fully Bayes posterior median estimator (PMed1) of Castillo and van der Vaart (2012). A few other methods have been considered in the literature recently, and some comments on why they are omitted from comparison here are given in Remark 3 below. Here, we find that our proposed empirical Bayes estimator is the top performer across all these settings.

Consider a single sample $X$ under the simulation setting described above, with $n=200$, where the first $s_{n}=10$ entries in $\theta^{\star}$ equal $A=7$, and the remaining entries are zero. For the given $X$, the Gibbs sampler is run to obtain a sample from our empirical Bayes posterior distribution of $\theta$. In Figure 2 we plot the posterior inclusion probability $\mathsf{P}(\theta_{i}\neq 0\mid X)$ as a function of the indices $i=1,\ldots,n$. It is evident that the empirical Bayes posterior is able to clearly identify the correct model.

As a second example, we reproduce a simulation study presented in Castillo and van der Vaart (2012). In this case, we look at $n=500$, $s_{n}=25,50,100$, and signals fixed at $A=3,4,5$. Table 2 displays estimates of the mean squared error based on 100 replications. This time, the methods are two fully Bayes posterior mean estimates (PM1 and PM2), two fully Bayes component-wise posterior medians (PMed1 and PMed2), Johnstone and Silverman (2004) empirical Bayes mean (EBM) and median (EBMed), and hard thresholding (HT) and hard thresholding oracle (HTO) rules. Our proposed empirical Bayes estimator, based on $\alpha=0.10$, is competitive when $A=4$, and clearly dominates when $A=5$, just like in the previous illustration. Interestingly, the empirical Bayes estimators are the better performers overall in this case.

One rather unusual observation is that some of the methods have, for given $s_{n}$, a mean square error increasing in the signal size $A$. We find this behavior to be counterintuitive, since it should be easier to detect stronger signals. The two thresholding estimators have decreasing mean square error as $A$ increases, as does our proposed estimator.

To follow up on the mean square error results in Table 2, we also display the posterior distribution of $\omega$ for two separate runs. As indicated from Theorem 3, the posterior distribution for $\omega$ should concentrate around $1-s_{n}n^{-1}$. For both cases in Figure 3, the posterior is concentrated exactly where we expect that it would be.

As a final example, consider a $n=1000$ dimensional mean vector, with the first 10 entries of $\theta^{\star}$ equal 10, the next 90 entries equal $A$, and the remaining 900 entries equal zero. Mean square errors for two Dirichlet–Laplace estimators in Bhattacharya et al. (2014) and our empirical Bayes estimator, based on $\alpha=0.05$, are displayed in Table 3. Here we consider a range of $A$, from $A=2$ to $A=7$. For the smaller signals, $A\leq 4$, the Dirichlet–Laplace estimator, with smaller prior weight $n^{-1}$ is the best, but our estimator is better for larger signals, $A>4$. The larger weight Dirichlet–Laplace prior estimator is dominated by our empirical Bayes estimator.

There are a number of existing methods available for this problem besides those included in our comparisons here. These include the lasso (Tibshirani 1996), the Bayesian lasso (Park and Casella 2008), the horseshoe prior estimator (Carvalho et al. 2010), the empirical Bayes estimators of Jiang and Zhang (2009), Brown and Greenshtein (2009), and, most recently, Koenker and Mizera (2014). Some of these methods, including a version of ours, are compared more extensively in Koenker (2014). Those estimators without minimax guarantees, such as the Koenker–Mizera estimator, can only be motivated by finite-sample simulation studies which, by necessity, are narrowly constructed. On the other hand, our estimator has the desired minimax property and also has the best overall finite-sample performance among those provably minimax competitors.

Discussion

The paper has considered a classical problem of estimating a sparse high-dimensional normal mean vector, and we have proposed a novel empirical Bayes solution. Though the stated prior itself may seem overly informative, we show that the prior induces a sort of shrinkage effect, preventing the posterior from tracking the data too closely. We go on to prove that the empirical Bayes posterior concentrates around $\theta^{\star}$ at the minimax rate, that its mean is an asymptotic minimax estimator, and that its effective dimension agrees with that of the true sparse mean vector.

In addition to the good large-sample properties, our empirical Bayes procedure is easy to compute, and, in a number of cases, the finite-sample performance of our empirical Bayes posterior mean is considerably better than that of existing methods with comparable large-sample properties (Remark 3). Since our method admits a full posterior distribution, any other feature, such as the inclusion probabilities displayed in Figure 2, useful in the signal detection problem, can be readily calculated.

A possible extension of the method presented herein is as follows. Suppose that each $X_{i}$ and $\theta_{i}$ are $r$-vectors, where $r=r_{n}$ possibly depends on $n$. Collecting some of the variables together in vectors introduce a group structure. This structure appears in a variety of applications, and this has motivated developments in model selection and estimation with grouped variables (e.g., Yuan and Lin 2006). Abramovich and Grinshtein (2013) prove asymptotic minimaxity of a Bayes method in this grouped setting, and we expect that similar results can be derived based on the ideas presented here.

Acknowledgements

The authors are thankful to Professor Roger Koenker who gave some helpful comments on an earlier draft as well as suggestions for improving our Gibbs sampler codes.

References