# Asymptotically optimal inference in sparse sequence models with a simple
data-dependent measure111This work is partially supported by the U.S. National
Science Foundation, DMS–1811802.

Ryan Martin222Department of Statistics, North Carolina State University,
<EMAIL_ADDRESS>
###### Abstract

For high-dimensional inference problems, statisticians have a number of
competing interests. On the one hand, procedures should provide accurate
estimation, reliable structure learning, and valid uncertainty quantification.
On the other hand, procedures should be computationally efficient and able to
scale to very high dimensions. In this note, I show that a very simple data-
dependent measure can achieve all of these desirable properties
simultaneously, along with some robustness to the error distribution, in
sparse sequence models.

_Keywords and phrases:_ high-dimensional inference; concentration rate;
structure learning; uncertainty quantification; variational approximation.

##  1 Introduction

Dating back at least to Stein, (1956, 1981), a fundamental problem in
statistical inference is that of estimating a high-dimensional mean vector
under additive noise. Specifically, suppose that the observable data is
$Y^{n}=(Y_{1},\ldots,Y_{n})^{\top}$ with posited model

$Y_{i}=\theta_{i}+Z_{i},\quad i=1,\ldots,n,$ (1)

where the errors, $Z_{1},\ldots,Z_{n}$, are independent and identically
distributed (iid) with mean zero and subgaussian tails; see Section 3 for
specifics. While I don’t assume the error distribution to be normal, I do
assume that it’s fully known, with no parameters to be estimated. The goal is
inference on $\theta=(\theta_{1},\ldots,\theta_{n})^{\top}$.

In practical applications using models like in (1), e.g., Efron, (2010), the
index $n$ is generally quite large, which makes this a high-dimensional
inference problem. As is often the case in such problems, certain structural
assumption on $\theta=(\theta_{1},\ldots,\theta_{n})^{\top}$ are necessary in
order to achieve accurate estimation. Here, the structural assumption I’ll
consider is sparsity—most of the entries in the $\theta$ vector are zero. This
notion of sparsity is consistent with the science: in genomics applications,
the biology says that only a relatively small number of genes would be
associated with a particular phenotype. Mathematically weaker notions of
sparsity are possible, ones that don’t assume any exact zeros, which I discuss
below in Section 5.

There now many different approaches to this problem. On the non-Bayesian side,
the focus has been on developing shrinkage and thresholding estimators, and
notable references include James and Stein, (1961), Efron and Morris, (1973),
Donoho and Johnstone, (1994), and Abramovich et al., (2006), among others. On
the Bayesian side, the strategy is to construct suitable sparsity-inducing
priors. This can be done using continuous shrinkage priors (e.g., Bhattacharya
et al., 2015; Bhadra et al., 2017; Carvalho et al., 2010), or with spike-and-
slab priors (e.g., Castillo and van der Vaart, 2012; Johnstone and Silverman,
2004; Martin and Walker, 2014). Inspired by the developments in Belitser,
(2017) and Belitser and Nurushev, (2020), here I want to take a different
approach. In particular, I construct a so-called data-dependent measure which
is based on a simple approximation to a more complicated (empirical) Bayes
formulation developed in Martin et al., (2017) and specialized to the present
context by Martin and Ning, (2020). Computation of this data-dependent measure
requires no Markov chain Monte Carlo (MCMC) or optimization—it has a simple,
intuitive, and explicit form. Most importantly, its theoretical properties are
optimal, or nearly so, in every relevant sense:

  * •

it concentrates asymptotically around the true sparse mean vector
$\theta^{\star}$ at the minimax optimal rate, adaptive to the unknown sparsity
level;

  * •

its mean vector is an asymptotically, adaptively minimax estimator;

  * •

it concentrates asymptotically on a subspace of $\mathbb{R}^{n}$ whose
dimension is roughly the same as the effective dimension of $\theta^{\star}$;

  * •

it consistently selects the non-zero entries of $\theta^{\star}$;

  * •

and it provides asymptotically valid uncertainty quantification.

The majority of these properties are familiar, perhaps only the last one
requires further explanation. From this data-dependent measure, one can
readily extract credible sets that capture a specified fraction of the
distribution’s mass. A question is if that credible set is also asymptotically
a confidence set in the sense that it covers the true $\theta^{\star}$ with
probability approximately equal to that same fraction. This has been a hot
topic in the Bayesian literature recently (e.g., van der Pas et al., 2017;
Szabó et al., 2015; Rousseau and Szabo, 2020; Castillo and Szabó, 2020), and
here I establish this property for my simple data-dependent measure. The ideas
and techniques developed in Belitser, (2017) and Belitser and Nurushev, (2020)
will prove to be useful in this regard.

The remainder of the paper is organized as follows. Section 2 gives some
background on an approach for high-dimensional inference that uses empirical
or data-driven priors. The data-dependent measure investigated here originally
arose as a variational approximation to this empirical Bayes-style posterior
distribution. Section 3 investigates the asymptotic concentration properties
of the data-dependent measure, providing justification for the claims made in
the bulleted list above. Some comments about computation of the data-dependent
measure are given in Section 4, along with an illustration to show that its
theoretical and computational simplicity don’t come at the cost of statistical
inefficiency in finite samples. Concluding remarks are given in Section 5 and
all the proofs are presented in Appendix A.

##  2 A simple data-dependent measure

###  2.1 Empirical prior and its posterior

While I’m not doing so here in this paper, it is common to assume that the
$Y_{i}$’s are independent and normally distributed, i.e., $Y_{i}\sim{\sf
N}(\theta_{i},\sigma^{2})$, for $i=1,\ldots,n$, with known variance
$\sigma^{2}$. This distributional assumption determines a likelihood function
and, in turn, a Bayes or empirical Bayes approach can be taken. One such
approach developed recently is that in Martin and Walker, (2014), which makes
us of an empirical or data-dependent prior, that is, a prior distribution that
directly depends on the data. The details here are based on the developments
in Martin et al., (2017) for the high-dimensional linear regression problem;
see Martin and Walker, (2019) for some general theory on empirical priors and
properties of the corresponding posterior distributions.

In the approach of Martin and Walker, (2019)—see Martin and Ning, (2020) for
the sequence model special case—the prior is formulated by first expressing
$\theta$ as $(S,\theta_{S})$, where $S\subseteq\\{1,2,\ldots,n\\}$ represents
the configuration of zeros and non-zeros and $\theta_{S}$ the configuration-
specific parameters. Then the (empirical) prior for $\theta$, denoted by
$\Pi_{n}$, with subscript “$n$” to indicate data dependence, is specified
hierarchically as follows.

  * •

Marginal prior for $S$: The the prior on the size/cardinality $|S|$ of $S$ has
a mass function $f_{n}$, and the conditional prior for $S$, given its
cardinality $|S|=s$ is uniform over the $\binom{n}{s}$ subsets of size $s$.
Several different forms of $f_{n}$ are discussed in Martin and Ning, (2020),
but I’ll not need these for what I plan to do here.

  * •

Conditional prior for $\theta_{S}$, given $S$: With the configuration $S$
given, it’s determined that $\theta_{i}=0$ for $i\not\in S$, so only a prior
for $\theta_{S}=\\{\theta_{i}:i\in S\\}$ is necessary. For that prior, the
choice in Martin and Ning, (2020) is

$(\theta_{S}\mid S)\sim{\sf N}_{|S|}(Y_{S},\sigma^{2}\gamma^{-1}I_{|S|}),$

where $\sigma^{2}$ is the assumed value of the variance in the normal model,
$\gamma>0$ is a prior hyperparameter, and $I_{|S|}$ is the $|S|\times|S|$
identity matrix. Note that this conditional prior depends on data through the
centering around $Y_{S}=\\{Y_{i}:i\in S\\}$.

With this empirical prior and the normal likelihood, Martin and Ning, (2020)
described construction of a corresponding posterior distribution as

$\Pi^{n}(d\theta)\propto L_{n}^{\alpha}(\theta)\,\Pi_{n}(d\theta),$

where
$L_{n}(\theta)=(2\pi\sigma^{2})^{-n/2}\exp\\{-\frac{1}{2\sigma^{2}}\|Y-\theta\|^{2}\\}$
is the likelihood function and $\alpha\in(0,1)$ is a fractional power. An
explanation of the role played by $\alpha$ is given in Martin and Walker,
(2019) and the references therein. Again, I’ll not be using this posterior so
it’s not necessary for me to give an explanation of $\alpha$ here.

###  2.2 An approximation

Computation of the posterior $\Pi^{n}$ requires MCMC, but this can be done
rather efficiently, as demonstrated in Martin and Ning, (2020) and van Erven
and Szabó, (2020), compared to the proposed Monte Carlo computations for the
horseshoe and other priors. The same is true for the high-dimensional
regression version; see, e.g., Martin et al., (2017) and Martin and Tang,
(2020). However, as shown in Ray and Szabó, (2020) and elsewhere, simple and
computationally efficient variational approximations of these high-dimensional
posterior are possible. This inspired Yang and Martin, (2020) to develop a
corresponding variational approximation for the empirical prior formulation.

The jumping off point for Yang and Martin, (2020) was the recognition that the
empirical prior described above can be written in a very simple form. Indeed,
the prior assumes that $\theta_{1},\ldots,\theta_{n}$ are independent and the
respective marginal distributions are

$\theta_{i}\sim\lambda_{n}\,{\sf
N}(Y_{i},\sigma^{2}\gamma^{-1})+(1-\lambda_{n})\,\delta_{0},\quad
i=1,\ldots,n,$ (2)

where $\lambda_{n}$ is the prior inclusion probability, which depends on the
sample size $n$ but not on the particular index $i$. It’s relatively easy to
show that $\lambda_{n}=n^{-1}\mathsf{E}|S|$, where $\mathsf{E}|S|$ is the
prior mean for $|S|$ under $f_{n}$. For the two kinds of prior mass function
$f_{n}$ they considered, it follows that, for a constant $a>0$,

$\lambda_{n}=O(n^{-(a+1)}),\quad n\to\infty.$ (3)

That $\lambda_{n}$ is vanishing with $n$—and faster than $n^{-1}$—is
consistent with the idea that $\theta$ is believed to be sparse. For
simplicity in what follows, I’ll take $\lambda_{n}=n^{-(1+a)}$.

When the prior can be expressed in the basic spike-and-slab form (2), it is
not too difficult to derive a variational approximation to the posterior
distribution. Indeed, consider a mean-field approximation family of the form

$\bigotimes_{i=1}^{n}\\{\phi_{i}\,{\sf
N}(\mu_{i},\tau_{i}^{2})+(1-\phi_{i})\,\delta_{0}\\}.$

This corresponds to independent components, each being a mixture of a normal
and a point mass at 0, but with component-specific parameters. The variational
approximation proceeds by finding the set of parameters
$\\{(\mu_{i},\tau_{i}^{2},\phi_{i}):i=1,\ldots,n\\}$ that minimizes the
Kullback–Leibler divergence of the posterior distribution $\Pi^{n}$ from the
above family. As Yang and Martin, (2020) show, there are closed-form
expressions—no update equations as is typical—for those “best” parameters:

$\displaystyle\mu_{i}$ $\displaystyle=y_{i}$ $\displaystyle\tau_{i}^{2}$
$\displaystyle=\sigma^{2}(\alpha+\gamma)^{-1}$ (4)
$\displaystyle\operatorname{logit}(\phi_{i})$
$\displaystyle=\operatorname{logit}(\lambda_{n})+\tfrac{1}{2}\log\tfrac{\gamma}{\alpha+\gamma}+\tfrac{\alpha}{2\sigma^{2}}y_{i}^{2}.$

Here, $\operatorname{logit}(\phi)=\log\\{\phi/(1-\phi)\\}$. Note that the
variance component $\tau_{i}^{2}$ is actually the same for each $i$. Moreover,
the weights $\phi_{i}$ that indicate whether a $\theta_{i}$ is zero or non-
zero are increasing in $y_{i}^{2}$, as one would expect.

###  2.3 My proposal

My proposal here in this paper is to define a data-dependent measure

$\Delta^{n}=\bigotimes_{i=1}^{n}\\{\phi_{i}\,{\sf
N}(\mu_{i},\tau_{i}^{2})+(1-\phi_{i})\,\delta_{0}\\},$ (5)

with the specific parameters (4) plugged in. This is different from the
perspective in Yang and Martin, (2020) because I’m not starting with a normal
model, constructing a posterior distribution based on an empirical prior, and
then developing a variational approximation. Instead, I’m directly defining
$\Delta^{n}$ as the data-dependent measure I intend to use for inference on
$\theta$. Consequently, there are no choices about priors to be explained or
claims of the variational approximation’s accuracy to be justified. Whether
$\Delta^{n}$ is a reasonable procedure to use rests entirely on what
properties it possesses and, as I show next in Section 3, its properties are
optimal in every practical respect.

##  3 Asymptotic properties

###  3.1 Concentration rates

To set the scene, recall that the errors $Z_{1},\ldots,Z_{n}$ are iid copies
of a random variable $Z$ whose distribution is known. Moreover, I’ll assume
that $Z$ is subgaussian in the sense that the moment-generating function of
$Z$ satisfies

$\mathsf{E}\exp(tZ)\leq\exp(\sigma^{2}t^{2}/2),\quad\text{all
$t\in\mathbb{R}$},$ (6)

where $\sigma$ is a scale parameter, often called the variance proxy, such
that $\sigma^{2}\geq\mathsf{V}(Z)$. Of course, this covers the Gaussian case,
but there are other examples too, including bounded random variables; see,
e.g., Boucheron et al., (2013). One key property of subgaussian random
variables is that they have exponential tail probability bounds; see Appendix
A.1 below. In addition, the tails are sufficiently thin to ensure that the
moment-generating function of $(Z/\sigma)^{2}$ exists in an interval that
contains the origin. What’s relevant to the analysis here is the upper
endpoint of that interval, which I’ll denote as $T>0$. The actual endpoint
depends on the specific form of the $Z$ distribution, e.g., if $Z$ is
Gaussian, then $T=1/2$. However, the largest bound that I’m aware of that
covers all subgaussian cases simultaneously is $T=1/4$ (Honorio and Jaakkola,
2014, Appendix B).

As indicated above, interest is in cases where the true mean vector
$\theta^{\star}$ is sparse, so I need to make this structural assumption
precise. For a generic vector $\theta\in\mathbb{R}^{n}$, let $S_{\theta}$
denote its configuration, i.e., $S_{\theta}=\\{i:\theta_{i}\neq 0\\}$. Let
$|S_{\theta}|$ denote the cardinality of the configuration. Of course,
$|S_{\theta^{\star}}|\leq n$ but I have in mind cases where the inequality is
strict, even/especially cases where $|S_{\theta^{\star}}|\ll n$, since these
are the only ones in which accurate estimation of $\theta^{\star}$ is
possible. To characterize this notion of accuracy, recall that the minimax
rate (e.g., Donoho et al., 1992) relative to $\ell_{2}$-error is

$\varepsilon_{n}^{2}(\theta^{\star})=|S_{\theta^{\star}}|\log(en/|S_{\theta^{\star}}|).$
(7)

That is, every estimator $\hat{\theta}$ satisfies

$\sup_{\theta^{\star}}\varepsilon_{n}^{-2}(\theta^{\star})\mathsf{E}_{\theta^{\star}}\|\hat{\theta}-\theta^{\star}\|^{2}\geq
1.$

So we say that an estimator is minimax optimal if equality holds up to a
constant,

$\sup_{\theta^{\star}}\varepsilon_{n}^{-2}(\theta^{\star})\mathsf{E}_{\theta^{\star}}\|\hat{\theta}-\theta^{\star}\|^{2}\lesssim
1,$ (8)

where “$\lesssim$” denotes inequality up to a universal constant. The
following two results show that the data-dependent measure $\Delta^{n}$ in
(5), with suitable choice of $\alpha$, has this same asymptotic concentration
rate property.

######  Theorem 1.

For the data-dependent measure $\Delta^{n}$ in (5), suppose that $\lambda_{n}$
and $\alpha$ in (4) satisfy, respectively, (3) and $\alpha<2T$, for $T$
determined by the subgaussian error distribution. For
$\varepsilon_{n}^{2}(\theta^{\star})$ defined in (7) and any sequence
$M_{n}>0$ with $M_{n}\to\infty$,

$\sup_{\theta^{\star}}\mathsf{E}_{\theta^{\star}}\Delta^{n}(\\{\theta\in\mathbb{R}^{n}:\|\theta-\theta^{\star}\|^{2}>M_{n}\varepsilon_{n}^{2}(\theta^{\star})\\})\to
0,\quad n\to\infty.$

######  Theorem 2.

Under the conditions of Theorem 1, the mean vector $\hat{\theta}$ derived from
the data-dependent measure $\Delta^{n}$ in (5) satisfies (8).

An important observation is that the data-dependent measure $\Delta^{n}$ is
not aware of the sparsity level $|S_{\theta^{\star}}|$ of the true
$\theta^{\star}$, and yet it concentrates at that specific optimal rate. This
feature is commonly referred to as adaptation, i.e., the concentration rate of
$\Delta^{n}$ is adaptive to the unknown sparsity level of $\theta^{\star}$
that determines the optimal rate.

Of course, the data-dependent measure’s (nearly) optimal concentration rate is
a plus, but this property alone doesn’t imply that $\Delta^{n}$ is learning
the low-dimensional structure in $\theta^{\star}$. The next result
demonstrates that indeed the data-dependent measure is learning that structure
in the sense that the dimension of the space on which $\Delta^{n}$
concentrates is roughly the same as the effective dimension
$|S_{\theta^{\star}}|$ of $\theta^{\star}$.

######  Theorem 3.

Under the conditions of Theorem 1, for any $M_{n}>1$ such that
$M_{n}\to\infty$, the data-dependent measure $\Delta^{n}$ in (5) satisfies

$\sup_{\theta^{\star}}\mathsf{E}_{\theta^{\star}}\Delta^{n}(\\{\theta\in\mathbb{R}^{n}:|S_{\theta}|>M_{n}|S_{\theta^{\star}}|\\})\to
0,\quad n\to\infty.$

###  3.2 Structure learning

Theorem 3 established that $\Delta^{n}$ concentrates on a subspace of roughly
the effective dimension of $\theta^{\star}$, but more can be said. The next
result shows that, asymptotically, $\Delta^{n}$ will not assign positive mass
to proper supersets of $S_{\theta^{\star}}$. To ensure that all the signals
are detectable, an additional assumption about the magnitude of those non-zero
$\theta_{i}^{\star}$ values is needed. Specifically, consider

$\min_{i\in S_{\theta^{\star}}}|\theta_{i}^{\star}|\geq
H:=\bigl{(}\tfrac{2\sigma^{2}K}{\alpha}\log n\bigr{)}^{1/2},\quad\text{for
some $K>2+a$},$ (9)

where $\sigma^{2}$ is the variance proxy of the error distribution, $\alpha$
is as in (4), and $a$ is as in (3). Up to constants, (9) is equivalent to the
“beta-min condition” common in the high-dimensional estimation literature.

######  Theorem 4.

Under the conditions of Theorem 1, the data-dependent measure $\Delta^{n}$
satisfies

$\sup_{\theta^{\star}}\mathsf{E}_{\theta^{\star}}\Delta^{n}(\\{\theta\in\mathbb{R}^{n}:S_{\theta}\supset
S_{\theta^{\star}}\\})\to 0,\quad n\to\infty.$

Moreover, if $\theta^{\star}$ is such that (9) holds, then

$\sup_{\theta^{\star}}\mathsf{E}_{\theta^{\star}}\Delta^{n}(\\{\theta\in\mathbb{R}^{n}:S_{\theta}\not\supseteq
S_{\theta^{\star}}\\})\to 0,\quad n\to\infty.$

If all the above conditions hold, then the two conclusions can be combined,
giving

$\mathsf{E}_{\theta^{\star}}\Delta^{n}(\\{\theta\in\mathbb{R}^{n}:S_{\theta}=S_{\theta^{\star}}\\})\to
1,\quad n\to\infty.$ (10)

This theorem says that, asymptotically, $\Delta^{n}$ will not support
configurations that contain zero entries in $\theta^{\star}$. Moreover, if the
non-zero entries in $\theta^{\star}$ are sufficiently large, in the sense of
(9), then $\Delta^{n}$ will not support configurations that miss any of those
non-zero entries either. In the latter case, the only option is that
$\Delta^{n}$ asymptotically supports the true configuration
$S_{\theta^{\star}}$, hence it effectively learns the low-dimensional
structure in $\theta^{\star}$. The result in (10) is often referred to as a
selection consistency property, since any reasonable selection procedure based
on $\Delta^{n}$, e.g.,

$\hat{S}=\arg\max_{S}\delta^{n}(S)\quad\text{or}\quad\hat{S}=\\{i:\phi_{i}>0.5\\},$

will, for large enough $n$, identify the correct $S_{\theta^{\star}}$.

###  3.3 Uncertainty quantification

An important question is if inferences derived from the data-dependent measure
are reliable in the sense that they control the frequency of errors, at least
asymptotically. There are a number of ways this can be assessed in the present
context. One is to consider certain one-dimensional summaries of the
$n$-dimensional vector $\theta$, in particular, linear functionals. Yang and
Martin, (2020) considered this when treating $\Delta^{n}$ as a variational
approximation of the full (empirical) Bayes posterior under a normal model.
Another angle is to consider a credible set for the full $n$-dimensional
vector. This approach has been considered for a variety of different kinds of
Bayes and empirical Bayes posterior distributions in the literature, e.g., van
der Pas et al., (2017), Belitser, (2017), Belitser and Nurushev, (2020), and
Belitser and Ghosal, (2019). Here I’m going to derive analogous properties for
the simple data-dependent measure $\Delta^{n}$ in (5) under the same general
subgaussian error structure as above.

Recall that $\hat{\theta}$ is the mean vector of the data-dependent measure
$\Delta^{n}$. Define a ball centered around $\hat{\theta}$ with radius $r>0$
as

$B_{n}(r)=\\{\theta\in\mathbb{R}^{n}:\|\theta-\hat{\theta}\|\leq r\\}.$

The goal is to choose a data-dependent radius $r=\hat{r}$ such that the ball
approximately achieves a target coverage probability $1-\zeta$ and has near-
optimal size. There are two natural strategies for selecting the radius
$\hat{r}$ based on $\Delta^{n}$. The first is motivated by ensuring the ball
has sufficient probability under $\Delta^{n}$, while the second is motivated
by achieving the optimal size.

  1. 1.

Quantile-based. Set
$\hat{r}=\inf\\{r:\Delta^{n}(\theta:\|\theta-\hat{\theta}\|\leq r)\geq
1-\zeta\\}$.

  2. 2.

Plug-in estimator-based. Define $\hat{S}=\\{i:\phi_{i}>\tfrac{1}{2}\\}$ and
set $\hat{r}^{2}=|\hat{S}|\log(en/|\hat{S}|)$.

The following theorem summarizes the coverage probability and size properties
of the two corresponding credible balls.

######  Theorem 5.

Assume the conditions of Theorem 4 hold. Also, let
$\Theta_{n}\subset\mathbb{R}^{n}$ denote the set where the condition (9) on
the minimum signal size holds. Fix a significance level
$\zeta\in(0,\frac{1}{2})$ and a threshold $\eta>0$.

  1. 1.

Let $\hat{r}$ denote the $\Delta^{n}$ quantile-based radius defined above.
Then there exists constants $L$ and $M$ such that, for all sufficiently large
$n$,

$\sup_{\theta^{\star}\in\Theta_{n}}\mathsf{P}_{\theta^{\star}}\\{B_{n}(Mg_{n}\hat{r})\not\ni\theta^{\star}\\}\leq\zeta\quad\text{and}\quad\sup_{\theta^{\star}}\mathsf{P}_{\theta^{\star}}\\{\hat{r}^{2}>L\varepsilon_{n}^{2}(\theta^{\star})\\}\leq\eta,$

where the inflation factor $g_{n}$ satisfies $g_{n}=\log(en)$.

  2. 2.

Let $\hat{r}$ denote the plug-in estimator-based radius defined above. Then
there exists constants $L$ and $M$ such that, for all sufficiently large $n$,

$\sup_{\theta^{\star}\in\Theta_{n}}\mathsf{P}_{\theta^{\star}}\\{B_{n}(M\hat{r})\not\ni\theta^{\star}\\}\leq\zeta\quad\text{and}\quad\sup_{\theta^{\star}}\mathsf{P}_{\theta^{\star}}\\{\hat{r}^{2}>L\varepsilon_{n}^{2}(\theta^{\star})\\}\leq\eta.$

More-or-less explicit expressions for the constants $(L,M)$ are given in the
proof, so one could technically use these values for practical implementation.
However, I make no claims that these constants are optimal, in fact, it’s
likely that they’re conservative. In any case, the point is simply to say that
the data-dependent measure spread is “right” in the sense that credible balls
with slightly larger than optimal size can achieve the nominal coverage
probability. The additional inflation factor $g_{n}$ in the quantile-based
credible ball is needed because, apparently, the quantile itself is too small
by a logarithmic factor. Similar inflation factors have been needed by other
authors proving analogous results (e.g., Belitser and Ghosal, 2019).

##  4 Numerical illustration

Computation of the data-dependent measure $\Delta^{n}$ in (5) is trivial and
fast. Virtually every summary has a closed-form expression, so it’s
straightforward to produce the mean $\hat{\theta}$, to select a set of
“active” variables via $\hat{S}=\\{i:\phi_{i}>\tfrac{1}{2}\\}$, and to extract
marginal credible intervals for each $\theta_{i}$. This can be done almost
instantaneously, far faster than the computations using the horseshoe package
in R (van der Pas et al., 2016) and at least most of those methods compared in
van Erven and Szabó, (2020).

There are a host of available methods that provide high-quality estimation and
structure learning. The theoretical support for uncertainty quantification
using the simple data-dependent measure is the chief novelty here, so that’s
what I’ll focus on here. Although the theory presented above is for the joint
credible ball, there is good reason (e.g., Yang and Martin, 2020; Martin and
Ning, 2020) to believe that the corresponding marginal credible intervals
would be approximately valid too. So my objective in this section is simply to
show that the very fast computations do not come at the expense of validity or
efficiency. That is, this simple data-dependent measure produces marginal
credible intervals which are as good or better than those from other methods
sharing the same theoretical guarantees but with heavier computational burden.

Specifically, I redo the simulation study presented in Martin and Ning, (2020)
comparing the coverage probability and mean length of the horseshoe and two
empirical prior-based credible intervals. Let $n=500$ and suppose that the
errors are iid standard normal. Similar to Section 2 of van der Pas et al.,
(2017), consider a case where the first five entries of $\theta^{\star}$ are
relatively large, i.e., $\theta_{1}^{\star}=\cdots=\theta_{5}^{\star}=7$, the
second five are intermediate, i.e.,
$\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.
Of interest is to see how large $\theta_{11}^{\star}$ needs to be in order for
the coverage probability to be approximately equal to 0.95, the nominal level.
Figure 1 plots the empirical coverage probability and mean lengths of the four
marginal credible intervals for $\theta_{11}$, as a function of the signal
size $\theta_{11}^{\star}$. Of course, the coverage probability will be low
when the signal size is small, so of primary interest is how quickly the
coverage probability climbs to near 0.95 as $\theta_{11}^{\star}$ increases.
It’s clear that the data-dependent measure and the empirical prior formulation
of Martin and Walker, (2014) perform comparably in the sense that both get to
the target coverage probability by around $\theta_{11}^{\star}\approx 6$,
before the other two. Interestingly, these two are no less efficient in terms
of interval lengths, since they’re all close to the optimal “$2\times
1.96=3.92$” length marked by the horizontal line on the plot.

Finally, the data-dependent measure can easily scale to far bigger $n$, e.g.,
$n\sim 10^{6}$, while the other methods would have serious difficulties with
problems of this size. A thorough comparison of the proposed data-dependent
measure with other fast methods for this problem (e.g., van Erven and Szabó,
2020; Ray and Szabó, 2020; Ročková and George, 2018) would be an interesting
direction to pursue.

(a) Coverage probability

(b) Mean length

Figure 1: Plots of the coverage probability and mean length of the marginal
credible intervals for $\theta_{11}$, as a function of $\theta_{11}^{\star}$,
based on the four methods: horseshoe (HS), a beta–binomial empirical prior
(EB1 Martin and Walker, 2014), a complexity-driven empirical prior (EB2,
Martin et al., 2017), and the data-dependent measure (DDM).

##  5 Conclusion

In this paper, I’ve considered inference on a sparse, high-dimensional mean
vector using a simple data-dependent measure. Assuming only subgaussianity of
the error distribution, I was able to show that the data-dependent measure has
optimal asymptotic convergence properties in virtually every respect. Most
notably, I was able to establish that the data-dependent measure provides
asymptotically valid uncertainty quantification in the sense that credible
balls centered around the data-dependent measure’s mean vector, with suitable
data-driven choices of radius, can achieve the nominal coverage probability
while maintaining roughly the optimal size. My proofs of the various results
presented herein are relatively straightforward thanks to the simple form of
the data-dependent measure under investigation. Moreover, the simple form
makes computation, even for very large-scale problems, fast and easy. The
numerical illustration in Section 4 shows, however, that the theoretical and
computational simplicity don’t come at the expense of statistical efficiency
or poor finite-sample performance.

One possible extension is to explore the asymptotic behavior of this data-
dependent measure under assumptions on $\theta^{\star}$ that are
mathematically weaker than my notion of sparsity here, e.g., under the so-
called excessive bias restriction in Belitser, (2017). I chose to work here
with sparsity because it’s a simpler and more intuitive condition, but I
expect that more general results are possible, even with only subgaussianity
assumptions.

Other kinds of low-dimensional structure in the mean vector can likely be
handled using similarly simple data-dependent measures. For example, in a
sequence model where the mean has a piecewise constant structure (e.g., van
der Pas and Rockova, 2017; Liu et al., 2020), sparsity shows up in the
successive differences, so things would not be too much different from the
case considered here. Clearly there are limitations to how far this kind of
simple approach could go, but it’s interesting and practically useful to find
where the boundary is. In the high-dimensional regression problem, for
example, the theoretical support available in, say, Ray and Szabó, (2020), for
mean-field variational approximations suggests that other simple and more
directly defined data-dependent measures could have similar properties. It’s
perhaps not surprising that posterior concentration rate results could be
achieved even if the inherent correlation in the full posterior distribution
is ignored, but it would be interesting to see if other structure learning or
uncertainty quantification properties were similarly unaffected.

It’s a mathematical fact that credible sets in high-dimensions can’t be both
valid confidence sets and of adaptively optimal size. The approach that most
investigations have taken, including mine here, is to start with a data-
dependent measure that achieves the optimal size property and show that its
credible sets approximately achieve the target coverage probability too. To
me, the most interesting take-away message from this paper is that apparently
very simple solutions can achieve this “optimal size and approximate coverage”
property. If both simple and not-so-simple solutions can achieve the same
properties, then arguably the standard is too low. How might we approach these
structured high-dimensional problems in a more discriminating way? One idea is
to turn the line of reasoning around, that is, to start with something that
achieves valid uncertainty quantification and think about how to introduce the
assumed structure, e.g., sparsity, in such a way that efficiency is gained but
validity isn’t lost. The approach I have in mind, with developments underway,
is to start with a valid inferential model (e.g., Martin and Liu, 2013, 2015),
treat the assumed structure as genuine but incomplete prior information
encoded as an imprecise probability, and combine with the inferential model
output in an appropriate way that preserves validity. The main difference
between this and the standard approach taken in this paper is that validity is
given higher priority than efficiency, which I believe to more appropriate for
scientific investigations.

##  Appendix A Proofs

###  A.1 Preliminary results

The only distributional assumption being made here is that the errors
$Z_{1},\ldots,Z_{n}$ in (1) are iid copies of a random variable $Z$ with
subgaussian tails. As the name suggests, this condition implies that $Z$ has
some Gaussian-like properties. Here I collect a few relevant facts about
subgaussian random variables that will be used in what follows.

  * •

It is well known that the square of a subgaussian random variable is
subexponential. I won’t need any specific properties of subexponential random
variables, so there’s no need to give a formal definition. All that matters
here is that subexponential random variables have a moment-generating function
in an interval that contains the origin. In particular,

$\mathsf{E}e^{t(Z/\sigma)^{2}}\lesssim 1,\quad t\in(0,T].$ (11)

Moreover, since translations don’t effect the tails of a distribution, the
moment-generating function of $(Z+u)^{2}$, for any $u$, also exists for some
arguments, in particular, when the argument is negative, I get

$\mathsf{E}e^{-t(Z+u)^{2}/\sigma^{2}}\lesssim e^{-tu^{2}/\sigma^{2}},\quad
t>0.$ (12)

  * •

An equivalent definition of subgaussian random variables is that they admit an
exponential tail probability bound just like the Gaussian. In particular,

$\mathsf{P}(|Z|>t)\leq 2e^{-t^{2}/2\sigma^{2}},\quad t>0.$ (13)

Next are two results that’ll be needed in the proofs of the theorems below.
These make use of the properties for subgaussian random variables described
above.

######  Lemma 1.

If $\lambda_{n}$ satisfies (3) and $\alpha<2T$, then

$\sum_{i=1}^{n}\mathsf{E}_{\theta_{i}^{\star}}\phi_{i}\leq|S^{\star}|+o(1)\lesssim|S^{\star}|.$

###### Proof.

First, split the sum as

$\sum_{i=1}^{n}\mathsf{E}_{\theta_{i}^{\star}}\phi_{i}=\sum_{i\in
S^{\star}}\mathsf{E}_{\theta_{i}^{\star}}\phi_{i}+\sum_{i\not\in
S^{\star}}\mathsf{E}_{\theta_{i}^{\star}}\phi_{i}.$

Since $\phi_{i}\leq 1$, the sum over $i\in S^{\star}$ is clearly
$\leq|S^{\star}|$. For the sum over $i\not\in S^{\star}$, note that all the
means are zero and, therefore, all the terms in the sum are the same, i.e.,

$\sum_{i\not\in
S^{\star}}\mathsf{E}_{\theta_{i}^{\star}}\phi_{i}=(n-|S^{\star}|)\mathsf{E}_{0}\phi_{i}.$

Again, since $\phi_{i}\leq 1$, it follows that

$\mathsf{E}_{0}\phi_{i}\leq\mathsf{E}_{0}e^{\operatorname{logit}(\phi_{i})}=\xi_{n}\,\mathsf{E}e^{\frac{\alpha}{2}(Z/\sigma)^{2}},$

where

$\xi_{n}=\exp\\{\operatorname{logit}(\lambda_{n})+\tfrac{1}{2}\log\tfrac{\gamma}{\alpha+\gamma}\\}.$

Since $\alpha<2T$ by assumption, it follows from (11) that

$\mathsf{E}_{0}\phi_{i}\lesssim\xi_{n}\lesssim\exp\\{\operatorname{logit}(\lambda_{n})\\}=n^{-(1+a)},$

which implies

$\sum_{i\not\in S^{\star}}\mathsf{E}_{\theta_{i}^{\star}}\phi_{i}\lesssim
n^{-a}=o(1),\quad\text{as $n\to\infty$}.$

Combining this with the bound from the sum over $i\in S^{\star}$ completes the
proof. ∎

######  Lemma 2.

If the $\lambda_{n}$ in (5) satisfies (3), then

$\mathsf{E}_{\theta^{\star}}\int\|\theta-\theta^{\star}\|^{2}\,\Delta^{n}(d\theta)\lesssim\varepsilon_{n}^{2}(\theta^{\star}).$

###### Proof.

By the definition of $\Delta^{n}$ in (5), it’s easy to check that

$\displaystyle\int\|\theta-\theta^{\star}\|^{2}\,\Delta^{n}(d\theta)$
$\displaystyle=\sum_{i=1}^{n}\int(\theta_{i}-\theta_{i}^{\star})^{2}\,\Delta^{n}(d\theta)$
$\displaystyle=\sum_{i=1}^{n}\Bigl{\\{}\phi_{i}\int(\theta_{i}-\theta_{i}^{\star})^{2}\,{\sf
N}(\theta_{i}\mid
Y_{i},\tau_{i}^{2})\,d\theta_{i}+(1-\phi_{i})\theta_{i}^{\star 2}\Bigr{\\}}$
$\displaystyle=\sum_{i=1}^{n}\bigl{[}\phi_{i}\\{\tau_{i}^{2}+(Y_{i}-\theta_{i}^{\star})^{2}\\}+(1-\phi_{i})\theta_{i}^{\star
2}\bigr{]}$ $\displaystyle=\sum_{i=1}^{n}\tau_{i}^{2}\phi_{i}+\sum_{i\not\in
S^{\star}}\phi_{i}Y_{i}^{2}+\sum_{i\in
S^{\star}}\\{\phi_{i}(Y_{i}-\theta_{i}^{\star})^{2}+(1-\phi_{i})\theta_{i}^{\star
2}\\}.$

Note that $\tau_{i}^{2}$ are constant in $i$ and do not depend on data, this
can come outside the same (and the following expectation). Taking expectation,
as using the fact that $\phi_{i}\leq 1$, the right-hand side is bounded by

$\tau^{2}\sum_{i=1}^{n}\mathsf{E}\phi_{i}+\sum_{i\not\in
S^{\star}}\mathsf{E}\phi_{i}Y_{i}^{2}+\sigma^{2}|S^{\star}|+\sum_{i\in
S^{\star}}\theta_{i}^{\star 2}\mathsf{E}(1-\phi_{i}).$

I’ll deal with each term in this sum separately. The first term is
$\lesssim|S^{\star}|$ by Lemma 1. Second, consider
$\mathsf{E}_{\theta_{i}^{\star}}\phi_{i}Y_{i}^{2}$ for $i\not\in S^{\star}$,
which means $\theta_{i}^{\star}=0$. For $x>0$ to be determined,

$\mathsf{E}_{0}\phi_{i}Y_{i}^{2}=\mathsf{E}_{0}\phi_{i}Y_{i}^{2}1_{|Y_{i}|\leq
x}+\mathsf{E}_{0}\phi_{i}Y_{i}^{2}1_{|Y_{i}|>x}.$

The first term on the right-hand side is $x^{2}\mathsf{E}_{0}\phi_{i}\lesssim
x^{2}n^{-(a+1)}$, as shown in the proof of Lemma 1. The second term is bounded
by $\mathsf{E}_{0}Y_{i}^{2}1_{|Y_{i}|>x}=\mathsf{E}Z^{2}1_{|Z|>x}$. For this,
we can use the tail probability bound (13) for $Z$ as follows:

$\displaystyle\mathsf{E}Z^{2}1_{|Z|>x}$
$\displaystyle\int_{0}^{\infty}\mathsf{P}(Z^{2}1_{|Z|>x}>t)\,dt$
$\displaystyle=\int_{0}^{\infty}\mathsf{P}\\{|Z|>\max(x,t^{1/2})\\}\,dt$
$\displaystyle=\int_{0}^{x^{2}}\mathsf{P}(|Z|>x)\,dt+\int_{x^{2}}^{\infty}\mathsf{P}(|Z|>t^{1/2})\,dt$
$\displaystyle\leq
2x^{2}e^{-x^{2}/2\sigma^{2}}+2\int_{x^{2}}^{\infty}e^{-t/2\sigma^{2}}\,dt$
$\displaystyle\lesssim(x^{2}+1)e^{-x^{2}/2\sigma^{2}}.$

Take $x=\\{2\sigma^{2}\log(n/|S^{\star}|)\\}^{1/2}$, so that

$\mathsf{E}_{0}\phi_{i}Y_{i}^{2}\lesssim
n^{-(a+1)}\log(n/|S^{\star}|)+|S^{\star}|n^{-1}\log(en/|S^{\star}|).$

Summing over $i\not\in S^{\star}$ gives

$\sum_{i\not\in S^{\star}}\mathsf{E}\phi_{i}y_{i}^{2}\lesssim
n^{-a}\log(n/|S^{\star}|)+|S^{\star}|\log(en/|S^{\star}|)\lesssim\varepsilon_{n}^{2}(\theta^{\star}).$
(14)

Lastly, for the third term, recall that

$1-\phi_{i}=1-\\{1+\xi_{n}^{-1}e^{-(\alpha/2\sigma^{2})Y_{i}^{2}}\\}^{-1},$

where $\xi_{n}\propto\exp\\{-\operatorname{logit}(\lambda_{n})\\}$. Since
$z\mapsto(1+z)^{-1}$ is convex, Jensen’s inequality says

$\mathsf{E}_{\theta_{i}^{\star}}(1-\phi_{i})\leq
1-\\{1+\xi_{n}^{-1}\mathsf{E}_{\theta_{i}^{\star}}e^{-(\alpha/2\sigma^{2})Y_{i}^{2}}\\}^{-1}.$

By (12), the expectation satisfies

$\mathsf{E}_{\theta_{i}^{\star}}e^{-(\alpha/2\sigma^{2})Y_{i}^{2}}=\mathsf{E}e^{-(\alpha/2)(Z+\theta_{i}^{\star})^{2}/\sigma^{2}}\leq
ce^{-(\alpha/2)\theta_{i}^{\star 2}/\sigma^{2}},$

for a constant $c>0$. Therefore,

$\mathsf{E}_{\theta_{i}^{\star}}(1-\phi_{i})\leq\frac{c\xi_{n}^{-1}e^{-k\theta_{i}^{\star
2}}}{1+c\xi_{n}^{-1}e^{-k\theta_{i}^{\star 2}}},$

where $k=\alpha/2\sigma^{2}$. Multiplying by $\theta_{i}^{\star 2}$ gives

$\theta_{i}^{\star
2}\,\mathsf{E}_{\theta_{i}^{\star}}(1-\phi_{i})\leq\frac{c\xi_{n}^{-1}\theta_{i}^{\star
2}e^{-k\theta_{i}^{\star 2}}}{1+c\xi_{n}^{-1}e^{-k\theta_{i}^{\star 2}}}.$

As a function of $\theta_{i}^{\star 2}$, this has the form of a gamma density
with shape parameter 2 and rate parameter $k$. Such a density has mode
$k^{-1}$. Plugging in that mode, what’s left is a bounded sequence in $n$, so
the right-hand side above is $\lesssim 1$, which implies

$\sum_{i\in S^{\star}}\theta_{i}^{\star
2}\,\mathsf{E}_{\theta_{i}^{\star}}(1-\phi_{i})\lesssim|S^{\star}|.$

Putting all the bounds together gives

$\mathsf{E}_{\theta^{\star}}\int\|\theta-\theta^{\star}\|^{2}\,\Delta^{n}(d\theta)\lesssim|S^{\star}|+\varepsilon_{n}^{2}(\theta^{\star})+1\lesssim\varepsilon_{n}^{2}(\theta^{\star}).\qed$

###  A.2 Proofs of Theorems 1–3

###### Proof of Theorem 1.

By Markov’s inequality,

$\Delta^{n}(\\{\theta:\|\theta-\theta^{\star}\|^{2}>M_{n}\varepsilon_{n}^{2}(\theta^{\star})\\})\leq\frac{1}{M_{n}\varepsilon_{n}^{2}(\theta^{\star})}\int\|\theta-\theta^{\star}\|^{2}\,\Delta^{n}(d\theta).$

Taking expectation and applying the bound in Lemma 2 gives

$\mathsf{E}_{\theta^{\star}}\Delta^{n}(\\{\theta:\|\theta-\theta^{\star}\|^{2}>M_{n}\varepsilon_{n}^{2}(\theta^{\star})\\})\lesssim
M_{n}^{-1},$

and since $M_{n}\to\infty$, the claim follows. ∎

###### Proof of Theorem 2.

By Jensen’s inequality,

$\|\hat{\theta}-\theta^{\star}\|^{2}\leq\int\|\theta-\theta^{\star}\|^{2}\,\Delta^{n}(d\theta).$

Then the claim follows by Lemma 2. ∎

###### Proof of Theorem 3.

By Markov’s inequality

$\Delta^{n}(\\{\theta:|S_{\theta}|>M_{n}|S_{\theta^{\star}}|\\})\leq\frac{1}{M_{n}|S_{\theta^{\star}}|}\sum_{i=1}^{n}\phi_{i},$

where the sum on the right-hand side is the expectation of $|S_{\theta}|$
under $\theta\sim\Delta^{n}$. Take expectation of both sides and apply the
bound in Lemma 1, to get

$\mathsf{E}_{\theta^{\star}}\Delta^{n}(\\{\theta:|S_{\theta}|>M_{n}|S_{\theta^{\star}}|\\})\lesssim
M_{n}^{-1},$

and, since $M_{n}\to\infty$, the claim follows. ∎

###  A.3 Proof of Theorem 4

Let $\delta^{n}$ denote the mass function of the marginal distribution of
$S_{\theta}$ under $\theta\sim\Delta^{n}$, i.e.,

$\delta^{n}(S)=\Delta^{n}(\\{\theta:S_{\theta}=S\\}),\quad
S\subseteq\\{1,2,\ldots,n\\}.$

Also, for the given $\theta^{\star}$, let $S^{\star}=S_{\theta^{\star}}$. From
the simple form of $\Delta^{n}$, it’s easy to check that

$\delta^{n}(S)=\prod_{i\in S}\phi_{i}\cdot\prod_{i\not\in S}(1-\phi_{i}).$

This leads to a convenient bound

$\delta^{n}(S)\leq\frac{\delta^{n}(S)}{\delta^{n}(S^{\star})}=\prod_{i\in
S\cap S^{\star c}}e^{\operatorname{logit}(\phi_{i})}\,\prod_{i\in S^{c}\cap
S^{\star}}e^{-\operatorname{logit}(\phi_{i})}.$

Since each $\phi_{i}$ only depends on $Y_{i}$, and these are independent, we
can interchange the order of expectation and product. Also, for those $i\in
S^{\star c}$, with $\theta_{i}^{\star}=0$, the $\phi_{i}$’s are iid, so each
term in that product has the same expectation. Therefore,

$\mathsf{E}_{\theta^{\star}}\delta^{n}(S)\leq\bigl{\\{}\mathsf{E}_{0}e^{\operatorname{logit}(\phi_{1})}\bigr{\\}}^{|S\cap
S^{\star c}|}\,\prod_{i\in S^{c}\cap
S^{\star}}\mathsf{E}_{\theta_{i}^{\star}}e^{-\operatorname{logit}(\phi_{i})}.$

By the moment-generating function bounds in (11) and (12),

$\displaystyle\mathsf{E}_{0}e^{\operatorname{logit}(\phi_{1})}$
$\displaystyle\leq c_{0}\xi_{n}$
$\displaystyle\mathsf{E}_{\theta_{i}^{\star}}e^{-\operatorname{logit}(\phi_{i})}$
$\displaystyle\leq c_{1}\xi_{n}^{-1}e^{-k\theta_{i}^{\star 2}},$

where $k=\alpha/2\sigma^{2}$,
$\xi_{n}=\exp\\{\operatorname{logit}(\lambda_{n})\\}\sim n^{-(1+a)}$, and
$c_{0}$ and $c_{1}$ are the hidden proportionality constants in (11) and (12),
respectively.

Note also that, by Theorem 3, the $\delta^{n}$-probability of the event
“$|S|>M_{n}|S^{\star}|$” has vanishing expectation for any $M_{n}\to\infty$.
Thanks to this, I can immediately restrict my attention to those $S$ such that
“$|S|\leq M_{n}|S^{\star}|$” in what follows.

Now consider two distinct cases separately, namely, $S\supset S^{\star}$ and
$S\not\supseteq S^{\star}$. First, for any $S\supset S^{\star}$, it follows
that $|S^{c}\cap S^{\star}|=0$. So

$\displaystyle\mathsf{E}_{\theta^{\star}}\delta^{n}(S:S\supset S^{\star})$
$\displaystyle\leq\sum_{S:S\supset S^{\star},|S|\leq
M_{n}|S^{\star}|}\\{\mathsf{E}_{0}e^{\operatorname{logit}(\phi_{1})}\\}^{|S\cap
S^{\star c}|}$
$\displaystyle=\sum_{t=1}^{M_{n}|S^{\star}|}\binom{n-|S^{\star}|}{t}\\{\mathsf{E}_{0}e^{\operatorname{logit}(\phi_{1})}\\}^{t}$
$\displaystyle\leq\sum_{t=1}^{M_{n}|S^{\star}|}\\{e(n-|S^{\star}|)\mathsf{E}_{0}e^{\operatorname{logit}(\phi_{1})}\\}^{t}$
$\displaystyle\lesssim(n-|S^{\star}|)\xi_{n}.$

By definition of $\xi_{n}$, the upper bound is vanishing as $n\to\infty$,
proving the first claim.

Next, for any $S\not\supseteq S^{\star}$, there must be at least one component
in $S^{\star}$ that is not included in $S$. So,

$\displaystyle\mathsf{E}_{\theta^{\star}}\delta^{n}(S:S\not\supseteq
S^{\star})$ $\displaystyle\leq\sum_{S:S\not\supseteq S^{\star},|S|\leq
M_{n}|S^{\star}|}\Bigl{[}\bigl{\\{}\mathsf{E}_{0}e^{\operatorname{logit}(\phi_{1})}\bigr{\\}}^{|S\cap
S^{\star c}|}\,\prod_{i\in S^{c}\cap
S^{\star}}\mathsf{E}_{\theta_{i}^{\star}}e^{-\operatorname{logit}(\phi_{i})}\Bigr{]}$
$\displaystyle\leq\sum_{S:S\not\supseteq S^{\star},|S|\leq
M_{n}|S^{\star}|}(c_{0}\xi_{n})^{|S\cap S^{\star
c}|}(c_{1}\xi_{n}^{-1}e^{-kH^{2}})^{|S^{c}\cap S^{\star}|}$
$\displaystyle=\sum_{s=0}^{M_{n}|S^{\star}|}\sum_{t=0}^{s\wedge(|S^{\star}|-1)}\binom{|S^{\star}|}{t}\binom{n-|S^{\star}|}{s-t}(c_{0}\xi_{n})^{s-t}(c_{1}\xi_{n}^{-1}e^{-kH^{2}})^{|S^{\star}|-t}$
$\displaystyle\leq\sum_{s=0}^{M_{n}|S^{\star}|}\sum_{t=0}^{s\wedge(|S^{\star}|-1)}\\{c_{0}(n-|S^{\star}|)\xi_{n}\\}^{s-t}\\{c_{1}|S^{\star}|\xi_{n}^{-1}e^{-kH^{2}}\\}^{|S^{\star}|-t}.$

(In the above derivation, $s$ represents $|S|$ and $t$ represents $|S\cap
S^{\star}|$, which implies $s-t=|S\cap S^{\star c}|$ and
$|S^{\star}|-t=|S^{c}\cap S^{\star}|$.) Note that $t<|S^{\star}|$ because
$S\not\supseteq S^{\star}$ implies that $S$ can’t include all the entries in
$S^{\star}$. This means that there is a constant factor

$|S^{\star}|\xi_{n}^{-1}e^{-kH^{2}},$ (15)

which goes to 0 as $n\to\infty$ if $H$ is as in (9). The terms that involve
$(n-|S^{\star}|)\xi_{n}$ also vanish as in the previous case above. So all the
terms in the sum are geometrically small, hence the sum is bounded. But since
the common factor (15) vanishes, the upper bound itself vanishes, proving the
second claim of the theorem.

###  A.4 Proof of Theorem 5

For any data-dependent radius $\hat{\rho}$,
“$B_{n}(\hat{\rho})\not\ni\theta^{\star}$” is equivalent to
“$\|\hat{\theta}-\theta^{\star}\|>\hat{\rho}$,” and the following
decomposition, which holds for any deterministic $R>0$, is helpful:

$\displaystyle\mathsf{P}_{\theta^{\star}}\\{\|\hat{\theta}-\theta^{\star}\|>\hat{\rho}\\}$
$\displaystyle\leq\mathsf{P}_{\theta^{\star}}\\{\|\hat{\theta}-\theta^{\star}\|>\hat{\rho},\,\hat{\rho}>R\\}+\mathsf{P}_{\theta^{\star}}\\{\hat{\rho}\leq
R\\}$
$\displaystyle\leq\mathsf{P}_{\theta^{\star}}\\{\|\hat{\theta}-\theta^{\star}\|^{2}>R^{2}\\}+\mathsf{P}_{\theta^{\star}}\\{\hat{\rho}^{2}\leq
R^{2}\\}.$ (16)

The first term in (16) has nothing to do with the radius, so this can be
approached the same way for both types of credible balls. Indeed, by Markov’s
inequality,

$\mathsf{P}_{\theta^{\star}}\\{\|\hat{\theta}-\theta^{\star}\|^{2}>R^{2}\\}\leq
R^{-2}\mathsf{E}_{\theta^{\star}}\|\hat{\theta}-\theta^{\star}\|^{2}.$

By Theorem 2, the expectation in the upper bound is no more than
$M^{\prime}\varepsilon_{n}^{2}(\theta^{\star})$ for some constant
$M^{\prime}>0$ and $\varepsilon_{n}^{2}(\theta^{\star})$ in (7). Therefore, if
$R^{2}$ is a suitable multiple of $\varepsilon_{n}^{2}(\theta^{\star})$, then
the first term in (16) can be made less than a fraction of $\zeta$. The
specific constants depend on how $\hat{r}$ is defined, and the details for
each case are presented below.

###### Proof for the quantile-based radius.

Start with a bound on the non-coverage probability for the quantile-based
radius. Here I’ll bound the non-coverage probability by a sum of three terms,
each will be bounded by $\zeta/3$, for sufficiently large $n$. The first of
these three terms comes from the above analysis, so we set
$R^{2}=(3M^{\prime}/\zeta)\varepsilon_{n}^{2}(\theta^{\star})$ and conclude

$\mathsf{P}_{\theta^{\star}}\\{\|\hat{\theta}-\theta^{\star}\|^{2}>R^{2}\\}\leq\zeta/3.$

Next, to bound the second term in (16), recall that
$\hat{\rho}^{2}=M^{2}g_{n}^{2}\hat{r}^{2}$, where $\hat{r}$ is based on the
$(1-\zeta)$-quantile of $\Delta^{n}$ and $M$ is some sufficiently large
constant yet to be determined. Define $\widetilde{R}=(Mg_{n})^{-1/2}R$. Then
by definition of $\hat{r}$, and Markov’s inequality (again),

$\displaystyle\mathsf{P}_{\theta^{\star}}\\{\hat{\rho}\leq R\\}$
$\displaystyle=\mathsf{P}_{\theta^{\star}}\\{\hat{r}\leq\widetilde{R}\\}$
$\displaystyle\leq\mathsf{P}_{\theta^{\star}}\\{\Delta^{n}(\|\theta-\hat{\theta}\|\leq\widetilde{R})\geq
1-\zeta\\}$
$\displaystyle\leq(1-\zeta)^{-1}\mathsf{E}_{\theta^{\star}}\Delta^{n}(\|\theta-\hat{\theta}\|\leq\widetilde{R}).$

So the second term in (16) can be upper-bounded if the above expectation can
be upper-bounded. Towards this, let $S^{\star}=S_{\theta^{\star}}$ and use the
total probability formula to write

$\displaystyle\Delta^{n}(\|\theta-\hat{\theta}\|\leq\widetilde{R})$
$\displaystyle=\sum_{S}\Delta^{n}(\|\theta-\hat{\theta}\|\leq\widetilde{R}\mid
S)\,\delta^{n}(S)$
$\displaystyle\leq\\{1-\delta^{n}(S^{\star})\\}+\Delta^{n}(\|\theta-\hat{\theta}\|\leq\widetilde{R}\mid
S^{\star}).$

Under the conditions of Theorem 4, the first term has expectation that
vanishes with $n$, hence is eventually smaller than $\zeta/3$. So it remains
to look at the second term in the above display. The conditional distribution
of $\theta$ under $\Delta^{n}$, given $S^{\star}$, is

$\theta\sim{\sf
N}_{|S^{\star}|}(Y_{S^{\star}},\tau^{2}I)\otimes\delta_{|S^{\star c}|},$

where $\tau^{2}=\sigma^{2}(\alpha+\gamma)^{-1}$ from (4) is deterministic. For
such a $\theta$,

$\displaystyle\|\theta-\hat{\theta}\|^{2}$
$\displaystyle\overset{d}{=}\sum_{i\in S^{\star}}\\{\tau
G_{i}+(1-\phi_{i})Y_{i}\\}^{2}+\sum_{i\not\in S^{\star}}(\phi_{i}Y_{i})^{2}$
$\displaystyle\geq\sum_{i\in S^{\star}}\\{\tau
G_{i}+(1-\phi_{i})Y_{i}\\}^{2},$

where $G_{S^{\star}}=\\{G_{i}:i\in S^{\star}\\}$ are iid ${\sf N}(0,1)$. Then

$\displaystyle\Delta^{n}(\|\theta-\hat{\theta}\|\leq\widetilde{R}\mid
S^{\star})$ $\displaystyle\leq\mathsf{P}\Bigl{\\{}\sum_{i\in S^{\star}}\\{\tau
G_{i}+(1-\phi_{i})Y_{i}\\}^{2}\leq\widetilde{R}^{2}\Bigr{\\}}$
$\displaystyle\leq\mathsf{P}\\{\tau^{2}\|G_{S^{\star}}\|^{2}\leq\widetilde{R}^{2}\\},$

where the second line follows by Anderson’s inequality. Note that
$\|G_{S^{\star}}\|^{2}\sim{\sf ChiSq}(|S^{\star}|)$, which means that
$\|G_{S^{\star}}\|^{2}$ scales like $|S^{\star}$. Consider two cases:
$|S^{\star}|=O(1)$ and $|S^{\star}|\to\infty$. In the former case,
$\|G_{S^{\star}}\|^{2}$ is stochastically bounded, so the upper bound in the
above display can be made less than $\zeta/3$ if $\widetilde{R}$ is small or,
equivalently, if $M$ is sufficiently large. For the latter case, with
$|S^{\star}|\to\infty$, the above probability can be bounded as

$\displaystyle\|G_{S^{\star}}\|^{2}\leq\tau^{-2}\widetilde{R}^{2}$
$\displaystyle\iff\|G_{S^{\star}}\|^{2}-|S^{\star}|\leq\tau^{-2}\widetilde{R}^{2}-|S^{\star}|$
$\displaystyle\iff\|G_{S^{\star}}\|^{2}-|S^{\star}|\leq-w|S^{\star}|,$

where $w=1-(\tau|S^{\star}|)^{-1}\widetilde{R}^{2}$. To see that $w>0$ or,
equivalently, that $\widetilde{R}^{2}<\tau|S^{\star}|$ for sufficiently large
$n$, plug in the definition of $\widetilde{R}$ to get

$\frac{(3M^{\prime}/\zeta)\varepsilon_{n}^{2}(\theta^{\star})}{Mg_{n}}<\tau|S^{\star}|\iff\frac{\log(en/|S^{\star}|)}{\log(en)}<\frac{M\tau}{3M^{\prime}/\zeta}.$

Note that the right-most inequality holds if $M>3M^{\prime}/(\tau\zeta)$.
Therefore,

$\mathsf{P}\\{\tau^{2}\|G_{S^{\star}}\|^{2}\leq\widetilde{R}^{2}\\}\leq\mathsf{P}\\{\|G_{S^{\star}}\|^{2}-|S^{\star}|\leq-w|S^{\star}|\\}.$

Using a standard tail probability bound for chi-square random variables (e.g.,
Laurent and Massart, 2000, Lemma 1), the probability in the above display is
bounded by

$\mathsf{P}\\{\tau\|G_{S}\|^{2}\leq\widetilde{R}^{2}\\}\leq\mathsf{P}\\{\|G_{S^{\star}}\|^{2}-|S^{\star}|\leq-2(|S^{\star}|x)^{1/2}\\}\leq
e^{-x},$

where $x=w^{2}|S^{\star}|/4$. Since $|S^{\star}|\to\infty$, this upper bound
will eventually be less than $\zeta/3$ for $M$ as above. Putting everything
together, all three terms upper bounding the non-coverage probability are less
than $\zeta/3$ for $n$ sufficiently large.

Now, for the size of the credible ball. Write
$\varepsilon_{n}^{2}=\varepsilon_{n}^{2}(\theta^{\star})$. By definition of
$\hat{r}$,

$\hat{r}^{2}>L\varepsilon_{n}^{2}\iff\Delta^{n}(\|\theta-\hat{\theta}\|\leq
L\varepsilon_{n}^{2})<1-\zeta\iff\Delta^{n}(\|\theta-\hat{\theta}\|>L\varepsilon_{n}^{2})>\zeta.$

For the constant $M^{\prime}$ used above from Theorem 2, and for the specified
threshold $\eta>0$, take $L\geq M^{\prime}\\{1+(\zeta\eta/2)^{-1/2}\\}^{2}$.
Then the triangle inequality implies

$\Delta^{n}(\|\theta-\hat{\theta}\|^{2}>L\varepsilon_{n}^{2})\leq\Delta^{n}(\|\theta-\theta^{\star}\|^{2}>M^{\prime}\varepsilon_{n}^{2})+1\\{\|\hat{\theta}-\theta^{\star}\|^{2}>(\zeta\eta/2)^{-1}M^{\prime}\varepsilon_{n}^{2}\\}.$

Using Markov’s inequality twice gives

$\displaystyle\mathsf{P}_{\theta^{\star}}\\{\hat{r}^{2}>L\varepsilon_{n}^{2}\\}$
$\displaystyle\leq\zeta^{-1}\mathsf{E}_{\theta^{\star}}\Delta^{n}(\|\theta-\hat{\theta}\|>L\varepsilon_{n}^{2})$
$\displaystyle\leq\zeta^{-1}\Bigl{\\{}\mathsf{E}_{\theta^{\star}}\Delta^{n}(\|\theta-\theta^{\star}\|^{2}>M^{\prime}\varepsilon_{n}^{2})+\frac{\zeta\eta\mathsf{E}_{\theta^{\star}}\|\hat{\theta}-\theta^{\star}\|^{2}}{2M^{\prime}\varepsilon_{n}^{2}}\Bigr{\\}}.$

The first term in the curly brackets is vanishing as $n\to\infty$ by Theorem 1
and, hence, will eventually be less than $\zeta\eta/2$. By Theorem 2 and the
definition of $M^{\prime}$, the second term in the curly brackets is no more
than $\zeta\eta/2$. Therefore, the upper bound is no more than $\eta$, which
proves the claim. ∎

###### Proof for the plug-in estimator-based radius.

Following the same argument as above, let
$R^{2}=(2M^{\prime}/\zeta)\varepsilon_{n}^{2}(\theta^{\star})$. Then Theorem 2
implies that the first term in (16) is no more than $\zeta/2$. For the second
term in (16), note that

$\text{$x\mapsto x\log(en/x)$ is increasing on $[0,n]$}.$ (17)

This implies that, for sufficiently small $c>0$,

$\displaystyle|\hat{S}|\log\tfrac{en}{|\hat{S}|}\leq
c|S^{\star}|\log\tfrac{en}{|S^{\star}|}$
$\displaystyle\implies|\hat{S}|\log\tfrac{en}{|\hat{S}|}\leq
c|S^{\star}|\log\tfrac{en}{c|S^{\star}|}$ $\displaystyle\iff|\hat{S}|\leq
c|S^{\star}|.$

Therefore, to get a bound on $\mathsf{P}_{\theta^{\star}}\\{\hat{\rho}^{2}\leq
R^{2}\\}$, with $\hat{\rho}=M\hat{r}$, it suffices to bound

$\mathsf{P}_{\theta^{\star}}\\{|\hat{S}|\leq c|S^{\star}|\\},$

where $c=2M^{\prime}/M\zeta$, which can be made small with suitable choice of
$M$. Note that $|\hat{S}|$ is a sum of independent but non-identically
distributed Bernoulli random variables, so its expected value is

$\mu^{\star}=\mathsf{E}_{\theta^{\star}}|\hat{S}|=\sum_{i\in
S^{\star}}\mathsf{P}_{\theta_{i}^{\star}}(\phi_{i}>\tfrac{1}{2})+(n-|S^{\star}|)\mathsf{P}_{0}(\phi_{i}>\tfrac{1}{2}).$

Both lower and upper bounds for $\mu^{\star}$ are needed. Towards this,

$\mathsf{P}_{0}(\phi_{i}>\tfrac{1}{2})=\mathsf{P}_{0}\\{\text{logit}(\phi_{i})>0\\}=\mathsf{P}\bigl{[}(Z/\sigma)^{2}>\tfrac{2\sigma^{2}}{\alpha}\\{\operatorname{logit}(\lambda_{n})+\tfrac{1}{2}\log\tfrac{\gamma}{\alpha+\gamma}\\}\bigr{]}.$

From the tail probability bound (13), it follows that

$\mathsf{P}_{0}(\phi_{i}>\tfrac{1}{2})\lesssim n^{-(1+a)},\quad i\not\in
S^{\star},$

and, therefore,

$\mu^{\star}\leq|S^{\star}|+n^{-(1+a)}(n-|S^{\star}|)=\\{1+o(n^{-a})\\}|S^{\star}|.$
(18)

For the lower bound on $\mu^{\star}$,

$\mu^{\star}\geq\sum_{i\in
S^{\star}}\mathsf{P}_{\theta_{i}^{\star}}(\phi_{i}>\tfrac{1}{2})\geq|S^{\star}|\mathsf{P}_{H}(\phi_{i}>\tfrac{1}{2}),$

where $H$ is the minimum (non-zero) signal size in (9). Then

$\displaystyle\mathsf{P}_{H}(\phi_{i}>\tfrac{1}{2})$
$\displaystyle=1-\mathsf{P}_{H}\\{\operatorname{logit}(\phi_{i})<0\\}$
$\displaystyle=1-\mathsf{P}_{H}\bigl{[}-\tfrac{\alpha}{2\sigma^{2}}Y_{i}^{2}>-\\{\operatorname{logit}(\lambda_{n})+\tfrac{1}{2}\log\tfrac{\gamma}{\alpha+\gamma}\\}\bigr{]}.$

Apply the exponential function to both sides of the inequality inside
$\mathsf{P}_{H}(\cdots)$ and then use Markov’s inequality, the bound in (12),
the size of $H$ in (9), and an argument like that at the end of the proof of
Theorem 4 to get

$\mathsf{P}_{H}\bigl{[}-\tfrac{\alpha}{2\sigma^{2}}Y_{i}^{2}>-\\{\operatorname{logit}(\lambda_{n})+\tfrac{1}{2}\log\tfrac{\gamma}{\alpha+\gamma}\\}\bigr{]}=o(1),\quad
n\to\infty.$

Therefore, $\mu^{\star}\geq\\{1-o(1)\\}|S^{\star}|$. Using the standard
Chernoff bounds for Bernoulli random variables, we get

$\mathsf{P}_{\theta^{\star}}\\{|\hat{S}|\leq
c|S^{\star}|\\}=\inf_{t>0}e^{t|S^{\star}|}\mathsf{E}_{\theta^{\star}}e^{-t|\hat{S}|}=\exp\Bigl{\\{}c|S^{\star}|\log\frac{\mu^{\star}}{c|S^{\star}|}+c|S^{\star}|-\mu^{\star}\Bigl{\\}}.$

Plug in the lower and upper bounds for $\mu^{\star}$ to get

$\mathsf{P}_{\theta^{\star}}\\{|\hat{S}|\leq
c|S^{\star}|\\}\leq\exp\Bigl{[}|S^{\star}|\Bigl{\\{}c\log\frac{1+o(n^{-a})}{c}-(1-c)+o(1)\Bigr{\\}}\Bigr{]}.$

The term inside $\\{\cdots\\}$ is negative for $c<1$ and $n$ sufficiently
large, so the right-hand side can be upper-bounded by $\zeta/2$ if $M$ is a
sufficiently large multiple of $2M^{\prime}/\zeta$.

Finally, for the size of the credible ball, by the same monotonicity property
(17) used above, it follows that

$\mathsf{P}_{\theta^{\star}}\\{\hat{r}>L\varepsilon_{n}^{2}\\}\leq\mathsf{P}_{\theta^{\star}}\\{|\hat{S}|>L|S^{\star}|\\}.$

By Markov’s inequality and the upper bound on
$\mu^{\star}=\mathsf{E}_{\theta^{\star}}|\hat{S}|$ in (18),

$\mathsf{P}_{\theta^{\star}}\\{|\hat{S}|>L|S^{\star}|\\}\leq\frac{\\{1+o(n^{-a})\\}|S^{\star}|}{L|S^{\star}|}.$

Therefore, if $L>\eta^{-1}$, then the size condition
$\mathsf{P}_{\theta^{\star}}\\{\hat{r}^{2}>L\varepsilon_{n}^{2}\\}\leq\eta$
holds. ∎

## References

  * Abramovich et al., (2006) Abramovich, F., Benjamini, Y., Donoho, D. L., and Johnstone, I. M. (2006).  Adapting to unknown sparsity by controlling the false discovery rate.  Ann. Statist., 34(2):584–653.
  * Belitser, (2017) Belitser, E. (2017).  On coverage and local radial rates of credible sets.  Ann. Statist., 45(3):1124–1151.
  * Belitser and Ghosal, (2019) Belitser, E. and Ghosal, S. (2019).  Empirical Bayes oracle uncertainty quantification.  Ann. Statist., to appear, http://www4.stat.ncsu.edu/~ghoshal/papers/oracle_regression.pdf.
  * Belitser and Nurushev, (2020) Belitser, E. and Nurushev, N. (2020).  Needles and straw in a haystack: Robust confidence for possibly sparse sequences.  Bernoulli, 26(1):191–225.
  * Bhadra et al., (2017) Bhadra, A., Datta, J., Polson, N. G., and Willard, B. (2017).  The horseshoe+ estimator of ultra-sparse signals.  Bayesian Anal., 12(4):1105–1131.
  * Bhattacharya et al., (2015) Bhattacharya, A., Pati, D., Pillai, N. S., and Dunson, D. B. (2015).  Dirichlet-Laplace priors for optimal shrinkage.  J. Amer. Statist. Assoc., 110(512):1479–1490.
  * Boucheron et al., (2013) Boucheron, S., Lugosi, G., and Massart, P. (2013).  Concentration Inequalities.  Oxford University Press, Oxford.  A Nonasymptotic Theory of Independence, With a foreword by Michel Ledoux.
  * Carvalho et al., (2010) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2010).  The horseshoe estimator for sparse signals.  Biometrika, 97(2):465–480.
  * Castillo and Szabó, (2020) Castillo, I. and Szabó, B. (2020).  Spike and slab empirical Bayes sparse credible sets.  Bernoulli, 26(1):127–158.
  * Castillo and van der Vaart, (2012) Castillo, I. and van der Vaart, A. (2012).  Needles and straw in a haystack: posterior concentration for possibly sparse sequences.  Ann. Statist., 40(4):2069–2101.
  * Donoho and Johnstone, (1994) Donoho, D. L. and Johnstone, I. M. (1994).  Minimax risk over $l_{p}$-balls for $l_{q}$-error.  Probab. Theory Related Fields, 99(2):277–303.
  * Donoho et al., (1992) Donoho, D. L., Johnstone, I. M., Hoch, J. C., and Stern, A. S. (1992).  Maximum entropy and the nearly black object.  J. Roy. Statist. Soc. Ser. B, 54(1):41–81.  With discussion and a reply by the authors.
  * Efron, (2010) Efron, B. (2010).  Large-Scale Inference, volume 1 of Institute of Mathematical Statistics Monographs.  Cambridge University Press, Cambridge.
  * Efron and Morris, (1973) Efron, B. and Morris, C. (1973).  Stein’s estimation rule and its competitors—an empirical Bayes approach.  J. Amer. Statist. Assoc., 68:117–130.
  * Honorio and Jaakkola, (2014) Honorio, J. and Jaakkola, T. (2014).  Tight Bounds for the Expected Risk of Linear Classifiers and PAC-Bayes Finite-Sample Guarantees.  In Kaski, S. and Corander, J., editors, Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, volume 33 of Proceedings of Machine Learning Research, pages 384–392, Reykjavik, Iceland. PMLR.
  * James and Stein, (1961) James, W. and Stein, C. (1961).  Estimation with quadratic loss.  In Proc. 4th Berkeley Sympos. Math. Statist. and Prob., Vol. I, pages 361–379. Univ. California Press, Berkeley, Calif.
  * Johnstone and Silverman, (2004) Johnstone, I. M. and Silverman, B. W. (2004).  Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences.  Ann. Statist., 32(4):1594–1649.
  * Laurent and Massart, (2000) Laurent, B. and Massart, P. (2000).  Adaptive estimation of a quadratic functional by model selection.  Ann. Statist., 28(5):1302–1338.
  * Liu et al., (2020) Liu, C., Martin, R., and Shen, W. (2020).  Empirical priors and posterior concentration in a piecewise polynomial sequence model.  arXiv:1712.03848.
  * Martin and Liu, (2013) Martin, R. and Liu, C. (2013).  Inferential models: a framework for prior-free posterior probabilistic inference.  J. Amer. Statist. Assoc., 108(501):301–313.
  * Martin and Liu, (2015) Martin, R. and Liu, C. (2015).  Inferential Models: Reasoning with Uncertainty, volume 147 of Monographs on Statistics and Applied Probability.  CRC Press, Boca Raton, FL.
  * Martin et al., (2017) Martin, R., Mess, R., and Walker, S. G. (2017).  Empirical Bayes posterior concentration in sparse high-dimensional linear models.  Bernoulli, 23(3):1822–1847.
  * Martin and Ning, (2020) Martin, R. and Ning, B. (2020).  Empirical priors and coverage of posterior credible sets in a sparse normal mean model.  Sankhyā A., 82:477–498.  Special issue in memory of Jayanta K. Ghosh.
  * Martin and Tang, (2020) Martin, R. and Tang, Y. (2020).  Empirical priors for prediction in sparse high-dimensional linear regression.  J. Mach. Learn. Res., 21(144):1–30.
  * Martin and Walker, (2014) Martin, R. and Walker, S. G. (2014).  Asymptotically minimax empirical Bayes estimation of a sparse normal mean vector.  Electron. J. Stat., 8(2):2188–2206.
  * Martin and Walker, (2019) Martin, R. and Walker, S. G. (2019).  Data-dependent priors and their posterior concentration rates.  Electron. J. Stat., 13(2):3049–3081.
  * Ray and Szabó, (2020) Ray, K. and Szabó, B. (2020).  Variational Bayes for high-dimensional linear regression with sparse priors.  J. Amer. Statist. Assoc., to appear; arXiv:1904.07150.
  * Rousseau and Szabo, (2020) Rousseau, J. and Szabo, B. (2020).  Asymptotic frequentist coverage properties of Bayesian credible sets for sieve priors.  Ann. Statist., 48(4):2155–2179.
  * Ročková and George, (2018) Ročková, V. and George, E. I. (2018).  The spike-and-slab LASSO.  J. Amer. Statist. Assoc., 113(521):431–444.
  * Stein, (1956) Stein, C. (1956).  Inadmissibility of the usual estimator for the mean of a multivariate normal distribution.  In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, 1954–1955, vol. I, pages 197–206, Berkeley and Los Angeles. University of California Press.
  * Stein, (1981) Stein, C. M. (1981).  Estimation of the mean of a multivariate normal distribution.  Ann. Statist., 9(6):1135–1151.
  * Szabó et al., (2015) Szabó, B., van der Vaart, A. W., and van Zanten, J. H. (2015).  Frequentist coverage of adaptive nonparametric Bayesian credible sets.  Ann. Statist., 43(4):1391–1428.
  * van der Pas and Rockova, (2017) van der Pas, S. and Rockova, V. (2017).  Bayesian dyadic trees and histograms for regression.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems 30, pages 2089–2099. Curran Associates, Inc.
  * van der Pas et al., (2016) van der Pas, S., Scott, J., Chakraborty, A., and Bhattacharya, A. (2016).  horseshoe: Implementation of the Horseshoe Prior.  R package version 0.1.0.
  * van der Pas et al., (2017) van der Pas, S., Szabó, B., and van der Vaart, A. (2017).  Uncertainty quantification for the horseshoe (with discussion).  Bayesian Anal., 12(4):1221–1274.  With a rejoinder by the authors.
  * van Erven and Szabó, (2020) van Erven, T. and Szabó, B. (2020).  Fast exact Bayesian inference for sparse signals in the normal sequence model.  Bayesian Anal., to appear; arXiv:1810.10883.
  * Yang and Martin, (2020) Yang, Y. and Martin, R. (2020).  Empirical priors and variational approximations of the posterior in high-dimensional linear models.  arXiv:2007.15930.