# Coverage of Credible Sets for Regression under Variable Selection

Samhita Pal, Subhashis Ghosal  Samhita Pal
Department of Statistics, North Carolina State University
and
Subhashis Ghosal
Department of Statistics, North Carolina State University<EMAIL_ADDRESS>
###### Abstract

We study the asymptotic frequentist coverage of credible sets based on a novel
Bayesian approach for a multiple linear regression model under variable
selection. We initially ignore the issue of variable selection, which allows
us to put a conjugate normal prior on the coefficient vector. The variable
selection step is incorporated directly in the posterior through a sparsity-
inducing map and uses the induced prior for making an inference instead of the
natural conjugate posterior. The sparsity-inducing map minimizes the sum of
the squared $\ell_{2}$-distance weighted by the data matrix and a suitably
scaled $\ell_{1}$-penalty term. We obtain the limiting coverage of various
credible regions and demonstrate that a modified credible interval for a
component has the exact asymptotic frequentist coverage if the corresponding
predictor is asymptotically uncorrelated with other predictors. Through
extensive simulation, we provide a guideline for choosing the penalty
parameter as a function of the credibility level appropriate for the
corresponding coverage. We also show finite-sample numerical results that
support the conclusions from the asymptotic theory. We also provide the
credInt package that implements the method in R to obtain the credible
intervals along with the posterior samples.

##  1 Introduction

The multiple linear regression model is one of the most useful tools for
analyzing data, expressing an overall relationship between a response variable
and a set of predictors as an affine function. Not all listed predictor
variables are often active in the regression, prompting the need to select the
relevant predictors for a better interpretable model and more precise
estimation and prediction. The issue of variable selection is especially vital
for data consisting of many predictors when a meaningful inference is possible
only by incorporating such a step in the analysis, but even for a fixed-
dimensional setting, a variable selection step is highly desirable. This leads
us to the model selection problem, which has been addressed by both Bayesian
and non-Bayesian methods. Classical methods of model selection include forward
selection and backward selection. The former begins with the null model and
keeps adding predictors in the model sequentially as long as the predictive
power goes up significantly, while the latter begins with the full model and
sequentially removes unneeded until a stopping criterion is met. Another
commonly employed approach involves penalization, which adds a penalty
function to the objective function. This encourages the minimizer to produce
solutions with greater sparsity than a minimizer without a penalty. The most
notable such method is the LASSO [Tibshirani, 1996], which penalizes the
$\ell_{1}$-norm of the coefficient vectors to obtain a sparse estimator.
Asymptotic properties such as consistency and the limiting distribution of the
LASSO were studied in Fu and Knight [2000] in the fixed dimensional setting.
If the dimension of the model is large, possibly larger than the sample size,
additional compatibility conditions are needed on the predictors and the true
value of the vector of regression coefficients to ensure identifiability. For
a detailed account of convergence and selection properties of LASSO-type
estimators, see Bühlmann and van de Geer [2011].

In the Bayesian domain, the emphasis has been on formulating a prior
distribution that encourages sparsity in the posterior model. The approach of
a spike-and-slab prior [Mitchell and Beauchamp, 1988, Ishwaran and Rao, 2005]
combines a point mass or an approximation to it (spike) at zero and a thick-
tailed continuous distribution (slab). Gibbs sampling and a data augmentation
technique are used to compute the posterior distribution; see George and
McCulloch [1993]. A typically faster alternative is to adopt a continuous-
shrinkage prior such as the horseshoe prior [Carvalho et al., 2010]. It
incorporates a single density function to emulate the concentration at zero
and the thickness of the tail in a spike-and-slab prior, simplifying the model
by using a one-component approach. Asymptotic concentration and selection
properties of the posterior based on a spike-and-slab prior were established
in Castillo et al. [2015] and analogous properties for a continuous shrinkage
prior in Song and Liang [2017]. Variational inference for the variable
selection problem was recently considered by some authors as an approximate
Bayesian inference method providing significantly faster computation; see Ray
and Szabó [2022], Ormerod et al. [2017], Huang et al. [2016], Mukherjee and
Sen [2022], Yang et al. [2020] Concentration properties for estimation and
selection consistency properties were studied by Zhang and Gao [2020], Han and
Yang [2019], Bai et al. [2020]. Yang and Martin [2020] showed asymptotic
frequentist coverage of a credible ball assuming an orthogonal design matrix.

While penalization methods for variable selection like the LASSO are hugely
popular and possess desirable convergence and model selection properties, they
do not naturally quantify uncertainty. Even in the fixed dimensional case,
because of the bias and degeneracy of the LASSO, the limiting distribution
obtained in Fu and Knight [2000] cannot be used to construct a confidence
region as the limiting distribution is dependent on the sparsity structure of
the true regression coefficients. By applying a debiasing technique on the
LASSO, [Zhang and Zhang, 2014, van de Geer et al., 2014, Javanmard and
Montanari, 2018] constructed confidence bands, but at the expense of losing
sparsity. Modifying the idea of the bootstrap LASSO in Fu and Knight [2000],
Chatterjee and Lahiri [2011] showed that the resulting confidence region has
good coverage in the fixed-dimensional setting. Confidence regions for the
regression vector in the high-dimensional linear regression setting under
sparsity appear to have been studied only in Nickl and van de Geer [2013] and
Cai and Guo [2017]. Bayesian methods automatically provide a measurement of
uncertainty through the posterior distribution. However, the Bayesian way of
quantifying the uncertainty may not match a frequentist assessment because
these two use very different notions of randomness. For smooth fixed-
dimensional parametric families, the Bernstein-von Mises theorem ensures the
approximate agreement between Bayesian and frequentist measures of
uncertainty, but this may not carry over to non-parametric and other settings;
see Cox [1993]. In smoothing regimes, asymptotic coverage may be assured only
by modifications such as undersmoothing or inflating a credible ball; see
Knapik et al. [2011], Castillo and Nickl [2013], Szabó et al. [2015], Yoo and
Ghosal [2016], Ray [2017], Sniekers and van der Vaart [2015], and others. For
a high-dimensional setting, the problem of frequentist validation of Bayesian
uncertainty quantification was addressed only by van der Pas et al. [2017],
Belitser and Nurushev [2019], Belitser and Ghosal [2020, 2024].

A novel and convenient Bayesian approach to inference on a restricted subspace
through the induced distribution by a restriction-complying “immersion” map,
which can often be a projection, from a simple unrestricted posterior
typically obtained through conjugacy. Lin and Dunson [2014] and Chakraborty
and Ghosal [2021a, b, 2022] used the technique for inference on monotone
functions; Bhaumik and Ghosal [2015, 2017], Bhaumik et al. [2022] for
differential equation models, and Wang and Ghosal [2023a, c, b] for problems
on multivariate monotonicity. The corresponding “projection-posterior”, or a
more general “immersion-posterior”, can be easily computed by conjugate
posterior sampling, typically followed by a simple optimization step. For
differential equation models, a Bernstein-von Mises theorem holds, justifying
the uncertainty quantification of the projection-posterior from the
frequentist principle. For shape-restricted models, the asymptotic coverage is
slightly higher than the credibility level, the opposite of the phenomenon
observed in Cox [1993] for smoothing problems. Interestingly, the relationship
between the limiting coverage and credibility does not depend on the true
function. Hence, a target exact asymptotic coverage can be obtained by
starting from a predetermined lower level (Chakraborty and Ghosal [2021b,
2022], Wang and Ghosal [2023a, b]). In a linear regression model with variable
selection, the need to address sparsity in the prior makes posterior
distributions harder to compute and analyze theoretically. This is
particularly apparent because results on coverage of Bayesian credible regions
are largely absent from the literature. The problem can be alleviated by using
an immersion posterior to make an inference. In the absence of sparsity, a
multivariate normal prior is conjugate for the linear regression problem,
leading to an explicit expression for the posterior distribution. An immersion
map given by the minimizer of the sum of a weighted squared Euclidean distance
and an $\ell_{1}$-penalty induces sparsity in the solution and is appropriate
for the problem.

This paper studies the construction and frequentist coverage properties of
credible sets for the coefficient vector using the sparse projection-posterior
method under the fixed dimension setting. We derive a weak limit of the
immersion-posterior distribution with an interesting structural similarity
with the distributional limit of the LASSO derived in Fu and Knight [2000].
Using the weak limit, we derive an expression for the limiting coverage of a
Bayesian credible ball. When a predictor variable is asymptotically
uncorrelated with the remaining predictors, we simplify the expression for the
coverage and compare it with the credibility level. While the limiting
coverage is lower than the credibility level if the coefficient is non-zero,
we can identify a higher credibility level, depending on the tuning parameter,
such that the limiting coverage coincides with the targeted coverage. The
corresponding limiting coverage is higher if the parameter assumes zero value,
provided that the intended coverage level is sufficiently high. Apart from
reconciling the Bayesian and frequentist notions of uncertainty quantification
as in the Bernsterin-von Mises theorem for smooth fixed-dimensional parametric
families, this result provides an automatic tuning mechanism.

The rest of the paper is organized as follows. The setup and the methodology
are formally described in the next section. In Section 3, we derive the
limiting results for establishing asymptotic coverage of the credible balls,
stating the necessary assumptions and confirming their pragmatism. Section 4
shows finite-sample results of the proven asymptotic guarantees and provides a
table for guidelines for choosing the correct penalty parameter. Finally, in
Section 5, we discuss the contribution of this work and the future directions,
followed by the proofs in the appendix.

##  2 Setup and Methodology

We consider the problem of model selection and uncertainty quantification of
the chosen variables in the context of the linear regression model

$\displaystyle\bm{Y}=\bm{X}\bm{\theta}+\bm{\varepsilon},\quad\bm{\varepsilon}\sim\mathrm{N}_{n}(\bm{0},\sigma^{2}\bm{I}_{n}),$
(2.1)

where $\bm{Y}\in\mathbb{R}^{n}$ is the response vector observed for $n$
samples, $\bm{X}\in\mathbb{R}^{n\times p}$ is a deterministic design matrix
and $\bm{\theta}\in\mathbb{R}^{p}$ is the vector of regression coefficients.
Regarding the methodology, the deterministic or random nature of the predictor
variables is immaterial, but it matters in the study of asymptotic properties.
If the predictors are random, the conditions imposed on the predictor
variables will have to be satisfied with probability tending to one.

Throughout the paper, the number of explanatory variables $p$ in this model is
assumed to be fixed and smaller than the sample size $n$. A smaller number $s$
of these $p$ predictors are active, but $s$ and the active predictor indices
are unknown. Although the regression vector may be estimated at the parametric
rate $n^{-1/2}$, identifying the model consisting only of the active
predictors is desirable since it gives a more interpretable model. Let $S_{0}$
consist of the indices corresponding to active predictor variables in the
model. We assume that the data $\bm{Y}=(Y_{1},\ldots,Y_{n})^{\mathrm{T}}$ is
generated through the process

$\displaystyle
Y_{i}=\bm{X}_{i}^{\mathrm{T}}\bm{\theta}^{0}+\varepsilon_{i},\quad\varepsilon_{1},\ldots,\varepsilon_{n}\mbox{
are i.i.d. with mean 0 and variance }\sigma_{0}^{2},$ (2.2)

where $\bm{X}_{i}^{\mathrm{T}}$ stands for the $i$th row of $\bm{X}$ (i.e.,
the vector of covariates for the $i$th observation). Thus, the true data
generating process follows the linear model (2.1) with the true value of the
regression coefficients $\bm{\theta}^{0}$ and the true value of the error
standard deviation $\sigma_{0}$, except that the error distribution need not
be normal.

Throughout the paper, we will use the following notations. Probability under
the true model will be denoted by $\mathbb{P}_{0}$ and the corresponding
expectation by $\mathbb{E}_{0}$. The normalized cross-products matrix is given
by $\bm{C}_{n}=n^{-1}{\bm{X}^{\mathrm{T}}\bm{X}}$. For notational convenience
in theoretical analysis, we shall place the $s_{0}$ active predictors
corresponding to non-zero regression coefficients at the beginning of the
arrangement: let $\bm{X}_{(1)}$ stand for the submatrix of the data matrix
$\bm{X}$ formed by the active predictors and $\bm{X}_{(2)}$ stand for the
submatrix corresponding to the remaining $p-s_{0}$ variables with zero true
coefficients. Then, the matrix $\bm{C}_{n}$ can be split into block matrices

$\displaystyle\begin{bmatrix}\bm{C}_{n(11)}&\bm{C}_{n(12)}\\\
\bm{C}_{n(21)}&\bm{C}_{n(22)}\end{bmatrix}=\frac{1}{n}\begin{bmatrix}\bm{X}_{(1)}^{\mathrm{T}}\bm{X}_{(1)}&\bm{X}_{(1)}^{\mathrm{T}}\bm{X}_{(2)}\\\
\bm{X}_{(2)}^{\mathrm{T}}\bm{X}_{(1)}&\bm{X}_{(2)}^{\mathrm{T}}\bm{X}_{(2)}\end{bmatrix}.$
(2.3)

The $q$-norm of a vector $\bm{a}\in\mathbb{R}^{d}$ is given by
$\left\|\bm{a}\right\|_{q}=(\sum_{i=1}^{d}|a_{i}|^{q})^{1/q}$. For $q=2$, this
corresponds to the Euclidean norm, which is denoted by $\|\cdot\|$. For a
matrix $\bm{M}\in\mathbb{R}^{d_{1}\times d_{2}}$, we denote by $\bm{M}_{A}$
the submatrix of containing the columns corresponding to a given set
$A\subset\\{1,2,\dots,d_{2}\\}$. We denote the zero vector by $\bm{0}$ and the
identity matrix of order $r$ by $\bm{I}_{r}$. We use the $\rightsquigarrow$
notation to denote weak convergence. For a sequence $a_{n}$ and positive
sequence $b_{n}$, let $a_{n}=\mathcal{O}(b_{n})$ mean that $|a_{n}|\leq
Cb_{n}$ for some constant $C>0$. For a sample space $\mathfrak{X}$, let
$\mathfrak{M}(\mathfrak{X})$ stand for the space of all probability measures
on $\mathfrak{X}$. The distributional law of a random element $T$ will be
denoted by $\mathcal{L}(T)$. The symbol $\mathrm{P}$ refers to a probability
statement with respect to a generic probability distribution. Let $\delta_{0}$
denote the Dirac delta measure at the point $0$. We shall write $\Phi$ for the
standard normal cumulative distribution function and
$z_{\alpha}=\Phi^{-1}(1-\alpha)$ for the $(1-\alpha)$-quantile of it.

The proposed immersion-posterior methodology for the variable selection
problem differs from the classical Bayesian approach of putting a prior taking
sparsity into account and then updating to the posterior distribution that
puts weights on different models corresponding to different sets of selected
predictors. In the proposed approach, we put a conditionally conjugate normal
prior on the parameter vector given, disregarding the issue of variable
selection at first, and then a conjugate inverse gamma prior for $\sigma^{2}$:
Further, we put an inverse gamma prior for $\sigma^{2}$, that is,

$\displaystyle\bm{\theta}|\sigma\sim\mathrm{N}_{p}(\textbf{0},\sigma^{2}a_{n}^{-1}\bm{I}_{p}),\qquad\sigma^{-2}\sim\mathrm{Ga}(b_{1},b_{2})$
(2.4)

for some sequence $a_{n}>0$ and constants $b_{1},b_{2}>0$. The non-informative
choice $b_{1}=b_{2}=0$, which corresponds to the density $\sigma^{-1}$ for
$\sigma$, may also be used. The choices of $b_{1}$ and $b_{2}$ will not matter
for the intended asymptotic study in this paper.

This results in the “unrestricted” posterior distribution given by

$\displaystyle\bm{\theta}|(\bm{Y},\sigma)\sim\mathrm{N}_{p}\big{(}\big{(}\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p}\big{)}^{-1}\bm{X}^{\mathrm{T}}\bm{Y},\sigma^{2}\big{(}\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p}\big{)}^{-1}\big{)}.$
(2.5)

Further, the marginal posterior for $\sigma$ is given by

$\displaystyle\sigma^{-2}|\bm{Y}\sim\mathrm{Ga}\big{(}b_{1}+\frac{n}{2},b_{2}+\frac{1}{2}\bm{Y}^{\mathrm{T}}(\bm{I}_{n}-\bm{X}(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}\bm{X}^{\mathrm{T}})\bm{Y}\big{)}.$
(2.6)

We denote the posterior mean
$\big{(}\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p}\big{)}^{-1}\bm{X}^{\mathrm{T}}\bm{Y}$
by $\hat{\bm{\theta}}^{\mathrm{R}}$ owing to its interpretation as the ridge-
regression estimator. The posterior distribution will be denoted by
$\Pi(\cdot|\bm{Y})$.

Note that the posterior distribution disregards variable selection because the
prior does not introduce such a mechanism. The issue is then addressed by
“correcting” the unrestricted posterior through an immersion map that
transforms dense vectors into sparse vectors, in that the distribution induced
by this map from the unrestricted posterior is used for inference. We choose
the sparsity-inducing immersion map given by

$\displaystyle\iota:\bm{\theta}\mapsto\bm{\theta}^{*}\coloneqq\operatorname*{arg\min}_{\bm{u}}\\{n^{-1}\left\|\bm{X}\bm{\theta}-\bm{X}\bm{u}\right\|^{2}+\lambda_{n}\left\|\bm{u}\right\|_{1}\\}.$
(2.7)

We shall call the map $\iota$ the sparse-projection operator. It depends on
the choice of the tuning parameter $\lambda_{n}$. The map is motivated by the
LASSO in that if $\bm{\theta}=\hat{\bm{\theta}}$, the least square estimator,
then $\bm{\theta}^{*}$ is the LASSO estimator $\hat{\bm{\theta}}^{\mathrm{L}}$
defined by

$\displaystyle\hat{\bm{\theta}}^{\mathrm{L}}$
$\displaystyle=\operatorname*{arg\min}_{\bm{u}}\\{n^{-1}\left\|\bm{Y}-\bm{X}\bm{u}\right\|^{2}+\lambda_{n}\left\|\bm{u}\right\|_{1}\\}$
$\displaystyle=\operatorname*{arg\min}_{\bm{u}}\\{n^{-1}\|\bm{X}\hat{\bm{\theta}}-\bm{X}\bm{u}\|^{2}+\lambda_{n}\left\|\bm{u}\right\|_{1}\\}.$
(2.8)

The distribution induced from the unrestricted posterior by the map $\iota$ is
an immersion posterior distribution in the terminology of Wang and Ghosal
[2023a]. In our context, it will be called the sparse-projection posterior
distribution and will be denoted by $\Pi^{*}(\cdot|\bm{Y})$, that is,
$\Pi^{*}(B|\bm{Y})=\Pi(\bm{\theta}^{*}\in B|\bm{Y})$ for
$B\subset\mathbb{R}^{p}$.

The immersion posterior is easy to compute by sampling. Having obtained a
dense posterior sample $\bm{\theta}$ from (2.5), we compute
$\bm{\theta}^{*}=\iota(\bm{\theta})$ and record that as a sample from the
immersion posterior distribution. We repeat the operation independently
sufficiently many times so that posterior probabilities and expectations can
be reliably calculated from sampling. In particular, the method also computes
model posterior probabilities. As the method does not use any Markov chain
Monte Carlo (MCMC) method, it promises faster computation. The method is also
amenable to parallelization since both $\bm{X}^{\mathrm{T}}\bm{X}$ and
$\bm{X}^{\mathrm{T}}\bm{Y}$ can be computed by dividing the dataset into
several parts and computing the sum of products at different machines,
followed by an aggregation step at the central server. This feature is
especially useful if the sample size is huge.

It may be noted that the immersion posterior method may also be used in
conjugation with other sparsity-inducing operators instead of (2.7). Other
penalties that may be used to define an immersion posterior include those
appearing in defining estimators alternative to the LASSO such as the Minimax
Concave Penalty (MCP) [Zhang, 2010], the Smoothly Clipped Absolute Deviation
(SCAD) [Fan and Li, 2001], the Dantzig selector [Candes and Tao, 2007], the
adaptive LASSO [Zou, 2006], the non-negative garrotte estimator [Breiman,
1995, Yuan and Lin, 2007] among others. In the present paper, we forgo other
possible penalty functions and study the asymptotic coverage posterior
credible regions corresponding only to the immersion map (2.7).

##  3 Main Results

Let the predictor dimension $p$ be fixed. We make the following assumptions:

######  Assumption 1 (Predictor).

The matrix $\bm{C}_{n}\to\bm{C}$ for some positive definite matrix $\bm{C}$
and $n^{-1}\max\\{\|\bm{X}_{i}\|^{2}:i=1,\ldots,n\\}\to 0$.

######  Assumption 2 (Tuning).

The tuning parameter $\lambda_{n}$ satisfies
$\lambda_{n}\sqrt{n}\to\lambda_{0}$ for some $\lambda_{0}\geq 0$.

Under the above conditions, Fu and Knight [2000] showed that the LASSO
estimator $\hat{\bm{\theta}}^{\mathrm{L}}$ is weakly consistent and the
limiting distribution of
$\sqrt{n}(\hat{\bm{\theta}}^{\mathrm{L}}-\bm{\theta}^{0})$ puts a positive
probability at zero for the components corresponding to the irrelevant
variables.

######  Theorem 1 (Fu and Knight [2000]).

Let $p$ be fixed, and Assumptions 1 and 2 hold. Then
$\hat{\bm{\xi}}_{n}:=\sqrt{n}(\hat{\bm{\theta}}^{\mathrm{L}}-\bm{\theta}^{0})\rightsquigarrow\bm{\xi}$,
where

$\displaystyle\bm{\xi}=\operatorname*{arg\min}_{\bm{v}\in\mathbb{R}^{p}}\\{\bm{v}^{\mathrm{T}}\bm{C}\bm{v}-2\sigma_{0}\bm{v}^{\mathrm{T}}\bm{C}^{1/2}\bm{\Delta}+\lambda_{0}\big{[}\sum_{j=1}^{s_{0}}v_{j}\mathrm{sign}(\theta^{0}_{j})+\sum_{j=s_{0}+1}^{p}|v_{j}|\big{]}\\},$
(3.1)

with $\bm{\Delta}\sim\mathrm{N}_{p}(\bm{0},\bm{I}_{p})$.

It also follows that under the stated conditions, the variance estimator
$\hat{\sigma}^{2}=n^{-1}\|\bm{Y}-\bm{X}\hat{\bm{\theta}}\|^{2}$ is
$\sqrt{n}$-consistent for $\sigma^{2}$, where $\hat{\bm{\theta}}$ is the least
square estimator or a ridge-regression estimator
$\hat{\bm{\theta}}^{\mathrm{R}}$; see Remark 3. Assumption 1 holds with
probability tending to one for random predictors if the values are sampled
independently from a fixed nonsingular distribution. The second part of the
assumption clearly holds if the predictors are uniformly bounded, a
restriction often imposed by scaling.

Under the same condition on the tuning parameter
${\lambda_{n}}{\sqrt{n}}\to\lambda_{0}$ for some constant $\lambda_{0}\geq 0$,
and assuming that $a_{n}/\sqrt{n}\to 0$, we obtain below a joint weak limit of
the posterior distribution of $\sqrt{n}(\bm{\theta}^{*}-\bm{\theta}^{0})$
given the data in the space $\mathfrak{M}(\mathbb{R}^{p})$, equipped with the
topology of weak convergence, and the distribution of the normalized LASSO
estimator in $\mathbb{R}^{p}$.

######  Theorem 2.

Let Conditions 1 and 2 hold, and $a_{n}/\sqrt{n}\to 0$. Then

$\displaystyle\big{(}\Pi(\sqrt{n}(\bm{\theta}^{*}-\bm{\theta}^{0})\in\cdot|\bm{Y}),\sqrt{n}(\hat{\bm{\theta}}^{\mathrm{L}}-\bm{\theta}^{0})\big{)}\rightsquigarrow\big{(}\mathcal{L}(\bm{T}^{*}\in\cdot|\bm{\Delta}),\bm{\xi}\big{)}$
(3.2)

on the product space $\mathfrak{M}(\mathbb{R}^{p})\times\mathbb{R}^{p}$, where

$\displaystyle\bm{T}^{*}=\operatorname*{arg\min}_{\bm{t}\in\mathbb{R}^{p}}\\{\bm{t}^{\mathrm{T}}\bm{C}\bm{t}-2\bm{t}^{\mathrm{T}}\bm{C}\bm{W}^{*}+\lambda_{0}\big{[}\sum_{j=1}^{s_{0}}t_{j}\textnormal{sign}(\theta^{0}_{j})+\sum_{j=s_{0}+1}^{p}|t_{j}|\big{]}\\},$
(3.3)

for
$\bm{W}^{*}|\bm{\Delta}\sim\mathrm{N}(\sigma_{0}\bm{C}^{-1/2}\bm{\Delta},\sigma_{0}^{2}\bm{C}^{-1})$
and $\bm{\xi}$ is defined in (3.1).

As shifting the random distribution by
$\hat{\bm{\xi}}_{n}=\sqrt{n}(\hat{\bm{\theta}}_{n}-\bm{\theta}^{0})$ is a
continuous operation in $\mathfrak{M}(\mathbb{R}^{p})$, we conclude the
following.

######  Corollary 1.

Under the conditions of Theorem 2,

$\displaystyle\big{(}\Pi(\sqrt{n}(\bm{\theta}^{*}-\hat{\bm{\theta}}^{\mathrm{L}})\in\cdot|\bm{Y}),\sqrt{n}(\hat{\bm{\theta}}^{\mathrm{L}}-\bm{\theta}^{0})\big{)}\rightsquigarrow\big{(}\mathcal{L}(\bm{T}^{*}-\bm{\xi}\in\cdot|\bm{\Delta}),\bm{\xi}\big{)}.$
(3.4)

We show in Lemma 3 that when $\theta^{0}_{j}=0$ for all $j=s_{0}+1,\dots,p$,
the minimizer $\bm{T}^{*}$ assumes the value $0$ for its last $p-s_{0}$
components with positive probability. The distribution has no closed-form
expression except for some special cases, and the statement is only marginally
useful. However, we can compute the limiting frequentist coverage of a
Bayesian credible region defined by the quantiles of the posterior
distribution of $\sqrt{n}(\bm{\theta}^{*}-\bm{\theta}^{0})$, which can be
obtained from Theorem 2.

For $\bar{\bm{\theta}}\in\mathbb{R}^{p}$ and $r>0$, let

$\displaystyle
B_{\mathcal{K}}(\bar{\bm{\theta}},r)=\\{\bm{\theta}:\|\sqrt{n}(\bm{\theta}-\bar{\bm{\theta}})\|_{\mathcal{K}}\leq
r\\},$ (3.5)

where the seminorm $\|\cdot\|_{\mathcal{K}}$ is defined by

$\displaystyle\|\bm{x}\|_{\mathcal{K}}=\inf\\{\lambda>0:\bm{x}/{\lambda}\in\mathcal{K}\\},$
(3.6)

the Minkowski functional of a convex set $\mathcal{K}$ containing $0$ in its
interior. Define a $(1-\alpha)$-credible set to be
$B_{\mathcal{K}}(\hat{\bm{\theta}}^{\mathrm{L}},r_{1-\alpha})$, where
$r_{1-\alpha}$ is the $(1-\alpha)$-quantile of the posterior distribution of
$\|\sqrt{n}(\bm{\theta}^{*}-\hat{\bm{\theta}}^{\mathrm{L}})\|_{\mathcal{K}}$,
that is

$\displaystyle\Pi^{*}(B_{\mathcal{K}}(\hat{\bm{\theta}}^{\mathrm{L}},r_{1-\alpha})|\bm{Y})=\Pi(\bm{\theta}^{*}:\|\sqrt{n}(\bm{\theta}^{*}-\hat{\bm{\theta}}^{\mathrm{L}})\|_{\mathcal{K}}\leq
r_{1-\alpha}|\bm{Y})=1-\alpha.$ (3.7)

If the equality cannot be obtained, we choose $r_{1-\alpha}$ the infimum of
all $r>0$ such that
$\Pi^{*}(B_{\mathcal{K}}(\hat{\bm{\theta}}^{\mathrm{L}},r)|\bm{Y})>1-\alpha$.
For different choices of $\mathcal{K}$, we can have coverage results for
different norms or pseudo-norms. For example, choosing
$\mathcal{K}=[-1,1]^{p}$, we get the max norm $\|\cdot\|_{\infty}$, choosing
$\mathcal{K}$ as a unit sphere, we get the unit norm $\|\cdot\|_{2}$, choosing
$\mathcal{K}$ as a unit diamond, we get the $\ell_{1}$-norm $\|\cdot\|_{1}$,
and choosing
$\mathcal{K}=\mathbb{R}\times\cdots\times\mathbb{R}\times[-1,1]\times\mathbb{R}\times\cdots\times\mathbb{R}$,
we get credible intervals for individual components.

######  Theorem 3.

For a given level of confidence $(1-\alpha)$ and
$\bm{\Delta}\sim\mathrm{N}_{p}(\bm{0},\bm{I}_{p})$, we have under Assumptions
1 and 2 and $a_{n}/\sqrt{n}\to 0$,

$\displaystyle\mathrm{P}\left(\mathrm{P}(\|\bm{T}^{*}-\bm{\xi}\|_{\mathcal{K}}\leq\|\bm{\xi}\|_{\mathcal{K}}\big{|}\bm{\Delta})<1-\alpha\right)$
(3.8)
$\displaystyle\leq\liminf_{n\to\infty}\mathbb{P}_{0}\big{(}\bm{\theta}^{0}\in
B_{\mathcal{K}}(\hat{\bm{\theta}}^{\mathrm{L}},r_{1-\alpha})\big{)}$
$\displaystyle\leq\limsup_{n\to\infty}\mathbb{P}_{0}\big{(}\bm{\theta}^{0}\in
B_{\mathcal{K}}(\hat{\bm{\theta}}^{\mathrm{L}},r_{1-\alpha})\big{)}$
$\displaystyle\leq\mathrm{P}\left(\mathrm{P}(\|\bm{T}^{*}-\bm{\xi}\|_{\mathcal{K}}\leq\|\bm{\xi}\|_{\mathcal{K}}\big{|}\bm{\Delta})\leq
1-\alpha\right),$

where $\bm{T}^{*}$ and $\bm{\xi}$ are as in Theorem 2.

The limiting coverage of a sparse-projection posterior credible region is thus
assured to be at least the probability on the right side, but it may not equal
the credibility level $1-\alpha$. The limiting coverage depends on $\bm{C}$,
$\sigma_{0}$, $\lambda_{0}$ and $\alpha$. In the following subsection, we
evaluate the limiting coverage explicitly for a credible interval for a
component under an asymptotic orthogonality setting.

###  3.1 Special Case: Asymptotically Uncorrelated Predictor

We now explicitly evaluate the limiting coverage of a credible interval of a
single regression coefficient $\theta_{j}$, assuming that the $j$th predictor
is asymptotically uncorrelated with the rest, for some $j=1,\ldots,p$.

######  Assumption 3 (Asymptotically uncorrelated predictor).

The $j$th predictor variable $X_{j}$ is asymptotically uncorrelated with the
remaining predictors $\\{X_{1},\ldots,X_{j-1},X_{j+1},\ldots,X_{p}\\}$, that
is, $\bm{C}=\lim_{n\to\infty}n^{-1}\bm{X}^{\mathrm{T}}\bm{X}$ is block-
diagonal with blocks $\\{j\\}$ and $[-j]:=\\{1,\ldots,j-1,j+1,\ldots,p\\}$.

Let $c_{j}$ stand for the $j$th diagonal entry of $\bm{C}$. Note that, under
the above assumption, all other entries of the $j$th row and $j$th column of
$\bm{C}$ are 0.

Although the assumption of asymptotic uncorrelatedness is very special, it may
naturally hold for predictors not related to other predictors. Moreover, the
condition can be enforced by an orthogonalization technique, but it then
changes the meaning of the coefficient.

Consider the credible interval
$[\hat{\theta}^{\mathrm{L}}_{j}-r_{1-\alpha,j},\hat{\theta}^{\mathrm{L}}_{j}+r_{1-\alpha,j}]$,
where $\hat{\theta}^{\mathrm{L}}_{j}$ is the $j$th coordinate of the LASSO
estimator $\hat{\bm{\theta}}^{\mathrm{L}}$ and $r_{1-\alpha,j}$ is the
$(1-\alpha)$-quantile of the induced posterior distribution of
$|\theta^{*}_{j}-\hat{\theta}_{j}^{\mathrm{L}}|$. This credible interval
corresponds to the credible region in (3.5) with
$\mathcal{K}=\mathbb{R}\times\cdots\times\mathbb{R}\times[-1,1]\times\mathbb{R}\times\cdots\times\mathbb{R}$,
where $[-1,1]$ appears at the $j$th coordinate. Under Assumption 3, it follows
from (3.1) that the $j$th component $\xi_{j}$ of $\bm{\xi}$ is given by

$\displaystyle\xi_{j}$
$\displaystyle=\mathrm{argmin}\\{c_{j}v_{j}^{2}-2\sigma_{0}\sqrt{c_{j}}\Delta_{j}v_{j}+\lambda_{0}[v_{j}\mathrm{sign}(\theta_{j}^{0})+|v_{j}|\mathbbm{1}(\theta_{j}=0)]:v_{j}\in\mathbb{R}\\}$
$\displaystyle=\begin{cases}\frac{\sigma_{0}}{\sqrt{c_{j}}}\Delta_{j}-\frac{\lambda_{0}}{2c_{j}},&\mbox{
if }\theta_{j}^{0}>0,\\\
\frac{\sigma_{0}}{\sqrt{c_{j}}}\Delta_{j}+\frac{\lambda_{0}}{2c_{j}},&\mbox{
if }\theta_{j}^{0}<0,\\\
(\frac{\sigma_{0}}{\sqrt{c_{j}}}\Delta_{j}-\frac{\lambda_{0}}{2c_{j}})\mathbbm{1}\\{\Delta_{j}>\frac{\lambda_{0}}{2\sigma_{0}\sqrt{c_{j}}}\\}\\\
+(\frac{\sigma_{0}}{\sqrt{c_{j}}}\Delta_{j}+\frac{\lambda_{0}}{2c_{j}})\mathbbm{1}\\{\Delta_{j}<-\frac{\lambda_{0}}{2\sigma_{0}\sqrt{c_{j}}}\\},&\mbox{
if }\theta_{j}^{0}=0.\end{cases}$ (3.9)

Thus, $\xi_{j}$ depends on $\bm{\Delta}$ only through its $j$ component
$\Delta_{j}$ and is independent of the remaining components.

Similarly, from (3.3), it follows that the $j$th coordinate $T_{j}^{*}$ of
$\bm{T}^{*}$ is given by

$\displaystyle T_{j}^{*}$
$\displaystyle=\mathrm{argmin}\\{c_{j}t_{j}^{2}-2c_{j}W_{j}^{*}t_{j}+\lambda_{0}[t_{j}\mathrm{sign}(\theta_{j}^{0})+|t_{j}|\mathbbm{1}(\theta_{j}=0)]:t_{j}\in\mathbb{R}\\}$
$\displaystyle=\begin{cases}W_{j}^{*}-\frac{\lambda_{0}}{2c_{j}},&\mbox{ if
}\theta_{j}^{0}>0,\\\ W_{j}^{*}+\frac{\lambda_{0}}{2c_{j}},&\mbox{ if
}\theta_{j}^{0}<0,\\\
(W_{j}^{*}-\frac{\lambda_{0}}{2c_{j}})\mathbbm{1}\\{W_{j}^{*}>\frac{\lambda_{0}}{2c_{j}}\\}\\\
+(W_{j}^{*}+\frac{\lambda_{0}}{2c_{j}})\mathbbm{1}\\{W_{j}^{*}<-\frac{\lambda_{0}}{2c_{j}}\\},&\mbox{
if }\theta_{j}^{0}=0,\end{cases}$ (3.10)

where $W_{j}^{*}$ is the $j$th component of $\bm{W}^{*}$, and that
$W_{j}^{*}|\Delta_{j}\sim\mathrm{N}(\Delta_{j},1)$.

We then introduce the notations

$\displaystyle h_{+}(\lambda_{0},\zeta)$
$\displaystyle=2\Phi(|\zeta-\lambda_{0}/2|)-1,$ $\displaystyle
h_{-}(\lambda_{0},\zeta)$
$\displaystyle=2\Phi(|\zeta+\lambda_{0}/2|)-1=h_{+}(\lambda_{0},-\zeta),$
$\displaystyle h_{0}(\lambda_{0},\zeta)$
$\displaystyle=\begin{cases}\Phi(\zeta-\lambda_{0}/2)-\Phi(-\zeta-\lambda_{0}/2),&\mbox{
if }\zeta>\lambda_{0}/2,\\\
\Phi(-\zeta+\lambda_{0}/2)-\Phi(\zeta+\lambda_{0}/2),&\mbox{ if
}\zeta<\lambda_{0}/2,\\\
\Phi(\zeta+\lambda_{0}/2)-\Phi(\zeta-\lambda_{0}/2),&\mbox{ if
}|\zeta|\leq\lambda_{0}/2.\\\ \end{cases}$ (3.11)

and

$\displaystyle\psi(\alpha,\lambda_{0})$
$\displaystyle=\Phi(\lambda_{0}/2+z_{\alpha/2})-\Phi(\lambda_{0}/2-z_{\alpha/2}),$
(3.12) $\displaystyle\psi_{0}(\alpha,\lambda_{0})$
$\displaystyle=\int_{-\infty}^{\infty}\mathbbm{1}\big{(}h_{0}(\lambda_{0},\zeta)\leq
1-\alpha\big{)}\phi(\zeta)d\zeta.$ (3.13)

The asymptotic coverage of the credible interval
$[\hat{\theta}^{\mathrm{L}}_{j}-r_{1-\alpha,j},\hat{\theta}^{\mathrm{L}}_{j}+r_{1-\alpha,j}]$
for $\theta_{j}$ is characterized in the following theorem.

######  Theorem 4.

If Assumptions 1, 2 and 3 hold and $a_{n}/\sqrt{n}\to 0$, then for a given
level of confidence $0<1-\alpha<1$,

$\displaystyle\mathbb{P}_{0}\big{(}\theta_{0j}\in[\hat{\theta}_{j}^{\mathrm{L}}-r_{1-\alpha,j},\hat{\theta}_{j}^{\mathrm{L}}+r_{1-\alpha,j}]\big{)}\to\begin{cases}\psi(\alpha,\lambda_{0}\sqrt{c_{j}}/\sigma_{0}),&\mbox{
if }\theta_{j}\neq 0,\\\
\psi_{0}(\alpha,\lambda_{0}\sqrt{c_{j}}/\sigma_{0}),&\mbox{ if
}\theta_{j}=0.\end{cases}$ (3.14)

######  Remark 1.

The original data generating process

$Y_{i}=\sum_{k=1}^{p}X_{ik}\theta_{k}+\varepsilon_{i}=X_{ij}\theta_{j}+\sum_{k=1:k\neq
j}^{p}X_{ik}\theta_{k}+\varepsilon_{i}$

can be rewritten as

$\displaystyle\tilde{Y}_{i}=\tilde{X}_{ij}\tilde{\theta}_{0j}+\sum_{k=1:k\neq
j}^{p}X_{ik}\frac{\theta_{0k}}{\sigma_{0}}+\tilde{\varepsilon}_{i},$ (3.15)

where $\tilde{Y}_{i}=Y_{i}/\sigma_{0}$, $\tilde{X}_{ij}=x_{ij}/\sqrt{c_{j}}$,
$\tilde{\varepsilon}_{i}=\varepsilon_{i}/\sigma_{0}$ and
$\tilde{\theta}_{j}^{0}=\theta_{j}^{0}\sqrt{c_{j}}/\sigma_{0}$. For the
transformed process, the $j$th predictor $\tilde{X}_{j}$ still remains
asymptotically uncorrelated with the remaining predictors, and the
corresponding $\bm{C}$-matrix has the $j$th diagonal entry 1. The true error
variance also reduces to 1. Thus, the asymptotic coverage for the credible
interval for $\theta_{j}$ using the penalty parameter $\lambda_{0}$ coincides
with that for $\tilde{\theta}_{j}=\theta_{j}\sigma_{0}/\sqrt{c_{j}}$ using the
penalty parameter $\lambda_{0}\sqrt{c_{j}}/\sigma_{0}$. Hence, to prove
Theorem 4, we can assume, without loss of generality, that $c_{j}=1$ and
$\sigma_{0}=1$.

Below, in Figure 1, we provide the plots of $h_{+}(\lambda_{0},\zeta)$ and
$h_{0}(\lambda_{0},\zeta)$ against $\zeta$ for different choices of the
limiting penalty $\lambda_{0}$ as well as the limiting coverage probabilities
against the confidence level for the signal and noise components. Since
$\psi(\alpha,\lambda_{0})\leq 1-\alpha$ with equality holding only for
$\lambda_{0}=0$, for the signal components, no non-zero value of the penalty
parameter $\lambda_{0}$ can lead to the exact $(1-\alpha)$ asymptotic
coverage, or more, of a $(1-\alpha)$ sparse-projection posterior credible
interval. However, by fixing a penalty parameter $\lambda_{0}$, we can ensure
that for some $\gamma$ depending on $\lambda_{0}$ and $\alpha$, a $(1-\gamma)$
credible interval will have exactly $(1-\alpha)$ asymptotic coverage if
$\theta^{0}_{j}\neq 0$. Moreover, in Figure 2, the function
$\psi(\alpha,\lambda_{0})$ is shown to be uniformly bounded above by
$\psi_{0}(\alpha,\lambda_{0})$ for a given $\alpha$, meaning that the limiting
coverage of a noise component will always be higher than that of a signal
component. Thus, owing to the fact that
$\psi_{0}(\alpha,\lambda_{0})\geq\psi(\alpha,\lambda_{0})$, as shown by
numerical evaluation in Figure 2, it follows that the asymptotic coverage of
the corresponding interval is at least $(1-\alpha)$.

Figure 1: The upper panel shows the graph of $h_{+}(\lambda_{0},\zeta)$ for
the positive coefficients (left), and that of $h(\lambda_{0},\zeta)$ for the
zero-valued components (right) as a function of $\zeta$ for various values of
$\lambda_{0}$. The lower panel shows coverage as a function of the credibility
level for non-zero (left) and zero (right) components for various values
$\lambda_{0}$. The black diagonal line corresponds to the exact coverage.
Figure 2: Plot of coverage probabilities against $\lambda$ for positive
components (dashed lines) and zero components (solid lines) for a few
interesting values of $1-\alpha$.

Based on the above discussion and the figures displayed above, we provide a
guideline for choosing the confidence level for a wide range of coverage
requirements. In Table 1, we give the value of $\gamma$ that attains the
limiting coverage $1-\alpha$ for any penalty parameter between 0.05 to 4. For
any value of the penalty parameter not listed in Table 1, the solve_gamma
function in the credInt R package [Pal, 2024] may be used. Moreover, from the
plot of $\psi_{0}(\alpha,\lambda)$ in Figure 1, we can ensure higher than
$(1-\alpha)$-coverage for the zero components for the chosen $(1-\gamma)$
level credible intervals. Thus, choosing $\gamma$ according to the table
attains the limiting coverage exactly $1-\alpha$ for a non-zero parameter
value and at least $1-\alpha$ for a zero value of the parameter.

$\lambda$ $1-\alpha$ | 0.9 | 0.925 | 0.95 | 0.975 | 0.99
---|---|---|---|---|---
0.05 | 0.9001 | 0.9251 | 0.9501 | 0.9751 | 0.9900
0.1 | 0.9004 | 0.9254 | 0.9503 | 0.9752 | 0.9901
0.15 | 0.9009 | 0.9258 | 0.9506 | 0.9754 | 0.9902
0.2 | 0.9017 | 0.9264 | 0.9511 | 0.9757 | 0.9904
0.25 | 0.9026 | 0.9272 | 0.9517 | 0.9761 | 0.9906
0.3 | 0.9037 | 0.9282 | 0.9525 | 0.9765 | 0.9909
0.35 | 0.9050 | 0.9293 | 0.9534 | 0.9771 | 0.9910
0.4 | 0.9066 | 0.9306 | 0.9543 | 0.9777 | 0.9914
0.45 | 0.9082 | 0.9320 | 0.9554 | 0.9784 | 0.9917
0.5 | 0.9100 | 0.9335 | 0.9566 | 0.9790 | 0.9920
0.55 | 0.9120 | 0.9352 | 0.9578 | 0.9798 | 0.9924
0.6 | 0.9141 | 0.9369 | 0.9591 | 0.9806 | 0.9927
0.65 | 0.9163 | 0.9387 | 0.9605 | 0.9814 | 0.9931
0.7 | 0.9186 | 0.9406 | 0.9619 | 0.9822 | 0.9934
0.75 | 0.9210 | 0.9426 | 0.9634 | 0.9830 | 0.9938
0.8 | 0.9235 | 0.9446 | 0.9649 | 0.9838 | 0.9942
0.85 | 0.9262 | 0.9468 | 0.9664 | 0.9846 | 0.9945
0.9 | 0.9288 | 0.9488 | 0.9679 | 0.9854 | 0.9948
0.95 | 0.9314 | 0.9510 | 0.9694 | 0.9862 | 0.9952
1 | 0.9340 | 0.9530 | 0.9708 | 0.9871 | 0.9955
1.1 | 0.9394 | 0.9572 | 0.9737 | 0.9885 | 0.9961
1.2 | 0.9446 | 0.9613 | 0.9765 | 0.9899 | 0.9966
1.3 | 0.9498 | 0.9652 | 0.9791 | 0.9912 | 0.9971
1.4 | 0.9546 | 0.9688 | 0.9815 | 0.9923 | 0.9976
1.5 | 0.9593 | 0.9722 | 0.9837 | 0.9934 | 0.9979
1.6 | 0.9636 | 0.9754 | 0.9857 | 0.9943 | 0.9982
1.7 | 0.9676 | 0.9783 | 0.9875 | 0.9951 | 0.9985
1.8 | 0.9713 | 0.9810 | 0.9891 | 0.9958 | 0.9987
1.9 | 0.9746 | 0.9833 | 0.9906 | 0.9964 | 0.9989
2 | 0.9776 | 0.9854 | 0.9918 | 0.9959 | 0.9991
2.2 | 0.9828 | 0.9889 | 0.9939 | 0.9978 | 0.9994
2.4 | 0.9869 | 0.9917 | 0.9956 | 0.9984 | 0.9996
2.6 | 0.9901 | 0.9939 | 0.9968 | 0.9988 | 0.9997
2.8 | 0.9927 | 0.9955 | 0.9977 | 0.9992 | 0.9998
3 | 0.9946 | 0.9967 | 0.9983 | 0.9994 | 0.9998
3.5 | 0.9976 | 0.9986 | 0.9993 | 0.9998 | 0.9999
4 | 0.9990 | 0.9994 | 0.9997 | 0.9999 | 1.000
Table 1: Calibration table for credibility level for an intended asymptotic
coverage.

######  Remark 2.

The component-wise credible intervals can be further used to get joint
credible hyper-rectangle. Fix $k\leq p$ components. Suppose we want a
$1-\beta$ credible region for those $k$ components. Then, we first build the
individual $1-\alpha$ credible intervals for each of those $k$ chosen
components using the sparse-projection posterior method with the corresponding
credibility levels calculated from Table 1. Using the asymptotic independence
resulting from the diagonal nature of the limiting covariance matrix, the
asymptotic coverage of the joint $k-$dimensional hyper-rectangle thus formed
is $(1-\alpha)^{k}$. Consequently, choosing $1-\alpha=(1-\beta)^{1/k}$, the
targeted asymptotic coverage $1-\beta$ for the joint credible region can be
achieved.

##  4 Numerical Results

We conduct a small-scale simulation to capture the performance of the credible
sets obtained by projecting the conjugate posterior samples to the sparse
subspace. We study several scenarios and report a complete numerical
performance encompassing the dimension, the sparsity, and the signal strength
of the underlying model. First, fixing the true sparsity at $s=5$ and
dimension at $p=20$, we first generate data using an independent design matrix
and signal strength $(-2,-1.5,0.5,1,2)$. We repeat this study with a
correlated design matrix where an autoregressive process with a correlation of
$0.7$ has been used. We study the coverage and length of the $(1-\gamma)$
credible intervals of the 5 signals and 5 randomly chosen noise variables.
Here, $\gamma$ is chosen in accordance with Table 1 to achieve our target
coverage of 0.95. We will need the value of the penalty parameter used in the
projection map to use the calibration information in Table 1. To specifically
know what $\gamma$ value returns the credible intervals with exact
$(1-\alpha)$ coverage, the solve_gamma function within the credInt R package
[Pal, 2024] can be used. For all the simulations, we implement the LASSO
regression of $\bm{Y}$ on $\bm{X}$ and use the cross-validated penalty
parameter in all the $R=5000$ posterior maps. These cases are observed under
three sample sizes, $n=500,\,1000$, and $5000$. The credInt function returns
the $(1-\gamma)$ credible intervals for all the regression coefficients. We
repeat the simulations $200$ times and report the average coverage values. The
goal is to see if the coverage of the component-wise credible intervals
increases and attains the desired level with the increase in sample size.

Besides this, we study the impact of sparsity on coverage for different
significance levels. We vary $s_{0}$ from 5 to 30 and consider three levels,
90%, 95%, and 99%, to see how coverage changes as sparsity increases and
whether the results improve when the sample size increases from 500 to 1000.
Note that the design matrix in this case is independent, and the error is
normally distributed with error variance 1.

Design | Model |  | $S_{1}$ | $S_{2}$ | $S_{3}$ | $S_{4}$ | $S_{5}$ | $N_{1}$ | $N_{2}$ | $N_{3}$ | $N_{4}$ | $N_{5}$
---|---|---|---|---|---|---|---|---|---|---|---|---
Independent | n = 500 | Coverage | 0.943 | 0.938 | 0.942 | 0.935 | 0.937 | 0.948 | 0.944 | 0.952 | 0.943 | 0.947
p = 20 | Length | 0.112 | 0.108 | 0.119 | 0.110 | 0.111 | 0.113 | 0.109 | 0.112 | 0.117 | 0.116
n = 1000 | Coverage | 0.949 | 0.945 | 0.947 | 0.951 | 0.948 | 0.955 | 0.953 | 0.954 | 0.950 | 0.957
p = 20 | Length | 0.117 | 0.114 | 0.114 | 0.113 | 0.115 | 0.113 | 0.114 | 0.116 | 0.113 | 0.114
n = 5000 | Coverage | 0.955 | 0.951 | 0.949 | 0.951 | 0.952 | 0.958 | 0.957 | 0.962 | 0.960 | 0.961
p = 20 | Length | 0.114 | 0.119 | 0.110 | 0.116 | 0.115 | 0.115 | 0.116 | 0.118 | 0.119 | 0.116
Correlated | n = 500 | Coverage | 0.933 | 0.937 | 0.931 | 0.933 | 0.932 | 0.945 | 0.944 | 0.944 | 0.941 | 0.942
p = 20 | Length | 0.119 | 0.121 | 0.117 | 0.118 | 0.119 | 0.119 | 0.115 | 0.116 | 0.118 | 0.120
n = 1000 | Coverage | 0.943 | 0.948 | 0.939 | 0.942 | 0.942 | 0.956 | 0.954 | 0.953 | 0.950 | 0.953
p = 20 | Length | 0.121 | 0.122 | 0.116 | 0.117 | 0.119 | 0.117 | 0.112 | 0.116 | 0.122 | 0.121
n = 5000 | Coverage | 0.951 | 0.952 | 0.947 | 0.949 | 0.950 | 0.957 | 0.958 | 0.958 | 0.957 | 0.960
p = 20 | Length | 0.118 | 0.117 | 0.115 | 0.120 | 0.118 | 0.130 | 0.127 | 0.123 | 0.119 | 0.125

Table 2: Coverage and length of $(1-\gamma)$ credible interval of $5$ signals
and randomly chosen $5$ noise variables studied over $3$ sample sizes, under
the independent and correlated design settings with target coverage $0.95$.

Table 2 provides the coverage and lengths of the component-wise intervals for
the cases ($p=20$ and $\bm{\theta}^{0}_{S_{0}}=(-2,-1.5,0.5,1,1.5)$) under the
two setups, namely, independent design and correlated design. We observe that
in both cases, there is a general pattern that coverage increases with an
increase in sample size and almost attains the desired level, which is 95% in
this case. The numerical results thus align with the theoretical findings and
provide enough support that the projection method leads to almost accurate
asymptotic coverage. It is noteworthy that although the theoretical results
for this method have only been proven under a limiting independence condition
on the covariates, the method is capable of performing well even when some
correlation is present between the explanatory variables in the model.

Figure 3, on the other hand, depicts the coverage property as the model
becomes sparser. We see the general tendency is that the coverage slightly
degrades as the number of truly active predictors decreases. However, the
solid line (showing the coverage for $n=1000$) stays uniformly above the
dashed line (showing coverage for $n=500$) for all three significance levels.
This guarantees that the asymptotic full coverage will hold as the sample size
grows. Moreover, the sparse projection-posterior method is much more flexible
in terms of coverage as the choice of the $\lambda_{0}$ values is user-
dependent.

Figure 3: The average signal coverage of $\tau$-credible regions for three
different levels $\tau=0.9,0.95,0.99$ (averaged over the $s$ signals) are
plotted against the degree of sparsity $s=30,25,20,15,10,5$ with two choices
of the sample size $n=500,1000$. The signal strength in all cases was
$\theta^{0}_{j}=1$ for all $j\in S=\\{1,\dots,s\\}$ and the results are based
on $200$ replications. Figure 4: Comparison of Horseshoe, Bayesian LASSO,
Bootstrap LASSO average signal (upper panel) and noise (lower panel) coverage
with dimension $p=50$, sparsity $s=20$ and signal intensity
$\beta_{j}^{0}=1.5$ for all $j=1,2,\dots,s$, corresponding to three sample
sizes, $n=500,1000$ and $2000$. Figure 5: Comparison of Horseshoe, Bayesian
LASSO, Bootstrap LASSO average signal (upper panel) and noise (lower panel)
coverage with dimension $p=100$, sparsity $s=20$ and signal intensity
$\beta_{j}^{0}=1.5$ for all $j=1,2,\dots,s$, corresponding to three sample
sizes, $n=500,1000$ and $2000$.

Next, in Figures 4 and 5, we compare the average signal and noise coverage and
length of credible intervals for the Bayesian methods of Bayesian LASSO,
Horseshoe, and the proposed sparse projection-posterior method, along with the
frequentist method of Bootstrap LASSO. We repeat the study for three levels of
significance $1-\alpha=0.9,0.95,0.99$, and for two dimensions $p=50,100$. We
observe whether there is any change in the coverage over three sample sizes
$n=500,1000,5000$, keeping the sparsity and the signal intensity fixed at
$s=20$ and $\beta_{j}^{0}=1.5,\enskip j=1,2,\dots,s$, respectively. Again, all
the projections computed used the LASSO cross-validated penalty parameter. We
drew $R=5000$ posterior samples. The corresponding $\gamma$ was chosen from
Table 1. Although the Bayesian LASSO does well in covering the true parameter
values for low dimensions, there are two more aspects to consider. Firstly,
the Bayesian LASSO does not have a model selection ability and always returns
the full feature space. This causes a problem in identifying weaker signals.
Secondly, the method is computationally infeasible when the dimension exceeds
100. On the other hand, the horseshoe method provides over-coverage for any
choice of $1-\alpha$ at the high cost of yielding very long credible
intervals. Accounting for these drawbacks, the sparse projection-posterior
method is the only Bayesian method capable of providing exact coverage with
reasonable lengths and comparable to the bootstrap LASSO technique of
Chatterjee and Lahiri [2011]. Its coverage is sometimes even better than the
bootstrap LASSO method, especially for larger sample sizes. Moreover, the
bootstrap LASSO often fails to cover the noise variables, but the sparse-
projection posterior method successfully provides coverage. As predicted by
the theory, the cost is a slightly elongated interval.

##  5 Conclusion

This study shows that the proposed Bayesian technique of projecting a dense
conjugate posterior sample to a sparser domain performs well, offering
adequate coverage by a modified credible under the fixed-dimensional setting.
The method’s ability to achieve reliable and consistent coverage holds
promising implications for applications in various fields, as it successfully
overcomes the drawbacks of existing methodologies to quantify the associated
uncertainties in sparse regression. This technique may be extended to address
the high dimensional setting and correlated predictors, and the method’s
robustness under different data distributions may be explored.

## Appendix: Proofs

######  Lemma 1.

Under Assumptions 1 and 2, and $a_{n}/\sqrt{n}\to 0$, we have that
$\sqrt{n}(\hat{\bm{\theta}}^{\mathrm{R}}-\bm{\theta}^{0})\rightsquigarrow\mathrm{N}_{p}(\bm{0},\sigma_{0}^{2}\bm{C}^{-1})$.

  * Proof.If $\bm{\varepsilon}$ is normally distributed, representing

$\hat{\bm{\theta}}^{\mathrm{R}}=(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}\bm{X}^{\mathrm{T}}\bm{X}\bm{\theta}^{0}+(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}\bm{X}^{\mathrm{T}}\bm{\varepsilon},$

we observe that $\hat{\bm{\theta}}^{\mathrm{R}}$ is also normally distributed
and hence it suffices to show that

$\displaystyle\sqrt{n}(\mathbb{E}_{0}(\hat{\bm{\theta}}^{\mathrm{R}})-\bm{\theta}^{0})$
$\displaystyle=\sqrt{n}((\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}\bm{X}^{\mathrm{T}}\bm{X}\bm{\theta}^{0}-\bm{\theta}^{0})\to\bm{0},$
$\displaystyle
n\mathbb{E}_{0}[(\hat{\bm{\theta}}^{\mathrm{R}}-\bm{\theta}^{0})(\hat{\bm{\theta}}^{\mathrm{R}}-\bm{\theta}^{0})^{\mathrm{T}}]$
$\displaystyle=n\sigma_{0}^{2}(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}\bm{X}^{\mathrm{T}}\bm{X}(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}$
$\displaystyle\to\sigma_{0}^{2}\bm{C}^{-1}.$

Since $n^{-1}\bm{X}^{\mathrm{T}}\bm{X}\to\bm{C}$ and $a_{n}/\sqrt{n}\to 0$, we
have
$(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}\bm{X}^{\mathrm{T}}\bm{X}=\bm{I}_{p}+\bm{o}(n^{-1/2})$;
here and elsewhere $\bm{o}(n^{-1/2})$ stands for a $p\times p$-matrix with all
entries $o(n^{-1/2})$. The first relation now follows. The second relation
holds even under the weaker condition $a_{n}/n\to 0$.

In general, when $\bm{\varepsilon}$ may not be normal, we verify, for any
fixed non-zero linear combination $\bm{b}^{\mathrm{T}}\hat{\bm{\theta}}$,
Lindeberg’s condition for the central limit theorem. The condition reduces to
$(\max_{1\leq i\leq
n}|\bm{b}^{\mathrm{T}}\bm{X}_{i}|^{2})/\sum_{i=1}^{n}|\bm{b}^{\mathrm{T}}\bm{X}_{i}|^{2}\to
0$, owing to the identical distribution of the error variables, invoke the
Cramér-Wold device. Since

$\max_{1\leq i\leq
n}n^{-1}|\bm{b}^{\mathrm{T}}\bm{X}_{i}|^{2}\leq\|\bm{b}\|^{2}\max_{1\leq i\leq
n}n^{-1}|\bm{X}_{i}\|^{2}\to
0,\;n^{-1}\sum_{i=1}^{n}|\bm{b}^{\mathrm{T}}\bm{X}_{i}|^{2}\to\bm{b}^{\mathrm{T}}\bm{C}\bm{b}>0$

by Assumption 1, Lindeberg’s condition is verified. As

$\mathrm{Cov}(n^{-1/2}\sum_{i=1}^{n}x_{ij}\varepsilon_{i},n^{-1/2}\sum_{i=1}^{n}x_{ik}\varepsilon_{i})=n^{-1}\sigma_{0}^{2}\sum_{i=1}^{n}x_{ij}x_{ik}\to\sigma_{0}^{2}c_{jk},$

the multivariate central limit theorem applies to
$(n^{-1/2}\sum_{i=1}^{n}x_{ij}\varepsilon_{i}:1\leq j\leq p)$ with mean the
zero vector and the dispersion matrix $\sigma_{0}^{2}\bm{C}$. Hence it follows
that
$(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}\bm{X}^{\mathrm{T}}\bm{\varepsilon}\rightsquigarrow\mathrm{N}_{p}(\bm{0},\sigma_{0}^{2}\bm{C}^{-1})$.
∎

######  Lemma 2.

Under Assumptions 1 and 2, and $a_{n}/\sqrt{n}\to 0$, the posterior
distribution of $\sigma$ contracts at $\sigma_{0}$ at the rate $n^{-1/2}$,
that is, for any $M_{n}\to\infty$

$\displaystyle\mathbb{E}_{0}\Pi(|\sigma-\sigma_{0}|>M_{n}n^{-1/2}|\bm{Y})\to
0.$ (.1)

  * Proof.It suffices to prove that the posterior for $\tau=\sigma^{-2}$ concentrates in $n^{-1/2}$-sized neighborhoods of $\tau_{0}=\sigma_{0}^{-2}$, which we verify from (2.6) using Chebyshev’s inequality. We show that $\mathbb{E}(\tau|\bm{Y})=\tau_{0}+O_{p}(n^{-1/2})$ and $\mathrm{Var}(\tau|\bm{Y})=O_{p}(n^{-1})$. Let $\bm{\delta}_{n}=n^{-1}\bm{X}^{\mathrm{T}}\bm{\varepsilon}$. Then, we can write $(b_{1}+n/2)[\mathbb{E}(\tau|\bm{Y})]^{-1}/(n/2)-2b_{2}/n$ as

$\displaystyle(\bm{\theta}^{0})^{\mathrm{T}}\bm{C}_{n}\bm{\theta}^{0}+2\bm{\delta}_{n}^{\mathrm{T}}\bm{\theta}^{0}+\|\bm{\varepsilon}\|^{2}/n-(\bm{C}_{n}\bm{\theta}^{0}+\bm{\delta}_{n})^{\mathrm{T}}(\bm{C}_{n}+n^{-1}a_{n}\bm{I}_{p})^{-1}(\bm{C}_{n}\bm{\theta}^{0}+\bm{\delta}_{n}).$

This reduces to

$\displaystyle\sigma_{0}^{2}$
$\displaystyle+O_{p}(n^{-1/2})-(\bm{\theta}^{0})^{\mathrm{T}}(\bm{C}_{n}+o(n^{-1/2}))\bm{\theta}^{0}$
$\displaystyle\quad-2\bm{\delta}_{n}^{\mathrm{T}}(\bm{I}_{p}+o(n^{-1/2}))\bm{\theta}^{0}-\bm{\delta}_{n}^{\mathrm{T}}(\bm{C}^{-1}+o(n^{-1/2}))\bm{\delta}_{n},$

which is $\sigma_{0}^{2}+O_{p}(n^{-1/2})$. This implies the first assertion.

The second assertion follows because
$\bm{Y}^{\mathrm{T}}[\bm{I}_{n}-\bm{X}(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1}\bm{X}]\bm{Y}$
grows at the rate $n$ in $\mathbb{P}_{0}$-probability. ∎

######  Remark 3.

It is well-known that the least square estimator $\hat{\bm{\theta}}$ is
$\sqrt{n}$-consistent for $\bm{\theta}$, and the variance estimator
$\hat{\sigma}^{2}=n^{-1}\|\bm{Y}-\bm{X}\hat{\bm{\theta}}\|^{2}$ is
$\sqrt{n}$-consistent for $\sigma^{2}$. We can write
$\hat{\bm{\theta}}^{\mathrm{R}}=(\bm{C}_{n}+n^{-1}a_{n}\bm{I}_{p})^{-1}\bm{C}_{n}\hat{\bm{\theta}}=(\bm{I}_{p}+\bm{o}(n^{-1/2}))\hat{\bm{\theta}}$.
Hence $\|\hat{\bm{\theta}}^{\mathrm{R}}-\hat{\bm{\theta}}\|=o_{p}(n^{-1/2})$,
and in particular, $\hat{\bm{\theta}}^{\mathrm{R}}$ is $\sqrt{n}$-consistent
for $\bm{\theta}$. Then it follows from the Cauchy-Schwarz inequality and
Assumption 1 that

$n^{-1}\|\bm{Y}-\bm{X}\hat{\bm{\theta}}^{\mathrm{R}}\|^{2}=n^{-1}[\|\bm{Y}-\bm{X}\hat{\bm{\theta}}\|^{2}-2(\bm{Y}-\bm{X}\hat{\bm{\theta}})^{\mathrm{T}}\bm{X}(\hat{\bm{\theta}}^{\mathrm{R}}-\hat{\bm{\theta}})+\|\bm{X}(\hat{\bm{\theta}}^{\mathrm{R}}-\hat{\bm{\theta}})\|^{2}]$

is also $\sqrt{n}$-consistent for $\sigma^{2}$.

######  Lemma 3.

Suppose $\theta_{j}^{0}=0$ for $j=s_{0}+1,\dots,p$. Then for any $\Delta$,
almost surely

$\mathrm{P}(\bm{T}^{*}_{S_{0}^{c}}=\bm{0}_{S_{0}^{c}})|\Delta)>0.$

  * Proof.The minimizer of (2) can be expressed as $\bm{T}^{*}=\begin{bmatrix}\bm{T}_{S_{0}}^{*}\\\ \bm{T}_{S_{0}^{c}}^{*}\end{bmatrix}$. Write

$\bm{C}=\begin{bmatrix}\bm{C}_{11}&\bm{C}_{12}\\\
\bm{C}_{21}&\bm{C}_{22}\end{bmatrix}\textnormal{ and
}\bm{W}^{*}=\begin{bmatrix}\bm{W}_{1}^{*}\\\ \bm{W}_{2}^{*}\end{bmatrix}.$

If $\bm{T}_{S_{0}^{c}}^{*}=\textbf{0}_{S_{0}^{c}}\in\mathbb{R}^{p-s_{0}}$,
then, the KKT conditions of (2) can be written as

$\displaystyle\bm{C}_{11}\bm{T}_{S_{0}}^{*}-(\bm{C}\bm{W}^{*})_{1}=-\frac{\lambda_{0}}{2}\begin{bmatrix}\text{sign}(\theta_{1}^{0})\\\
\vdots\\\ \text{sign}(\theta_{s_{0}}^{0})\end{bmatrix}$ (.2)

and

$\displaystyle-\frac{\lambda_{0}}{2}\bm{1}\leq\bm{C}_{21}\bm{T}_{S_{0}}^{*}-(\bm{C}\bm{W}^{*})_{2}\leq\frac{\lambda_{0}}{2},$
(.3)

where
$(\bm{C}\bm{W}^{*})_{1}=\bm{C}_{11}\bm{W}^{*}_{1}+\bm{C}_{12}\bm{W}^{*}_{2}$,
$(\bm{C}\bm{W}^{*})_{2}=\bm{C}_{21}\bm{W}^{*}_{1}+\bm{C}_{22}\bm{W}^{*}_{2}$
and $\bm{1}$ is a $(p-s_{0})$ dimensional vector of 1’s. Then, solving for
$\bm{T}_{S_{0}}^{*}$ from (.2) and substituting in (.3), we have

$-\frac{\lambda_{0}}{2}\bm{1}\leq\bm{C}_{21}\bm{C}_{11}\big{(}(\bm{C}\bm{W}^{*})_{1}-\frac{\lambda_{0}}{2}\begin{bmatrix}\text{sign}(\theta_{1}^{0})\\\
\vdots\\\
\text{sign}(\theta_{s_{0}}^{0})\end{bmatrix}\big{)}-(\bm{C}\bm{W}^{*})_{2}\leq\frac{\lambda_{0}}{2}\bm{1}.$

Now, recalling that
$\bm{W}^{*}|\bm{\Delta}\sim\mathrm{N}(\sigma_{0}\bm{C}^{-1/2}\bm{\Delta},\sigma_{0}^{2}\bm{C}^{-1})$,
we have the result. ∎

  * Proof of Theorem 2.In view of Lemma 2, it suffices to restrict $\sigma$ to a shrinking neighborhood $\mathcal{U}_{n}$ of $\sigma_{0}$ and derive the result by conditioning on $\sigma$ to a value $\sigma_{n}$, uniformly in $\sigma_{n}\in\mathcal{U}_{n}$. We can write $\bm{\theta}^{*}=\operatorname*{arg\min}_{\bm{u}}M_{n}(\bm{u})$, where

$\displaystyle
M_{n}(\bm{u})=\operatorname*{arg\min}_{\bm{u}}\\{(\bm{u}-\bm{\theta})^{\mathrm{T}}\bm{X}^{\mathrm{T}}\bm{X}(\bm{u}-\bm{\theta})+{\lambda_{0}}{\sqrt{n}}\|\bm{u}\|_{1}\\}$

Then, defining $\bm{t}=\sqrt{n}(\bm{u}-\bm{\theta}^{0})$, we have

$\bm{T}^{*}_{n}=\sqrt{n}({\bm{\theta}}^{*}-\bm{\theta}^{0})=\operatorname*{arg\min}_{\bm{t}}\big{\\{}M_{n}(\bm{\theta}^{0}+n^{-1/2}{\bm{t}})-M_{n}(\bm{\theta}^{0})\big{\\}}.$

Because
$\bm{\theta}|(\bm{Y},\sigma=\sigma_{n})\sim\mathrm{N}(\hat{\bm{\theta}}^{\mathrm{R}},\sigma_{n}^{2}(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1})$,
we can rewrite
$\bm{\theta}=\hat{\bm{\theta}}^{\mathrm{R}}-\sigma_{n}(\bm{X}^{\mathrm{T}}\bm{X}+a_{n}\bm{I}_{p})^{-1/2}\bm{U}_{n}^{*}$,
where $\bm{U}_{n}^{*}\sim\mathrm{N}_{p}(\bm{0},\bm{I}_{p})$. Then, writing
$\bm{\Delta}_{n}=-\sqrt{n}\bm{C}_{n}^{1/2}(\hat{\bm{\theta}}^{\mathrm{R}}-\bm{\theta}^{0})/\sigma_{n}$,
we have

$\displaystyle M_{n}(\bm{\theta}^{0}+n^{-1/2}{\bm{t}})-M_{n}(\bm{\theta}^{0})$
$\displaystyle=$
$\displaystyle(\bm{\theta}^{0}+n^{-1/2}{\bm{t}}-\bm{\theta})^{\mathrm{T}}\bm{X}^{\mathrm{T}}\bm{X}(\bm{\theta}^{0}+n^{-1/2}{\bm{t}}-\bm{\theta})+\lambda_{0}\sqrt{n}\|\bm{\theta}^{0}+n^{-1/2}{\bm{t}}\|_{1}$
$\displaystyle-(\bm{\theta}^{0}-\bm{\theta})^{\mathrm{T}}\bm{X}^{\mathrm{T}}\bm{X}(\bm{\theta}^{0}-\bm{\theta})-\lambda_{0}\sqrt{n}\|\bm{\theta}^{0}\|_{1}$
$\displaystyle=$
$\displaystyle\bm{t}^{\mathrm{T}}\bm{C}_{n}\bm{t}+2\bm{t}^{\mathrm{T}}\bm{C}_{n}\sqrt{n}(\bm{\theta}-\bm{\theta}^{0})+\lambda_{0}\sqrt{n}(\|\bm{\theta}^{0}+n^{-1/2}{\bm{t}}\|_{1}-\|\bm{\theta}^{0}\|_{1})$
$\displaystyle=$
$\displaystyle\bm{t}^{\mathrm{T}}\bm{C}_{n}\bm{t}+2\bm{t}^{\mathrm{T}}\bm{C}_{n}\sqrt{n}(\hat{\bm{\theta}}^{\mathrm{R}}-\sigma_{n}(n\bm{C}_{n}+a_{n}\bm{I})^{-1/2}\bm{U}_{n}^{*}-\bm{\theta}^{0})$
$\displaystyle+\lambda_{0}\sqrt{n}(\|\bm{\theta}^{0}+n^{-1/2}{\bm{t}}\|_{1}-\|\bm{\theta}^{0}\|_{1})$
$\displaystyle=$
$\displaystyle\bm{t}^{\mathrm{T}}\bm{C}_{n}\bm{t}-2\sigma_{n}\bm{t}^{\mathrm{T}}\bm{C}_{n}^{1/2}\bm{\Delta}_{n}-2\sigma_{n}\bm{t}^{\mathrm{T}}\bm{C}_{n}(n\bm{C}_{n}+a_{n}\bm{I})^{-1/2}\bm{U}_{n}^{*}$
$\displaystyle+\lambda_{0}\sqrt{n}(\|\bm{\theta}^{0}+n^{-1/2}{\bm{t}}\|_{1}-\|\bm{\theta}^{0}\|_{1}).$

In view of Lemma 1, $\bm{\Delta}_{n}$ converges weakly to
$\bm{\Delta}\sim\mathrm{N}_{p}(\bm{0},\bm{I}_{p})$. Writing $\bm{U}^{*}$ for a
$\mathrm{N}_{p}(\bm{0},\bm{I}_{p})$ variable independent of $\bm{\Delta}$, it
follows that the weak limit of
$M_{n}(\bm{\theta}^{0}+n^{-1/2}{\bm{v}})-M_{n}(\bm{\theta}^{0})$, conditional
on the data and $\sigma=\sigma_{n}$, is the stochastic process

$\displaystyle\bm{t}^{\mathrm{T}}\bm{C}\bm{t}-2\sigma_{0}\bm{t}^{\mathrm{T}}\bm{C}^{1/2}(\bm{U}^{*}+\bm{\Delta})+\lambda_{0}\big{[}\sum_{j=1}^{s_{0}}t_{j}\mathrm{sign}(\theta^{0}_{j})+\sum_{j=s_{0}+1}^{p}|t_{j}|\big{]}.$
(.4)

The limiting process has a unique separated maximum. Hence, by the Argmax
Theorem (Theorem 3.2.2 of van der Vaart and Wellner [2023]), the conditional
distribution of $\bm{T}_{n}^{*}$ given the data converges to the distribution
of

$\bm{T}^{*}=\operatorname*{arg\min}_{\bm{t}}\\{\bm{t}^{\mathrm{T}}\bm{C}\bm{t}-2\sigma_{0}\bm{t}^{\mathrm{T}}\bm{C}^{1/2}(\bm{U}^{*}+\bm{\Delta})+\lambda_{0}\big{[}\sum_{j=1}^{s_{0}}t_{j}\mathrm{sign}(\theta^{0}_{j})+\sum_{j=s_{0}+1}^{p}|t_{j}|\big{]}\\}$

given $\bm{\Delta}$, in the space of probability measures
$\mathfrak{M}(\mathbb{R}^{p})$ under the weak topology. Writing
$\bm{W}^{*}=\sigma_{0}\bm{C}^{-1/2}(\bm{U}^{*}+\bm{\Delta})$, the
distributional representation in (3.3) is obtained.

Let $\Xi_{n}$ stand for the random probability measure
$\mathcal{L}(\bm{T}_{n}^{*}\in\cdot|\bm{Y},\sigma=\sigma_{n})$. Writing $\Xi$
for the random probability measure $\mathcal{L}(\bm{T}\in\cdot|\bm{\Delta})$,
we have $\Xi_{n}\rightsquigarrow\Xi$ in the space
$\mathfrak{M}(\mathbb{R}^{p})$ with respect to the topology of weak
convergence.

We also note that the conditional distribution of the process
$M_{n}(\bm{\theta}^{0}+n^{-1/2}{\bm{t}})-M_{n}(\bm{\theta}^{0})$ given the
data and $\sigma=\sigma_{n}$ depends on the data through $\bm{\Delta}_{n}$,
which is a continuous function of $\bm{X}^{\mathrm{T}}\bm{\varepsilon}$. Also,
$\sqrt{n}(\bm{\theta}^{\mathrm{L}}-\bm{\theta}_{0})$ is a continuous function
of $\bm{X}^{\mathrm{T}}\bm{\varepsilon}$ too. Further
$n^{-1/2}\bm{X}^{\mathrm{T}}\bm{\varepsilon}\rightsquigarrow\bm{\Delta}$ on
$\mathbb{R}^{p}$. Hence, combining with Theorem 3.1, we can conclude that
$(\Xi_{n},\bm{\xi}_{n})\rightsquigarrow(\Xi,\bm{\xi})$; that is, the joint
weak convergence holds in the space
$\mathfrak{M}(\mathbb{R}^{p})\times\mathbb{R}^{p}$ with respect to the product
of the weak topology and the Euclidean topology. ∎

  * Proof of Theorem 3.Since a norm is a continuous functional, it follows from (3.4) that

$\displaystyle\big{(}\mathcal{L}(\|\bm{T}_{n}^{*}-\hat{\bm{\xi}}_{n}\|_{\mathcal{K}}\in\cdot|\bm{Y}),\|\hat{\bm{\xi}}_{n}\|_{\mathcal{K}}\big{)}\rightsquigarrow\big{(}\mathcal{L}(\|\bm{T}-\bm{\xi}\|_{\mathcal{K}}\in\cdot|\bm{\Delta}),\|\xi\|_{\mathcal{K}}\big{)}$
(.5)

in the space $\mathfrak{M}(\mathbb{R}^{p})\times\mathbb{R}^{p}$. To obtain the
limiting coverage of the obtained credible region, note that
$\\{\bm{\theta}^{0}\in
B_{r}\\}=\\{\|\sqrt{n}(\bm{\theta}^{0}-\hat{\bm{\theta}}^{\mathrm{L}})\|_{\mathcal{K}}\leq
r_{1-\alpha,n}\\}$, and

$\displaystyle\\{\Pi(\|\sqrt{n}(\bm{\theta}^{*}-\hat{\bm{\theta}}^{\mathrm{L}})\|_{\mathcal{K}}\leq\|\sqrt{n}(\bm{\theta}^{0}-\hat{\bm{\theta}}^{\mathrm{L}})\|_{\mathcal{K}}|\bm{Y})<1-\alpha\\}$
$\displaystyle\subset\\{\bm{\theta}^{0}\in B_{r}\\}$
$\displaystyle\subset\\{\Pi(\|\sqrt{n}(\bm{\theta}^{*}-\hat{\bm{\theta}}^{\mathrm{L}})\|_{\mathcal{K}}\leq\|\sqrt{n}(\bm{\theta}^{0}-\hat{\bm{\theta}}^{\mathrm{L}})\|_{\mathcal{K}}|\bm{Y})\leq
1-\alpha\\}.$

The assertion is now immediate. ∎

  * Proof of Theorem 4.Let the limit $\bm{C}$ of $n^{-1}\bm{X}^{\mathrm{T}}\bm{X}$ be a block-diagonal matrix with blocks $\\{j\\}$ and $[-j]$. As argued in Remark 1, it suffices to prove the result under the additional assumption that $c_{j}=1$ and $\sigma_{0}=1$. Then $\xi_{j}$ reduces to $\operatorname*{arg\min}\\{v_{j}^{2}-2\Delta_{j}v_{j}+\lambda_{0}[v_{j}\text{sign}(\theta^{0}_{j})+\mathbbm{1}(\theta_{j}^{0}=0)|v_{j}|]\\}$ For a non-zero-component, the minimizer of $v_{j}^{2}-2\Delta_{j}v_{j}+\lambda_{0}v_{j}\text{sign}(\theta^{0}_{j})$ is

$\displaystyle\xi_{j}=\begin{cases}\Delta_{j}-{\lambda_{0}}/{2},&\text{if
}\theta_{j}^{0}>0,\\\ \Delta_{j}+{\lambda_{0}}/{2},&\text{if
}\theta_{j}^{0}<0,\end{cases}$ (.6)

while for a zero-component, the minimizer of
$v_{j}^{2}-2\Delta_{j}v_{j}+\lambda_{0}|v_{j}|$ is

$\displaystyle\xi_{j}=\begin{cases}\Delta_{j}-{\lambda_{0}}/{2},&\text{if
}\Delta_{j}>{\lambda_{0}}/{2},\\\ \Delta_{j}+{\lambda_{0}}/{2},&\text{if
}\Delta_{j}<-{\lambda_{0}}/{2},\\\ 0,&\text{if
}|\Delta_{j}|\leq{\lambda_{0}}/{2}.\end{cases}$ (.7)

In both cases, we note that $\xi_{j}$ is a continuous function of
$\Delta_{j}$. In the first case, $\xi_{j}$ is a continuous random variable,
while in the second case, $\xi_{j}$ has a point mass at $0$. Proceeding in the
same way, we have for a non-zero component

$\displaystyle
T_{j}^{*}=\operatorname*{arg\min}\\{t_{j}^{2}-2W_{j}^{*}t_{j}+\lambda_{0}t_{j}\textnormal{sign}(\theta_{j}^{0})\\}=\begin{cases}W_{j}^{*}-{\lambda_{0}}/{2},&\text{if
}\theta_{j}^{0}>0,\\\ W_{j}^{*}+{\lambda_{0}}/{2},&\text{if
}\theta_{j}^{0}<0.\end{cases}$ (.8)

and for a zero-component

$\displaystyle
T_{j}^{*}=\operatorname*{arg\min}\\{t_{j}^{2}-2W_{j}^{*}t_{j}+\lambda_{0}|t_{j}|\\}=\begin{cases}W_{j}^{*}-{\lambda_{0}}/{2},&\text{if
}W_{j}^{*}>{\lambda_{0}}/{2},\\\ W_{j}^{*}+{\lambda_{0}}/{2},&\text{if
}W_{j}^{*}<-{\lambda_{0}}/{2},\\\ 0,&\text{if
}|W_{j}^{*}|\leq{\lambda_{0}}/{2},\end{cases}$ (.9)

where $W_{j}^{*}|\Delta_{j}\sim\mathrm{N}(\Delta_{j},1)$. Writing
$Z_{j}=W_{j}^{*}-\Delta_{j}$, and noting that
$Z_{j}|\Delta_{j}\sim\mathrm{N}(0,1)$, we obtain that for $\theta_{j}^{0}>0$,

$\displaystyle\mathrm{P}(|T_{j}^{*}-\xi_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})=\mathrm{P}(|Z_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})=2\Phi(|\xi_{j}|)-1=h_{+}(\lambda_{0},\Delta_{j}),$
(.10)

because
$T_{j}^{*}-\xi_{j}=(W_{j}^{*}-\lambda_{0}/2)-(\Delta_{j}-\lambda_{0}/2)=Z_{j}$
in this case.

Similarly, for $\theta_{j}<0$,
$T_{j}^{*}-\xi_{j}=(W_{j}^{*}+\lambda_{0}/2)-(\Delta_{j}+\lambda_{0}/2)=Z_{j}$
again, and hence

$\displaystyle\mathrm{P}(|T_{j}^{*}-\xi_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})=\mathrm{P}(|Z_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})=2\Phi(|\xi_{j}|)-1=h_{-}(\lambda_{0},\Delta_{j}).$
(.11)

Since $2\Phi(|\xi_{j}|)-1$ is a continuous random variable, the coverage
probability
$\mathbb{P}_{0}(\theta_{j}^{0}\in[\hat{\theta}_{j}^{\mathrm{L}}-r_{1-\alpha,j},\hat{\theta}_{j}^{\mathrm{L}}+r_{1-\alpha,j}])$
converges to

$\displaystyle\mathrm{P}(2\Phi(|\xi_{j}|)-1\leq
1-\alpha)=\mathrm{P}(|\xi_{j}|\leq
z_{\alpha/2})=\Phi(\lambda_{0}/2+z_{\alpha/2})-\Phi(\lambda_{0}/2-z_{\alpha/2})$

whenever $\theta^{0}_{j}\neq 0$.

To compute $\mathrm{P}(|T_{j}^{*}-\xi_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})$ for
$\theta_{j}^{0}=0$, we consider the three cases separately. For
$\Delta_{j}>\lambda_{0}/2$, we have $\xi_{j}=\Delta_{j}-\lambda_{0}/2>0$ so
that from (.9) and (.7),

$\displaystyle\mathrm{P}(|T_{j}^{*}-\xi_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})$
$\displaystyle=\mathrm{P}(|W_{j}^{*}-\lambda_{0}/2-\xi_{j}|\leq\xi_{j},W_{j}^{*}>\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle\quad+\mathrm{P}(|W_{j}^{*}+\lambda_{0}/2-\xi_{j}|\leq\xi_{j},W_{j}^{*}<-\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle\quad+\mathrm{P}(|0-\xi_{j}|\leq\xi_{j},|W_{j}^{*}|\leq\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle=\mathrm{P}(|Z_{j}|\leq\Delta_{j}-\lambda_{0}/2,Z_{j}>-\Delta_{j}+\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle\quad+\mathrm{P}(|Z_{j}+\lambda_{0}|\leq\Delta_{j}-\lambda_{0}/2,Z_{j}<-\Delta_{j}-\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle\quad+\mathrm{P}(|-\Delta_{j}+\lambda_{0}/2|\leq\Delta_{j}-\lambda_{0}/2,|Z_{j}+\Delta_{j}|\leq\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle=\mathrm{P}(-\Delta_{j}+\lambda_{0}/2\leq
Z_{j}\leq\Delta_{j}-\lambda_{0}/2|\Delta_{j})+0$
$\displaystyle\quad+\mathrm{P}(-\Delta_{j}-\lambda_{0}/2\leq
Z_{j}\leq-\Delta_{j}+\lambda_{0}/2|\Delta_{j})$
$\displaystyle=\Phi(\Delta_{j}-\lambda_{0}/2)-\Phi(-\Delta_{j}-\lambda_{0}/2)$

For the case where $\Delta_{j}<-\lambda_{0}/2$, we have $|\xi_{j}|=-\xi_{j}$,
so that

$\displaystyle\mathrm{P}(|T_{j}^{*}-\xi_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})$
$\displaystyle=\mathrm{P}(|Z_{j}-\lambda_{0}|\leq-\Delta_{j}-\lambda_{0}/2,Z_{j}>-\Delta_{j}+\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle\quad+\mathrm{P}(|Z_{j}|\leq-\Delta_{j}-\lambda_{0}/2,Z_{j}<-\Delta_{j}-\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle\quad+\mathrm{P}(|\Delta_{j}+\lambda_{0}/2|\leq|\Delta_{j}+\lambda_{0}/2|,|Z_{j}+\Delta_{j}|\leq\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle=\Phi(-\Delta_{j}+\lambda_{0}/2)-\Phi(\Delta_{j}+\lambda_{0}/2).$

Finally, when $|\Delta_{j}|\leq\lambda_{0}/2$, $\xi_{j}=0$ so that

$\displaystyle\mathrm{P}(|T_{j}^{*}-\xi_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})$
$\displaystyle=\mathrm{P}(T_{j}^{*}=0\big{|}\Delta_{j})$
$\displaystyle=\mathrm{P}(|Z_{j}+\Delta_{j}|\leq\lambda_{0}/2\big{|}\Delta_{j})$
$\displaystyle=\Phi(\Delta_{j}+\lambda_{0}/2)-\Phi(\Delta_{j}-\lambda_{0}/2).$

Thus in all cases for $\theta_{j}^{0}=0$,
$\mathrm{P}(|T_{j}^{*}-\xi_{j}|\leq|\xi_{j}|\big{|}\Delta_{j})=h_{0}(\lambda_{0},\Delta_{j})$,
which is a continuous random variable because $\zeta\mapsto
h_{0}(\lambda_{0},\zeta)$ is continuous and not flat, and
$\Delta_{j}\sim\mathrm{N}(0,1)$ is continuous. Therefore, it follows that

$\displaystyle\mathbb{P}_{0}(\theta_{j}^{0}\in[\hat{\theta}_{j}^{\mathrm{L}}-r_{1-\alpha,j},\hat{\theta}_{j}^{\mathrm{L}}+r_{1-\alpha,j}])\to\mathrm{P}(h_{0}(\lambda_{0},\Delta_{j})\leq
1-\alpha)=\psi_{0}(\alpha,\lambda_{0}),$

when $\theta_{j}^{0}=0$. ∎

##  Appendix A Funding

Research is partially supported in part by ARO grant number 76643MA 2020-0945.

## References

  * Bai et al. [2020] Jincheng Bai, Qifan Song, and Guang Cheng.  Nearly optimal variational inference for high dimensional regression with shrinkage priors.  _arXiv preprint arXiv:2010.12887_ , 2020.
  * Belitser and Ghosal [2020] Eduard Belitser and Subhashis Ghosal.  Empirical Bayes oracle uncertainty quantification for regression.  _The Annals of Statistics_ , 48(6):3113–3137, 2020.
  * Belitser and Ghosal [2024] Eduard Belitser and Subhashis Ghosal.  Bayesian uncertainty quantification and structure detection for multiple change points models.  _Bernoulli_ , 2024.
  * Belitser and Nurushev [2019] Eduard Belitser and Nurzhan Nurushev.  General framework for projection structures.  _arXiv preprint arXiv:1904.01003_ , 2019.
  * Bhaumik and Ghosal [2015] Prithwish Bhaumik and Subhashis Ghosal.  Bayesian two-step estimation in differential equation models.  _Electronic Journal of Statistics_ , 9(2):3124–3154, 2015.
  * Bhaumik and Ghosal [2017] Prithwish Bhaumik and Subhashis Ghosal.  Efficient Bayesian estimation and uncertainty quantification in ordinary differential equation models.  _Bernoulli_ , 23(4B):3537–3570, 2017.
  * Bhaumik et al. [2022] Prithwish Bhaumik, Wenli Shi, and Subhashis Ghosal.  Two-step Bayesian methods for generalized regression driven by partial differential equations.  _Bernoulli_ , 28(3):1625–1647, 2022.
  * Breiman [1995] Leo Breiman.  Better subset regression using the nonnegative garrote.  _Technometrics_ , 37(4):373–384, 1995.
  * Bühlmann and van de Geer [2011] Peter Bühlmann and Sara van de Geer.  _Statistics for High-dimensional Data: Methods, Theory and Applications_.  Springer Science & Business Media, 2011.
  * Cai and Guo [2017] T Tony Cai and Zijian Guo.  Confidence intervals for high-dimensional linear regression: Minimax rates and adaptivity.  _The Annals of Statistics_ , 45(2):615–646, 2017.
  * Candes and Tao [2007] Emmanuel Candes and Terence Tao.  The Dantzig selector: Statistical estimation when p is much larger than n.  _The Annals of Statistics_ , 35(6):2313–2351, 2007.
  * Carvalho et al. [2010] Carlos M Carvalho, Nicholas G Polson, and James G Scott.  The horseshoe estimator for sparse signals.  _Biometrika_ , 97(2):465–480, 2010.
  * Castillo and Nickl [2013] Ismaël Castillo and Richard Nickl.  Nonparametric bernstein–von mises theorems in gaussian white noise.  2013\.
  * Castillo et al. [2015] Ismaël Castillo, Johannes Schmidt-Hieber, and Aad van der Vaart.  Bayesian linear regression with sparse priors.  _The Annals of Statistics_ , 43(5):1986–2018, 2015.
  * Chakraborty and Ghosal [2021a] Moumita Chakraborty and Subhashis Ghosal.  Convergence rates for Bayesian estimation and testing in monotone regression.  _Electronic Journal of Statistics_ , 15(1):3478–3503, 2021a.
  * Chakraborty and Ghosal [2021b] Moumita Chakraborty and Subhashis Ghosal.  Coverage of credible intervals in nonparametric monotone regression.  _The Annals of Statistics_ , 49(2):1011–1028, 2021b.
  * Chakraborty and Ghosal [2022] Moumita Chakraborty and Subhashis Ghosal.  Rates and coverage for monotone densities using projection-posterior.  _Bernoulli_ , 28(2):1093–1119, 2022.
  * Chatterjee and Lahiri [2011] Arindam Chatterjee and Soumendra Nath Lahiri.  Bootstrapping lasso estimators.  _Journal of the American Statistical Association_ , 106(494):608–625, 2011.
  * Cox [1993] Dennis D Cox.  An analysis of bayesian inference for nonparametric regression.  _The Annals of Statistics_ , pages 903–923, 1993.
  * Fan and Li [2001] Jianqing Fan and Runze Li.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _Journal of the American Statistical Association_ , 96(456):1348–1360, 2001.
  * Fu and Knight [2000] Wenjiang Fu and Keith Knight.  Asymptotics for lasso-type estimators.  _The Annals of Statistics_ , 28(5):1356–1378, 2000.
  * George and McCulloch [1993] Edward I George and Robert E McCulloch.  Variable selection via Gibbs sampling.  _Journal of the American Statistical Association_ , 88(423):881–889, 1993.
  * Han and Yang [2019] Wei Han and Yun Yang.  Statistical inference in mean-field variational bayes.  _arXiv preprint arXiv:1911.01525_ , 2019.
  * Huang et al. [2016] Xichen Huang, Jin Wang, and Feng Liang.  A variational algorithm for bayesian variable selection.  _arXiv preprint arXiv:1602.07640_ , 2016.
  * Ishwaran and Rao [2005] Hemant Ishwaran and J Sunil Rao.  Spike and slab variable selection: frequentist and Bayesian strategies.  _The Annals of Statistics_ , 33(2):730–773, 2005.
  * Javanmard and Montanari [2018] Adel Javanmard and Andrea Montanari.  Debiasing the lasso: Optimal sample size for Gaussian designs.  _The Annals of Statistics_ , 46(6A):2593–2622, 2018.
  * Knapik et al. [2011] Bartek T Knapik, Aad W van der Vaart, and J Harry van Zanten.  Bayesian inverse problems with gaussian priors.  2011\.
  * Lin and Dunson [2014] Lizhen Lin and David B Dunson.  Bayesian monotone regression using Gaussian process projection.  _Biometrika_ , 101(2):303–317, 2014.
  * Mitchell and Beauchamp [1988] Toby J Mitchell and John J Beauchamp.  Bayesian variable selection in linear regression.  _Journal of the American Statistical Association_ , 83(404):1023–1032, 1988.
  * Mukherjee and Sen [2022] Sumit Mukherjee and Subhabrata Sen.  Variational inference in high-dimensional linear regression.  _Journal of Machine Learning Research_ , 23(304):1–56, 2022.
  * Nickl and van de Geer [2013] Richard Nickl and Sara van de Geer.  Confidence sets in sparse regression.  _The Annals of Statistics_ , 41(6):2852–2876, 2013.
  * Ormerod et al. [2017] John T Ormerod, Chong You, and Samuel Müller.  A variational bayes approach to variable selection.  2017\.
  * Pal [2024] Samhita Pal.  credint.  https://github.com/SamhitaPal3/credInt, 2024.
  * Ray [2017] Kolyan Ray.  Adaptive bernstein–von mises theorems in gaussian white noise.  2017\.
  * Ray and Szabó [2022] Kolyan Ray and Botond Szabó.  Variational bayes for high-dimensional linear regression with sparse priors.  _Journal of the American Statistical Association_ , 117(539):1270–1281, 2022.
  * Sniekers and van der Vaart [2015] Suzanne Sniekers and Aad van der Vaart.  Adaptive bayesian credible sets in regression with a gaussian process prior.  _Electronic Journal of Statistics_ , 9(2):2475–2527, 2015.
  * Song and Liang [2017] Qifan Song and Faming Liang.  Nearly optimal Bayesian Shrinkage for High Dimensional Regression.  _Science China Mathematics_ , 66(2):409–442, 2017.
  * Szabó et al. [2015] Botond Szabó, Aad W van der Vaart, and JH van Zanten.  Frequentist coverage of adaptive nonparametric bayesian credible sets.  2015\.
  * Tibshirani [1996] Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 58(1):267–288, 1996.
  * van de Geer et al. [2014] Sara van de Geer, Peter Bühlmann, Ya’acov Ritov, and Ruben Dezeure.  On asymptotically optimal confidence regions and tests for high-dimensional models.  _The Annals of Statistics_ , 42(3):1166–1202, 2014.
  * van der Pas et al. [2017] Stéphanie van der Pas, Botond Szabó, and Aad van der Vaart.  Uncertainty quantification for the horseshoe (with discussion).  2017\.
  * van der Vaart and Wellner [2023] A W van der Vaart and Jon A Wellner.  _Weak Convergence and Empirical Processes: With Applications to Statistics_.  Springer, second edition, 2023.
  * Wang and Ghosal [2023a] Kang Wang and Subhashis Ghosal.  Coverage of credible intervals in bayesian multivariate isotonic regression.  _The Annals of Statistics_ , 51:1376–1400, 2023a.
  * Wang and Ghosal [2023b] Kang Wang and Subhashis Ghosal.  Bayesian inference for multivariate monotone densities.  _arXiv preprint arXiv:2306.05202_ , 2023b.
  * Wang and Ghosal [2023c] Kang Wang and Subhashis Ghosal.  Posterior contraction and testing for multivariate isotonic regression.  _Electronic Journal of Statistics_ , 17(1):798–822, 2023c.
  * Yang and Martin [2020] Yue Yang and Ryan Martin.  Variational approximations of empirical bayes posteriors in high-dimensional linear models.  _arXiv preprint arXiv:2007.15930_ , 2020.
  * Yang et al. [2020] Yun Yang, Debdeep Pati, and Anirban Bhattacharya.  $\alpha$-variational inference with statistical guarantees.  _The Annals of Statistics_ , 48(2):886–905, 2020.
  * Yoo and Ghosal [2016] William Weimin Yoo and Subhashis Ghosal.  Supremum norm posterior contraction and credible sets for nonparametric multivariate regression.  2016\.
  * Yuan and Lin [2007] Ming Yuan and Yi Lin.  On the non-negative garrotte estimator.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 69(2):143–161, 2007.
  * Zhang [2010] Cun-Hui Zhang.  Nearly unbiased variable selection under minimax concave penalty.  _The Annals of Statistics_ , 38(2):894–942, 2010.
  * Zhang and Zhang [2014] Cun-Hui Zhang and Stephanie S Zhang.  Confidence intervals for low dimensional parameters in high dimensional linear models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(1):217–242, 2014.
  * Zhang and Gao [2020] Fengshuo Zhang and Chao Gao.  Convergence rates of variational posterior distributions.  _The Annals of Statistics_ , 48(4):2180–2207, 2020.
  * Zou [2006] Hui Zou.  The adaptive lasso and its oracle properties.  _Journal of the American Statistical Association_ , 101(476):1418–1429, 2006.