Empirical priors and coverage of posterior credible sets in a sparse normal mean model

TL;DR

This paper examines the frequentist validity of Bayesian credible sets in sparse normal mean models, demonstrating that empirical priors can produce credible intervals with accurate coverage under mild conditions.

Contribution

It introduces empirical priors that ensure credible sets attain nominal coverage, even with weaker conditions than traditional selection consistency and Bernstein--von Mises theorems.

Findings

Credible intervals achieve nominal frequentist coverage.

Empirical Bayes method outperforms other fully Bayes methods in coverage.

Conditions for coverage are weaker than those for selection consistency.

Abstract

Bayesian methods provide a natural means for uncertainty quantification, that is, credible sets can be easily obtained from the posterior distribution. But is this uncertainty quantification valid in the sense that the posterior credible sets attain the nominal frequentist coverage probability? This paper investigates the frequentist validity of posterior uncertainty quantification based on a class of empirical priors in the sparse normal mean model. In particular, we show that our marginal posterior credible intervals achieve the nominal frequentist coverage probability under conditions slightly weaker than needed for selection consistency and a Bernstein--von Mises theorem for the full posterior, and numerical investigations suggest that our empirical Bayes method has superior frequentist coverage probability properties compared to other fully Bayes methods.

Figures

Equations

$$Y_i\sim \mathsf{N}({\theta }_i,{\sigma }^2),i=1,\dots ,n,$$

$$S\sim \pi (S)\text{and}({\theta }_S\mid S)\sim {\mathsf{N}}_{|S|}(Y_S,{\sigma }^2{\tau }^{-1}{I}_{|S|}),$$

$$f_n(s)\propto (c{n}^a{)}^{-s},s=0,1,\dots ,n,$$

$$f_n(s)=b_n\left(\genfrac{}{}{0ex}{}{n}{s}\right){\int }_0^1{w}^{n+b_n-s-1}(1-w{)}^{s}dw,$$

$${\Pi }^n(d\theta )\propto L_n(\theta {)}^{\alpha }{\Pi }_n(d\theta ),$$

$${\pi }^n(S)\propto \pi (S)(1+\alpha {\tau }^{-1}{)}^{-|S|/2}{e}^{-\frac{\alpha }{2{\sigma }^2}\Vert Y_{S^c}{\Vert }^2}.$$

$$K_1{n}^{-a_1}\le \frac{f_n(s)}{f_n(s-1)}\le K_2{n}^{-a_2},s=1,\dots ,n,$$

$$\underset{{\theta }^{\star }:|S_{{\theta }^{\star }}|=s^{\star }}{\sup }{\mathsf{E}}_{{\theta }^{\star }}{\Pi }^n(\{\theta \in {\mathbb{R}}^n:\Vert \theta -{\theta }^{\star }{\Vert }^2>M^{\prime }{s}^{\star }\log (n/s^{\star })\})\to 0,n\to \infty .$$

$$\rho_{n}=\Bigl{\{}\frac{2\sigma^{2}(1+\alpha)M\log n}{\alpha}\Bigr{\}}^{1/2}$$

$$\mathsf{E}_{\theta^{\star}}d_{\text{\sc tv}}\bigl{(}\Pi^{n},\mathsf{N}_{|S^{\star}|}(Y_{S^{\star}},v_{\alpha}I_{|S^{\star}|})\otimes\delta_{S^{\star c}}\bigr{)}\to 0$$

$$D_{n}=\sup_{A}\bigl{|}\Pi^{n}(A)-\{\mathsf{N}_{|S^{\star}|}(Y_{S^{\star}},v_{\alpha}I_{|S^{\star}|})\otimes\delta_{S^{\star c}}\}(A)\bigr{|},$$

$${\Pi }^n(A)=\sum _S{\pi }^n(S)\{{\mathsf{N}}_{|S^{\star }|}(Y_{S^{\star }},v_{\alpha }{I}_{|S^{\star }|})\otimes {\delta }_{S^{\star c}}\}(A),$$

$$\sum_{S}\pi^{n}(S)\bigl{|}\{\mathsf{N}_{|S|}(Y_{S},v_{\alpha}I_{|S|})\otimes\delta_{S^{c}}\}(A)-\{\mathsf{N}_{|S^{\star}|}(Y_{S^{\star}},v_{\alpha}I_{|S^{\star}|})\otimes\delta_{S^{\star c}}\}(A)\bigr{|}.$$

$$D_n\le 2\sum _{S\ne S^{\star }}{\pi }^n(S)=2\{1-{\pi }^n(S^{\star })\}.$$

$$H_{n}(t)=\sum_{S}\pi^{n}(S)\,\Pi^{n}(\{\theta:x_{S}^{\top}\theta_{S}\leq t\}\mid S)=\sum_{S}\pi^{n}(S)\,\Phi\Bigl{(}\frac{t-\hat{\psi}_{S}}{v_{\alpha}^{1/2}\|x_{S}\|}\Bigr{)},$$

$$H_{n}(t)=(1-p_{k}^{n})\,1_{[0,\infty)}(t)+p_{k}^{n}\,\Phi\bigl{(}v_{\alpha}^{-1/2}(t-Y_{k})\bigr{)},$$

$$t_{\gamma }=\inf \{t:H_n(t)\ge 1-\gamma \}.$$

$$\mathsf{P}_{\theta^{\star}}\bigl{(}t_{\gamma}\geq x^{\top}\theta^{\star}\bigr{)}\geq 1-\gamma+o(1),\quad n\to\infty.$$

$$\mathsf{E}_{\theta^{\star}}d_{\text{\sc tv}}\bigl{(}\Pi_{\psi}^{n},\mathsf{N}(\hat{\psi}_{S^{\star}},v_{\alpha}\|x_{S^{\star}}\|^{2})\bigr{)}\to 0,$$

$${\mathsf{P}}_{{\theta }^{\star }}(t_{\gamma }\ge 0)$$

$$\ge {\mathsf{P}}_{{\theta }^{\star }}(1-p_k^n\ge 1-\gamma )$$

$$\ge {\gamma }^{-1}\{{\mathsf{E}}_{{\theta }^{\star }}(1-p_k^n)-(1-\gamma )\}$$

$$={\gamma }^{-1}\{\gamma -{\mathsf{E}}_{{\theta }^{\star }}(p_k^n)\}.$$

$$H_n({\theta }_k^{\star })=(1-p_k^n)+p_k^n{V}_k,$$

$${\mathsf{P}}_{{\theta }^{\star }}\{H_n({\theta }_k^{\star })>1-\gamma \}={\mathsf{P}}_{{\theta }^{\star }}\{p_k^n(1-V_k)<\gamma \}.$$

$${\mathsf{P}}_{{\theta }^{\star }}\{p_k^n(1-V_k)\le \gamma \}$$

$$\le {\mathsf{P}}_{{\theta }^{\star }}\{1-V_k<\frac{\gamma }{b}\}+{\mathsf{P}}_{{\theta }^{\star }}\{p_k^n\le b\}$$

$$\le \frac{\gamma }{b}+{\mathsf{P}}_{{\theta }^{\star }}\{p_k^n\le b\}.$$

$${\mathsf{P}}_{{\theta }^{\star }}\{p_k^n\le b\}={\mathsf{P}}_{{\theta }^{\star }}\{1-p_k^n\ge 1-b\}\le (1-b{)}^{-1}{\mathsf{E}}_{{\theta }^{\star }}(1-p_k^n).$$

$$\limsup_{n\to\infty}\,\bigl{[}\text{non-coverage probability}\bigr{]}\leq\tfrac{\gamma}{b}.$$

$$\zeta_{n}^{2}=\tfrac{2\sigma^{2}(1+\alpha)}{\alpha}\bigl{\{}a_{1}\log n+\log\bigl{(}\tfrac{n-s^{\dagger}}{s^{\dagger}+1}\bigr{)}+(s^{\star}-s^{\dagger}-1)\log 2\bigr{\}},$$

$$1-p_k^n={\pi }^n(S\not\ni k)=\sum _{S\in {\mathcal{S}}^0:S\not\ni k}{\pi }^n(S).$$

$${\pi }^n(S)\le \frac{{\pi }^n(S)}{{\pi }^n(S\cup \{k\})}=\frac{\pi (S)}{\pi (S\cup \{k\})}(1+\alpha {\tau }^{-1}{)}^{1/2}{e}^{-\frac{\alpha }{2}|Y_k{|}^2},$$

Full text

Keywords and phrases: Bayesian inference; Bernstein–von Mises theorem; concentration rate; high-dimensional model; uncertainty quantification.

1 Introduction

Consider the high-dimensional normal mean model

[TABLE]

where $Y=(Y_{1},\ldots,Y_{n})^{\top}$ are independent, the variance $\sigma^{2}>0$ is known, and inference on the unknown mean vector $\theta=(\theta_{1},\ldots,\theta_{n})^{\top}$ is desired. We call this “high-dimensional” because there are $n$ unknown parameters but only $n$ data points, with $n\to\infty$. It is virtually hopeless to make quality inference on a high-dimensional parameter without some additional structure to effectively reduce the dimension, and the structure we assume here is sparsity, i.e., most $\theta_{i}$’s are zero. Work on the sparse normal mean problem goes back at least to Donoho and Johnstone, (1994), and it is now a canonical example in the high-dimensional inference literature. Some recent efforts on constructing estimators that achieve the asymptotically optimal minimax rate can be found in Johnstone and Silverman, (2004), Jiang and Zhang, (2009), Castillo and van der Vaart, (2012), Martin and Walker, (2014), van der Pas et al., (2014); van der Pas et al., 2017a , and Ghosh and Chakrabarti, (2015). Extensions to the high-dimensional regression setting can be found in Castillo et al., (2015), Martin et al., (2017), and Ning and Ghosal, (2018).

The first author (RM) learned about and came to appreciate the intricacies of this sparse normal means problem from his primary advisor, Professor Jayanta K. Ghosh, or JKG for short. His book (Ghosh et al., 2006, Chapter 9) gives an excellent introduction to this problem, including the classical James–Stein and full/empirical Bayes solutions, and much of JKG’s later work was devoted to the large-scale significance testing version of this problem, e.g., Bogdan et al., (2008), Bogdan et al., (2011), and Datta and Ghosh, (2013), among others, and there have been many subsequent developments based on his work. Sadly, JKG passed away on September 30th, 2017, before we had a chance to talk with him about our empirical priors. But it is an honor present some new developments on one of JKG’s favorite problems in this paper dedicated to him.

An advantage of a Bayesian approach to these problems is that it yields a posterior distribution which can be used to quantify uncertainty. Computational issues aside, it is easy to construct a posterior credible set for $\theta$. The question is: does the credibility level of the aforementioned set equal the frequentist coverage probability of the set? If so, then we say that the posterior provides valid uncertainty quantification. It is not clear, however, whether one should expect this property to hold in problems where the dimension is increasing with the sample size. In fact, there are negative results (e.g., Li, 1989) which say, roughly, credible sets based on posterior distributions that adaptively achieve the optimal concentration rate cannot simultaneously provide valid uncertainty quantification uniformly over the parameter space. That is, for any posterior with the optimal concentration rate, there must be true parameter values that the $100(1-\gamma)$% posterior credible sets cover with probability (much) less than $1-\gamma$. Consequently, there is interest in identifying these troublemaker parameter values. Recent efforts along these lines in the sparse normal mean model include Belitser, (2017), Belitser and Nurushev, (2017), van der Pas et al., 2017b , Castillo and Szabó, (2019), and Nurushev and Belitser, (2019); for details beyond the normal mean model, see, e.g., Szabó et al., (2015), Belitser and Ghosal, (2019), and Martin and Tang, (2019).

In this paper we pursue the question of valid posterior uncertainty quantification in the sparse normal mean model in several new ways. First, we focus our investigation on a relatively new type of posterior distribution based on a suitable empirical prior. This approach differs considerably from classical empirical Bayes and has been shown to have strikingly good practical and theoretical performance; see Martin and Walker, (2019) and the references in Section 2. That this method would also provide valid uncertainty quantification was conjectured by Martin, (2017) and our work here confirms that. Second, while previous investigations into posterior uncertainty quantification have focused primarily on credible $\ell_{2}$-balls for the full $\theta$ vector, we consider an arguably more practical question of marginal credible intervals for certain features of $\theta$. From a general Bernstein–von Mises result for the full posterior in Theorem 3, we give sufficient conditions for valid uncertainty quantification about a general linear functional of $\theta$ in Corollary 1. For a specific linear functional, however, such as $\theta_{1}$, one might ask if the conditions can be weakened. To this end, we show that the nominal coverage can be achieved even if the true $|\theta_{1}^{\star}|$ is slightly smaller than what is needed for selection consistency or the Bernstein–von Mises theorem; see Theorem 4. To our knowledge, there are no results—at least none that are positive—concerning the coverage probability properties of marginal credible intervals for the horseshoe or any other prior distributions. Third, we give numerical results to compare the performance of credible sets based on our empirical priors and those based on the horseshoe prior. These results confirm what the theory suggests, namely, that there is a non-trivial range of non-zero $\theta_{1}^{\star}$ values that the former can properly cover while the latter cannot. Finally, in Section 5, we provide some concluding remarks and directions for future research.

2 An empirical Bayes model

2.1 Prior and posterior construction

Under assumption of sparsity, it makes sense to re-express the vector $\theta$ as a pair $(S,\theta_{S})$, where $S\subseteq\{1,2,\ldots,n\}$, the configuration of $\theta$, indicates which coordinates are non-zero, and $\theta_{S}$ is an $|S|$-vector that contains $S$-specific non-zero values; here $|S|$ is the cardinality of $S$. With this decomposition, a hierarchical prior is natural, i.e., a prior for the configuration $S$ and then, given $S$, a prior for the $S$-specific parameter $\theta_{S}$. The sparsity assumption provides some relevant prior information to help with the choice of the prior on $S$. On the other hand, no prior information is available for $\theta_{S}$, so a “default” prior is recommended. It turns out that certain features of this default prior—in particular, its tails—can significantly affect the properties of the corresponding posterior. Indeed, using normal priors for $\theta_{S}$ can lead to suboptimal posterior concentration properties, whereas heavier-tailed priors like Laplace have much better performance (e.g., Castillo and van der Vaart, 2012). Unfortunately, while the normal priors are conjugate and lead to relatively simple posterior computations, the heavy-tailed priors make computation more difficult. So one faces a dilemma: use a simple conjugate prior that is easy to compute but may have suboptimal theoretical performance, or pay the non-trivial computational price for the use of a theoretically justified prior. Is such a choice really necessary?

Martin and Walker, (2014) argued that the prior tails cannot have much of an effect if its center is appropriately chosen. In other words, an appropriately-centered conjugate normal prior should not suffer the same suboptimality property as if it had a fixed center. Of course, an “appropriate” center can only originate from the data, hence an empirical prior. The particular form of our empirical prior, $\Pi_{n}$, for $\theta=(S,\theta_{S})$ is

[TABLE]

where $\pi(S)$ is a marginal prior for the configuration $S$ to be specified below, $Y_{S}$ is the sub-vector of $Y$ corresponding to configuration $S$, and $\tau>0$ is a prior precision factor which, as we discuss later, will be taken relatively small. Again, the idea is to properly center the prior based on the data, so that the thin normal tails will not affect the asymptotic concentration properties. A regression version of this formulation is given in Martin et al., (2017), and Martin and Walker, (2019) describe a general empirical prior framework.

Before getting into the specifics of the empirical prior formulation, we pause here to say a few words about the philosophy behind such an approach. As we see it, in modern Bayesian analysis, especially that of high-dimensional problems, the prior is treated primarily as an input that can be tuned and tweaked in order to obtain a posterior distribution which is “good” in one ore more senses. For example, the horseshoe prior of Carvalho et al., (2010) and its various extensions are popular nowadays, not because they capture genuine subjective beliefs, but because the corresponding posterior is (believed to be) fast to compute and achieves desirable asymptotic concentration rates and uncertainty quantification properties.333In fact, it is not uncommon in seminar talks or less formal discussions to hear one motivate the construction of a new prior by saying that existing priors “don’t work” and/or the new prior “works better,” not that it more accurately reflects subjective prior beliefs, etc. If the prior is chosen primarily for the frequentist properties that the corresponding posterior satisfies, then there is no reason not to consider using a data-dependent prior if that too can produce a posterior with as good or better properties. A goal of this paper, and the related literature on empirical priors, is to show that, indeed, they achieve these desirable properties and, at least in some cases, outperform the state-of-the-art Bayesian methods based on fixed priors.

Returning to the construction, one has some flexibility in the choice of marginal prior for the configuration and, in our numerical investigations in Section 4, we will consider two such priors. Both decompose the prior $\pi$ for $S$ into a marginal prior, $f_{n}(s)$, for the size, $s=|S|$, of the configuration, and a uniform prior over the space of configurations of a given size. Below are the respective mass functions, $f_{n}$, for the two priors.

• •

The complexity prior, described in Castillo and van der Vaart, (2012) and used in Martin et al., (2017) for regression, has mass function $f_{n}$ given by

[TABLE]

for constants $a,c>0$. This is just a truncated geometric prior with success probability equal to $1-(cn^{a})^{-1}$ so, depending on $a$, it strongly penalizes complex configurations, i.e., those models with a large number of non-zero means.

• •

The beta–binomial prior in Martin and Walker, (2014) is of the form

[TABLE]

where $b_{n}\to\infty$ as $n\to\infty$. This corresponds to a beta prior, ${\sf Beta}(b_{n},1)$, on a latent variable $W$ and, given $W=w$, a binomial prior, $\mathsf{Bin}(n,1-w)$, for $|S|$. The mean of $W$ is $b_{n}(b_{n}+1)^{-1}$, which makes the prior mean of $|S|$ close to zero so, like (1), the prior in (2) also severely penalizes complex configurations.

Finally, the empirical Bayes posterior distribution for $\theta$, denoted by $\Pi^{n}$, is obtained by combining the prior with a fractional power of the likelihood function according to Bayes’s theorem. That is,

[TABLE]

where $\alpha\in(0,1)$ and $L_{n}(\theta)\propto\exp\{-\frac{1}{2\sigma^{2}}\|Y-\theta\|^{2}\}$, with $\|\cdot\|$ the $\ell_{2}$-norm. The power $\alpha$ does not affect the normality of the likelihood, $L_{n}(\theta)$, it only inflates the variance, hence making its contours wider. In fact, given the normal form of the conditional priors, the posterior actually takes a relatively simple form; see Section 2.2 below.

As for the fractional power, the results in Belitser and Nurushev, (2017) and Belitser and Ghosal, (2019) suggest that taking $\alpha=1$ might be possible, perhaps with some adjustments elsewhere, but we believe there are reasons to retain this flexibility, especially in the context of uncertainty quantification. In particular, existing results on coverage probability of credible sets require a blow-up factor—e.g., the factor $L$ in Equation (3) of van der Pas et al., 2017b —to inflate the credible set beyond the size determined by the posterior distribution itself. As is well-known in the generalized Bayes community (e.g., Grünwald and van Ommen, 2017; Syring and Martin, 2019; Holmes and Walker, 2017), the inclusion of $\alpha$ makes the posterior distribution wider and, therefore, also has a beneficial blow-up effect on the credible sets; see Remark 1.

2.2 Posterior computation

As advertised, the empirical prior leads to relatively simple posterior computation. In general, since the prior for $\theta_{S}$, given $S$, is normal, it is possible to analytically integrate out $\theta_{S}$ to obtain a marginal posterior distribution444The expression for $\pi^{n}(S)$ given in Section 4.1 of Martin et al., (2017) has a typo, but the correct formula is given in the supplement at https://arxiv.org/abs/1406.7718. for the configuration $S$, i.e.,

[TABLE]

In situations where only characteristics of the marginal posterior of $S$ are needed, e.g., for feature selection, certain credible sets (see Section 3), etc, then the above formula can be used to design a straightforward Metropolis–Hastings algorithm to sample $S$; Liu et al., (2018) adopt a shotgun stochastic search based on (4). Since $\theta_{S}$, given $S$, is normal, the Metropolis–Hastings procedure can be easily augmented to sample both $S$ and $\theta_{S}$.

With the beta–binomial prior, there is additional structure that can be used to design a posterior sampling algorithm. That is, by introducing that latent variable “$W$,” the full conditionals—$\theta\mid W$ and $W\mid\theta$—are available in closed-form, hence a Gibbs sampler, as in Martin and Walker, (2014), can be readily employed.

2.3 Asymptotic concentration properties

Let $\theta^{\star}\in\mathbb{R}^{n}$ denote the true mean vector, with configuration $S^{\star}=S_{\theta^{\star}}$, where the notation “$S_{\theta}$” refers to the configuration of the vector $\theta$. We assume that $\theta^{\star}$ is sparse and, in particular, $s^{\star}=|S^{\star}|$ is $o(n)$ as $n\to\infty$. The results reviewed below show that the posterior, $\Pi^{n}$, is able to optimally identify both $\theta^{\star}$ and $S^{\star}$ as $n\to\infty$.

The only aspect of the model described above that is not completely determined is the prior $f_{n}$ on the configuration size. The key is that the tails of $f_{n}$ have a certain rate of decay, as described in Equation (2.2) of Castillo and van der Vaart, (2012) and Assumption 1 of Castillo et al., (2015). In particular, we assume here that

[TABLE]

where $(K_{1},K_{2})$ and $(a_{1},a_{2})$ are suitable constants; in particular $a_{1}>a_{2}$. This implies that $f_{n}(s+t)/f_{n}(s)$ is lower and upper bounded by $(K_{j}n^{-a_{j}})^{t}$ for $j=1,2$, respectively, which is critical for the existing methods of proof. Clearly, the complexity prior (1) satisfies this and, moreover, any prior that does so must be similar to that in (1). It turns out that the beta–binomial prior also satisfies (5) if $b_{n}\propto n^{\xi}$ for some $\xi>1$.

Theorem 1**.**

If the prior for the configuration size satisfies (5), then the posterior $\Pi^{n}$ adaptively attains the minimax concentration rate, i.e., for a constant $M^{\prime}>0$,

[TABLE]

Proof.

Follows from the arguments in Martin et al., (2017). ∎

Besides concentration of the posterior around $\theta^{\star}$, it is interesting to consider the posterior for the configuration. Ideally, it would put all of its mass on $S^{\star}$, asymptotically, but this requires some conditions on the magnitudes of the non-zero means; see Result 3 in Theorem 2 below. These selection results will be used in Section 3.

Theorem 2**.**

For the pair $a_{1}>a_{2}$ in (5), assume that $s^{\star}n^{-a_{2}}\to 0$

No proper supersets.* $\mathsf{E}_{\theta^{\star}}\Pi^{n}(\{\theta:S_{\theta}\supset S^{\star}\})\to 0$.* 2. 2.

No important misses.* Define the threshold*

[TABLE]

where $M=\max\{M^{\prime},1+a_{1}\}$ for the constant $M^{\prime}$ as in Theorem 1. Define $S^{\dagger}=\{i:|\theta_{i}^{\star}|>\rho_{n}\}\subseteq S^{\star}$; then $\mathsf{E}_{\theta^{\star}}\Pi^{n}(\{\theta:S_{\theta}\not\supseteq S^{\dagger}\})\to 0$. 3. 3.

Selection consistency.* If $S^{\dagger}=S^{\star}$, i.e., if all the non-zero means have magnitude larger than $\rho_{n}$ in (6), then $\mathsf{E}_{\theta^{\star}}\pi^{n}(S^{\star})\to 1$.*

Proof.

Part 1 is a consequence of the analysis in Martin et al., (2017); a proof of Part 2 is given in Appendix A below. Part 3 follows immediately from Parts 1–2. ∎

3 Uncertainty quantification

3.1 Bernstein–von Mises phenomenon

The key to proving that the posterior uncertainty quantification is valid is establishing some suitable approximation of the posterior in terms of familiar quantities whose distribution is known. A Bernstein–von Mises theorem, for example, would show that the posterior is approximately normal. Therefore, posterior credible sets must closely resemble the corresponding sets for the normal approximation, and properties of the mean and variance of the normal distribution can be used to obtain coverage results. In our present case, under suitable conditions, a Bernstein–von Mises theorem is relatively easy to establish, and it says that the posterior $\Pi^{n}$ can be approximated, asymptotically, by a product of normal distribution on the true configuration, $S^{\star}$, and a point mass at 0 on its complement, $S^{\star c}$. Compare this result to Theorem 6 in Castillo et al., (2015).

Theorem 3**.**

Let $d_{\text{\sc tv}}$ denote the total variation distance and $S^{\star}=S_{\theta^{\star}}$. Then

[TABLE]

for all $\theta^{\star}$ such that $\mathsf{E}_{\theta^{\star}}\pi^{n}(S^{\star})\to 1$, where $v_{\alpha}=\sigma^{2}(\alpha+\tau)^{-1}$ and $\delta_{S^{\star c}}$ denotes the point mass distribution for $\theta_{S^{\star c}}$ concentrated at the origin.

Proof.

Let $D_{n}=D_{n}(Y)$ denote the total variation distance in (7). By definition,

[TABLE]

where the supremum is over all Borel measurable subsets of $\mathbb{R}^{n}$. Next, again by definition,

[TABLE]

so we can immediately upper-bound the absolute value on the right-hand side of (8) by

[TABLE]

The absolute difference above is 0 when $S=S^{\star}$ and bounded by 2 otherwise, so

[TABLE]

After taking expectation of both sides, if $\theta^{\star}$ is such that $\mathsf{E}_{\theta^{\star}}\pi^{n}(S^{\star})\to 1$, then the upper bound vanishes, thus proving the theorem. ∎

Remark 1*.*

Recall that $v_{\alpha}=\sigma^{2}(\alpha+\tau)^{-1}$. Therefore, in order for the variance in the posterior distribution to be at least as large as the variance of $Y_{S^{\star}}$, we need $v_{\alpha}\geq\sigma^{2}$ or, equivalently, $\alpha+\tau\leq 1$. The other properties of the posterior require $\alpha<1$ and $\tau>0$ so we do not have the option to choose $\alpha=1$ and $\tau=0$ at this point. We do have the option to take $\tau=1-\alpha$, in which case $v_{\alpha}=\sigma^{2}$. However, in the spirit of conservatism and validity, we prefer to take $\alpha+\tau$ slightly less than 1. This agrees with Grünwald and Mehta, (2017) who say that Bayesian inference is “safe” if the learning rate—their version of $\alpha$—is strictly less than a critical threshold. In particular, for our simulations in Section 4, we take $\alpha=0.95$ and $\tau=0.025$, so that $\alpha+\tau=0.975<1$.

3.2 Credible sets and their coverage properties

The literature on Bayesian uncertainty quantification in sparse high-dimensional problems has focused primarily on the coverage of $\ell_{2}$-balls centered at the posterior mean. Here we consider a different and arguably more practical question concerning the coverage of marginal credible intervals for certain one-dimensional features of $\theta$. In particular, we consider a general linear functional $\psi=x^{\top}\theta$, where $x\in\mathbb{R}^{n}$ is a specified vector. This includes, as an important special case, inference about an individual mean, say, $\theta_{1}$, and we expect that other functionals of $\theta$, which are only approximately linear, could likely be handled in a similar way.

Without loss of generality, we consider coverage probability of upper credible bounds for $\psi$. Towards this, we need its marginal posterior distribution function, call it $H_{n}(t)=\Pi^{n}(\{\theta:x^{\top}\theta\leq t\})$, which is given by

[TABLE]

where $\Phi$ is the standard normal distribution function, $x_{S}$ is the $|S|$-vector corresponding to configuration $S$, $\hat{\psi}_{S}=x_{S}^{\top}Y_{S}$, and $v_{\alpha}=\sigma^{2}(\alpha+\tau)^{-1}$ as before.

For the special case where $x=e_{k}$, the $k^{\text{th}}$ standard basis vector, this reduces to

[TABLE]

where $p_{k}^{n}=\pi^{n}(S\ni k)$ is the posterior inclusion probability of index $k$, i.e., the posterior probability that the configuration contains the particular index $k$, a quantity often used in Bayesian feature selection problems (e.g., Barbieri and Berger, 2004). The above expression reveals that the marginal posterior of $\theta_{k}$ is, as expected, a mixture of a point mass at zero and a normal distribution centered at $Y_{k}$. The only setback is that the inclusion probability, $p_{k}^{n}$, is a complicated function of the full data.

Towards uncertainty quantification, fix a nominal confidence level $\gamma\in(0,\frac{1}{2})$. Then the $100(1-\gamma)$% credible upper bound for $\psi=x^{\top}\theta$, denoted by $t_{\gamma}$, is given by

[TABLE]

The somewhat complicated expression here is to handle the discontinuity in $H_{n}$ at the origin. That is, if the inclusion probability is small, there will be a large jump at the origin, in which case, $H_{n}(t)=1-\gamma$ has no solution, but the infimum ensures credible upper bound will be $t_{\gamma}=0$. Otherwise, if there is no large jump at the origin, then $H_{n}(t)=1-\gamma$ has a unique solution, and that would be $t_{\gamma}$. Note that this credible bound is precisely the one that practitioners would read off from the corresponding quantiles of the posterior samples, different from the marginal posterior mean plus an inflated standard error as studied in van der Pas et al., 2017b and elsewhere.

Of course, the credible bound $t_{\gamma}=t_{\gamma}(Y)$ depends on the data, so, for the posterior uncertainty quantification to be valid in a frequentist sense, the relevant property is that the coverage probability is close to the nominal $1-\gamma$ level, i.e.,

[TABLE]

It turns out that this coverage probability property is a more-or-less immediate consequence of the Bernstein–von Mises phenomenon.

Corollary 1**.**

If $\alpha+\tau\leq 1$, so that $v_{\alpha}\geq\sigma^{2}$ then (7) implies (11) and, therefore, valid uncertainty quantification.

Proof.

Intuitively, if $\theta\sim\Pi^{n}$ is approximately normal, then $\psi=x^{\top}\theta$ is approximately normal too, with mean $\hat{\psi}_{S^{\star}}$ and variance $v_{\alpha}\|x_{S^{\star}}\|^{2}$. More formally, an argument identical to that in the proof of Theorem 3 gives

[TABLE]

where $\Pi_{\psi}^{n}$ is the derived posterior distribution of $\psi=x^{\top}\theta$. This convergence in total variation implies that the posterior quantile, $t_{\gamma}$, asymptotically agrees with the corresponding quantile for $\mathsf{N}(\hat{\psi}_{S^{\star}},v_{\alpha}\|x_{S^{\star}}\|^{2})$. The latter quantile and its distributional properties are familiar and, therefore, (7) implies (11) if $v_{\alpha}\geq 1$. ∎

To summarize so far, if the posterior distribution can correctly identify $S^{\star}$, the true configuration, i.e., if $\pi^{n}(S^{\star})\to 1$, then the corresponding posterior uncertainty quantification about any linear functional will be valid in a frequentist sense. We know, from Theorem 2, that if all the non-zero means have magnitude exceeding the threshold $\rho_{n}$ in (6), then $\pi^{n}(S^{\star})\to 1$. But it is worth asking if this can be established under weaker conditions, at least in certain cases. One such case, which is practically relevant, is a marginal credible interval for $\psi=\theta_{k}$, say, and we ask how large does $|\theta_{k}^{\star}|$ need to be in order for the marginal posterior uncertainty quantification to be valid. This is directly related to the numerical experiments considered in Section 4 below.

At least intuitively, for the marginal posterior for $\psi=\theta_{k}$, whether the uncertainty quantification is valid hinges entirely on the behavior of the inclusion probability. For example, if $\theta_{k}^{\star}\neq 0$, and if $p_{k}^{n}\to 1$, then the coverage probability of the $100(1-\gamma)$% marginal posterior credible upper bound should be approximately $1-\gamma$. The following theorem confirms this intuition.

Lemma 1**.**

The $100(1-\gamma)$% marginal posterior credible interval for $\theta_{k}$ is valid, i.e., the frequentist coverage probability is at least $1-\gamma$ for sufficiently large $n$, if $\theta_{k}^{\star}=0$ and $\mathsf{E}_{\theta^{\star}}(p_{k}^{n})\to 0$ or if $\theta_{k}^{\star}\neq 0$ and $\mathsf{E}_{\theta^{\star}}(1-p_{k}^{n})\to 0$.

Proof.

Start with the case $\theta_{k}^{\star}=0$. If $t_{\gamma}$ is the corresponding upper credible limit, then the relevant coverage probability can be bounded as follows:

[TABLE]

The first inequality above is due to the fact that $H_{n}(0)\geq 1-p_{k}^{n}$; the second is by the reverse Markov inequality. By assumption, $\mathsf{E}_{\theta^{\star}}(p_{k}^{n})\to 0$ and, therefore, the coverage probability converges to 1, as expected.

Next, without loss of generality, we consider $\theta_{k}^{\star}>0$. In this case, the relevant coverage event is $\{H_{n}(\theta_{k}^{\star})\leq 1-\gamma\}$ where, in this case,

[TABLE]

with $V_{k}=\Phi\{v_{\alpha}^{-1/2}(Y_{k}-\theta_{k}^{\star})\}$. Note, if $\alpha+\tau\leq 1$, then $v_{\alpha}\geq\sigma^{2}$, which implies that $V_{k}$ has tails no thinner than $\mathsf{Unif}(0,1)$, the uniform distribution between 0 and 1. Below we show that the non-coverage probability is upper-bounded by $\gamma$ as $n\to\infty$.

The non-coverage probability is

[TABLE]

Intuitively, since $p_{k}^{n}$ will be close to 1, and $1-V_{k}$ has tails no heavier than $\mathsf{Unif}(0,1)$, we expect the right-most probability above to be less than $\gamma$. To make this formal, introduce a $b\in(0,1)$ and write

[TABLE]

By Markov’s inequality,

[TABLE]

If $\mathsf{E}_{\theta^{\star}}(1-p_{k}^{n})\to 0$, as we assumed, then we get

[TABLE]

Since this holds for all $b\in(0,1)$, it holds for the infimum over $b$ so, therefore, the limiting non-coverage probability is no bigger than $\gamma$, which proves the claim. ∎

To control the inclusion probability in the case $\theta_{k}^{\star}\neq 0$, it will be necessary to distinguish the means whose magnitudes exceed the cutoff $\rho_{n}$ in (6) and those that are just simply non-zero, so recall, from Section 2.3, $S^{\dagger}=\{i:|\theta_{i}^{\star}|>\rho_{n}\}$ and let $s^{\dagger}=|S^{\dagger}|$. Here we are specifically considering the possibility that $|\theta_{k}^{\star}|$ is less than $\rho_{n}$, so we have $s^{\dagger}<s^{\star}$. The following theorem gives a lower bound on the magnitude of $\theta_{k}^{\star}$ such that its inclusion probability, $p_{k}^{n}$, converges to 1, and that bound depends on how large $s^{+}$ is compared to $s^{\star}$. And, in light of Lemma 1, this establishes how large the mean needs to be in order for the marginal posterior distribution to provide valid uncertainty quantification.

Theorem 4**.**

If the conditions of Theorem 2 hold, and if $|\theta_{k}^{\star}|>\zeta_{n}$, where

[TABLE]

and $s^{\dagger}<s^{\star}$ is the number of non-zero means have magnitude exceeding the threshold in (6), then the inclusion probability satisfies $\mathsf{E}_{\theta^{\star}}(p_{k}^{n})\to 1$ as $n\to\infty$ and, consequently, for $\alpha+\tau\leq 1$, the marginal posterior for $\theta_{k}$ provides valid uncertainty quantification.

Proof.

Without loss of generality, take $\sigma^{2}=1$. By Theorem 2, we know that those $s^{\dagger}$ means are sufficiently large that the posterior will not exclude them. Moreover, we know that no non-zero means will be included. So the only configurations that need to be considered are those that contain all of the $s^{\dagger}$ large means and any of the $s^{\star}-s^{\dagger}$ non-zero but not-large means. Let $\mathcal{S}^{0}$ denote this subset of configurations. Then the exclusion probability for $\theta_{k}$ is asymptotically equal to

[TABLE]

For a generic $S$ in this sum, we can bound the posterior probability,

[TABLE]

where the equality is based on plugging in the expression for $\pi^{n}$ in (4). Taking expectation and using the bound on $\theta_{k}^{\star}$ in the lemma’s statement, we get

[TABLE]

Plugging in the definition of the prior $\pi$ for $S$, using the assumed bound in (5), and writing $s=|S|$, we can simplify the right-hand side above:

[TABLE]

Summing over all $S\in S^{+}$ with $S\not\ni k$, we get the upper bound

[TABLE]

It is easy to check that the remaining sum is upper-bounded by

[TABLE]

Finally, we get

[TABLE]

and it can be immediately seen that if $\zeta_{n}$ is as in (12), then the exclusion probability for $\theta_{k}$ vanishes. Moreover, the claimed frequentist validity of the marginal posterior credible interval follows from Lemma 1. ∎

One special case of the above setting is where all the $s^{\star}$ non-zero means exceed the threshold in (6) except for $\theta_{k}^{\star}$, so that $s^{\dagger}=s^{\star}-1$. In that case, the threshold (12) is

[TABLE]

Since $\log(\tfrac{n-s^{\star}+1}{s^{\star}})$ is strictly less than $\log n$, and the $M$ in (6) is larger than $1+a_{1}$, there is some improvement in the sense that $|\theta_{k}^{\star}|$ can actually be smaller than the threshold from Theorem 2 and still achieve the desired coverage probability. However, since $s^{\star}=o(n)$, the difference is not substantial enough to impact the threshold’s rate, apparently the improvement only shows up in the constant. But compare this to Theorem 2.1 in van der Pas et al., 2017b , where they only have negative results about the coverage probability for marginal credible sets under the horseshoe prior when the non-zero means are smaller than a threshold like in (6), even with a so-called “blow-up factor.” But, most importantly, we show numerically in Section 4, that this difference between the horseshoe and our empirical priors shows up in finite-samples.

3.3 Extensions

There are two immediate directions in which one might like to generalized the above results: to the regression setting, where the mean vector is $\theta=X\beta$, with $X$ a $n\times p$ matrix of covariates and $\beta$ a sparse $p$-vector of coefficients, possibly with $p\gg n$, and/or the error variance, $\sigma^{2}$, is unknown. For the former case, Martin et al., (2017) establish posterior concentration rates for a version of the empirical prior (on $\beta$) considered here, along with selection consistency. For this regression problem, the key technical obstacle is separating $\beta$ from $X\beta$, which requires use of a so-called “sparse singular value” of $X$; see Equation (13) in Martin et al., (2017), which is very similar to Equation (11) in Arias-Castro and Lounici, (2014) and Definition 2.3 in Castillo et al., (2015). These results for the regression problem where used by Martin and Tang, (2019)—much like we used Theorems 1–2 here—to establish a Bernstein–von Mises theorem and to investigate coverage properties of marginal credible intervals, in particular, for the purpose of out-of-sample prediction. For the second direction, unknown $\sigma^{2}$, so far there has been no theoretical investigation into the performance of these posteriors based on an empirical prior. The simulations presented in Martin and Tang, (2019), using a either a plug-in estimator of or a prior for $\sigma^{2}$, show good empirical performance. So we fully expect that the theoretical properties of the posterior distribution would not be affected by unknown $\sigma^{2}$ either, but this still needs to be worked out. We believe that this is doable because, at least with a conjugate prior for $\sigma^{2}$, there would be a closed-form expression for the marginal posterior distribution of $S$, which is what our theoretical analysis here focused on. So, extending the results here would boil down to analyzing a different version of the “$\pi^{n}(S)$” expression considered here.

4 Numerical investigations

In this section we investigate the numerical performance of the posterior distributions based on the two styles of empirical priors described in Section 2.1 compared to those based on the horseshoe prior introduced in Carvalho et al., (2010). In particular, we consider the following three methods:

• HS.

For the horseshoe method we employ the horseshoe package in R (van der Pas et al., 2016) where maximum marginal likelihood, based on their function HS.MMLE, is used to estimate the global scale parameter.

• EB1.

The approach in Martin and Walker, (2014) based on the beta–binomial prior $f_{n}$. We use exactly the settings in that paper and the R code they provided.

• EB2.

The approach in Martin et al., (2017) based on the complexity prior $f_{n}$. We use the R code they provided, specialized from the regression to the means problem, and we take $\alpha=0.95$ and $\tau=0.025$, so that $\alpha+\tau<1$; see Remark 1.

Here we assume that $\sigma^{2}=1$ and all three methods use this known value.

Our goal here is to compare the performance of marginal credible intervals for the three approaches described above. The specific setting we consider here, similar to that in Section 2 of van der Pas et al., 2017b , assumes the first five means are relatively large, namely, $\theta_{1}^{\star}=\cdots=\theta_{5}^{\star}=7$, the second five means are of intermediate magnitude, namely, $\theta_{6}^{\star}=\cdots=\theta_{10}^{\star}=2$, $\theta_{11}^{\star}$ will vary, and the remaining $\theta_{12}^{\star},\cdots,\theta_{n}^{\star}$ are 0. We want to know how large does $\theta_{11}^{\star}$ have to be in order for the (two-sided) marginal credible interval for $\theta_{11}$ to attain the nominal frequentist coverage probability, which we take as 95%, i.e., $\gamma=0.05$. Figure 1 plots the empirical coverage probability and average lengths of the marginal credible intervals for $\theta_{11}$, as a function of the signal size $\theta_{11}^{\star}$, for two different values of $n$, namely, $n=200$ and $n=500$. These are Monte Carlo estimates based on 500 data sets at each value of $\theta_{11}^{\star}$ along a grid from 0 to 10.

The primary question we seek to answer here is for what range of $\theta_{11}^{\star}$ values is the nominal 95% coverage rate attained? The theoretical results in Section 3 indicated that the credible sets based on the empirical prior might have an advantage over the horseshoe-based credible sets in this respect. So our conjecture is that the former will have a slightly broader range of nominal-coverage-attainment than the latter. The results in Figure 1 confirm this conjecture, that is, both the EB1 and EB2 coverage curves reach the nominal level well before HS for $n=200$ and $n=500$. An interesting surprise is that EB1, based on the beta–binomial prior in Martin and Walker, (2014) performs considerably better in terms of coverage than EB2 based on the complexity prior. In our opinion, the Gibbs sampler for EB1 is a more elegant approach than the Metropolis–Hastings implementation of EB2, but it is not clear why the former would outperform the latter in terms of coverage probabilities. And the length curves in Panels (b) and (d) indicate that the EB gains over HS in coverage are not due to excessive length.

Finally, despite being based on a two-groups or spike-and-slab formulation, both EB implementations were much faster to compute than HS based on the horseshoe package. In fact, we also ran the above experiment for $n=1000$ and could easily produce results for EB1 and EB2, which look similar to those in Figure 1, but the HS computations were prohibitively slow and regularly produced errors, hence not reported here.

5 Conclusion

In this paper, we focus on uncertainty quantification derived from posterior distributions based on a class of empirical priors. This is certainly a relevant question to ask about such methods, given their double-use of the data. That is, while estimation may not be negatively affected by the use of data in both the prior and the likelihood, it is possible that such an “informative” prior would make the posterior spread too narrow and, consequently, the coverage probability would fall below the nominal level. Our investigation here revealed that this actually is not the case; in fact, in a certain sense, the empirical prior leads to improved uncertainty quantification compared to a more traditional approach where the prior is free of data. Follow-up investigations should look into other aspects of uncertainty quantification, such as credible balls for the full $\theta$ vector with respect to various distances, and extension of these results to the regression setting and other structured high-dimensional problems.

A unique feature of the empirical priors in this normal mean model is that the conditional prior variance for $\theta_{S}$, given $S$, is constant, it does not depend on data or even the sample size. In general, however, the empirical priors constructed in Martin and Walker, (2019) would have prior spread vanishing at some rate with $n$. For example, in the regression setting, after writing the coefficient vector as $\beta=(S,\beta_{S})$, the conditional empirical prior for $\beta_{S}$, given $S$, is $\mathsf{N}_{|S|}(\hat{\beta}_{S},\sigma^{2}\tau^{-1}(X_{S}^{\top}X_{S})^{-1})$, where $X_{S}$ is the $n\times|S|$ sub-matrix of the full $n\times p$ matrix $X$, and $\hat{\beta}_{S}$ is the least-squares estimate corresponding to $X_{S}$. Note that each entry in this variance is $O(n^{-1})$. Despite the relatively tight prior concentration, the method in Martin et al., (2017) performs well in terms of estimation and variable selection. That the prior variance is consistent with that of the sampling distribution of $\hat{\beta}_{S}$, the prior mean, suggests that uncertainty quantification will not be impacted by the prior concentration. Indeed, in Martin and Tang, (2019), the focus was on posterior credible sets for out-of-sample prediction in the sparse high-dimensional regression problem and, in addition to the theoretical results, their simulation study reveals properties similar those demonstrated here, namely, that the empirical Bayes credible sets achieve the desired coverage and are no less efficient than the horseshoe, even with unknown error variance. Moreover, in a monotone density estimation context, Martin, (2018) constructed an empirical prior, whose spread was vanishing with $n$, and his simulation experiments revealed that the coverage of posterior credible sets was not negatively affected by the prior concentration and, in fact, the coverage performance was better in some cases than the proper Bayesian posterior credible sets from the Dirichlet process-based formulation in Salomond, (2014). Full justification of these claims in this and other nonparametric problems, however, requires further numerical and theoretical investigations.

Acknowledgments

The authors thank the editors of the special issue of Sankhya A dedicated to Jayanta K. Ghosh for the invitation to contribute, and the anonymous reviewers for their helpful suggestions that improved both our results and presentation. This work is partially supported by the National Science Foundation, DMS–1737933.

Appendix A Proof of Theorem 2

The proof strategy here closely follows that of Theorem $5^{\prime}$ in the supplement to Martin et al., (2017) presented in the most recent arXiv version (arXiv:1406.7718). To fix notation, let $S^{\star}$ be the true configuration of size $s^{\star}=|S^{\star}|$, and let $S^{\dagger}\subseteq S^{\star}$ be the set of all $i$ such that $|\theta_{i}^{\star}|\geq\rho_{n}$, where $\rho_{n}$ is as in (6), and write $s^{\dagger}=|S^{\dagger}|$.

Based on Theorem 2 in Martin et al., (2017), we can restrict to configurations $S$ such that $|S|\leq Cs^{\star}$, where $s^{\star}=|S^{\star}|$ and $C$ is a large constant. Take such an $S$ that also satisfies $S\not\supseteq S^{\dagger}$. Then $\pi^{n}(S)$ can be bounded as follows:

[TABLE]

where $z=(1+\alpha\tau^{-1})^{-1/2}<1$. A key observation is that

[TABLE]

and the latter two terms are independent since they depend on disjoint sets of $Y_{i}$’s. Therefore, using this independence and the familiar central and non-central chi-square moment generating functions, we get

[TABLE]

By definition of $S^{\dagger}$, and the fact that $1+\alpha>1$, the above expectation can be upper-bounded by

[TABLE]

Putting the pieces together we have

[TABLE]

We want to sum this over all $S\not\supseteq S^{\dagger}$ but, since it only involves size of $S$, we only need to sum over sizes. Indeed, after plugging in the definition of $\pi(S)$ we get

[TABLE]

For the binomial coefficient ratio we have the following simplification and bound:

[TABLE]

Next, to bound the double-sum, split it into two parts:

[TABLE]

We need to show that both parts on the right-hand side above vanish as $n\to\infty$. For the first double-sum we have

[TABLE]

Since $M>1+a_{1}$, the inner sum is $O(1)$ and the outer sum—because there is a common $n^{-(M-a_{1}-1)}$ factor—is $o(1)$ as $n\to\infty$. Similarly, for the second double-sum we have

[TABLE]

The inner sum is $O(1)$ and, since $a_{2}<a_{1}$, the outer sum is upper-bounded by $O(s^{\star}n^{-a_{2}})$ which goes to 0 by assumption. Both terms in the double-sum above vanish as $n\to\infty$, thus proving the claim.