Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
91,267 characters · 16 sections · 0 citation commands
Nonparametric Empirical Bayes Confidence Intervals
\lstset{ basicstyle=, numbers=left, keywordstyle= \color{ blue!70},commentstyle=\color{red!50!green!50!blue!50} }
\noindentKeywords: Empirical Bayes, Individual heterogeneity, Uncertainty quantification, Nonparametric methods, Ill-posed inverse problems \\[4pt] \noindentJEL classification codes: C11, C14, C13
\thispagestyle{empty}
Characterizing individual heterogeneity has become increasingly prominent in applied economics. Much modern empirical work treats individual effects, rather than only their average, as the primary object of interest.\footnote{See \citeA{walters2024empirical,bonhomme2024estimating} for recent overviews. A growing body of empirical work studies unit-specific parameters, including neighborhood effects (\citeNP{chetty2018impactsI}), workplace heterogeneity (\citeNP{card2013workplace}), teacher value-added (\citeNP{chetty2014measuringI}), manager effects (\citeNP{fenizia2022managers}), geographic variation in health care (\citeNP{finkelstein2016sources}), bank-specific credit-supply effects (\citeNP{amiti2018much}), employer-specific discrimination (\citeNP{kline2022systemic,kline2024discrimination}), police-officer heterogeneity (\citeNP{goncalves2021few}), among many others. } This heterogeneity is often substantively important in its own right, as it describes how these economically relevant effects vary across units. Yet individual effects are unobserved and must be inferred from noisy data, often with heteroskedastic variances. Credible documentation of such heterogeneity therefore requires valid and informative uncertainty quantification for the individual effects themselves.
We consider the problem of constructing confidence intervals for unit-specific effects in a heteroskedastic normal means model. For each unit $i$, let $Y_i$ be a noisy estimate of the unobservable individual effect $\theta_i$ and satisfy $Y_i|\theta_i\sim N(\theta_i,\sigma_i^2)$ with known sampling variance $\sigma_i^2$. The naive $z$-interval $Y_i\pm z_{1-\alpha/2}\sigma_i$ has exact frequentist coverage for $\theta_i$, but it can be quite wide when $\sigma_i^2$ is large.\footnote{Here $z_{\alpha}$ is the $\alpha$-quantile of the standard normal distribution. } Motivated by this concern, \citeA{cox1975prediction} and \citeA{morris1983parametric} introduced parametric empirical Bayes confidence intervals (EBCIs) constructed from the posterior distribution of $\theta_i|Y_i$ under a common Gaussian prior for all $\theta_i$. Under this normal-normal model, the resulting parametric EBCIs borrow strength across units and can be much shorter than the naive $z$-interval (\citeNP{morris1983parametric}). However, when the true prior is non-Gaussian, they can substantially undercover.
\citeA[hereafter AKP]{armstrong2022robust} are the first to propose EBCIs that are robust to the failure of the Gaussian prior assumption. Like Cox-Morris EBCIs, their intervals are centered at linear shrinkage estimators, but use worst-case critical values calibrated to guarantee average coverage over a class of priors characterized by certain moments of $\theta_i$. Because the full prior distribution is nonparametrically identified in the normal means model, this moment-based calibration leaves some available information unused. As a result, their robust EBCIs can be conservative when the true prior is far from the least favorable one in that moment class.
This paper proposes nonparametric empirical Bayes confidence intervals (NP-EBCIs) that fully utilize the entire prior distribution. Following \citeA{cox1975prediction,morris1983parametric}, we construct these intervals from the lower and upper quantiles of the posterior distribution of $\theta_i|Y_i$. Rather than imposing a Gaussian prior, we assume $\theta_i\sim G$ and allow the distribution $G$ to be fully nonparametric. We target the oracle posterior credible intervals under the prior $G$ and construct their feasible counterparts by estimating the posterior quantiles nonparametrically. This contrasts with \citeA{armstrong2022robust}, whose robust EBCIs are calibrated using only moment information about $\theta_i$.
The NP-EBCI construction makes posterior quantiles central to our analysis, and our first contribution is to characterize the minimax optimal rate for estimating them. We establish this rate by deriving a minimax lower bound for the estimation risk and constructing a one-step kernel-based estimator that achieves it. This optimal rate is inherently logarithmic, reflecting the fundamental limits of empirical Bayes posterior quantile estimation. This result highlights a sharp contrast with posterior mean estimation (\citeNP{jiang2009general,brown2009nonparametric}). In the normal means model, recovering the unknown prior $G$ is a severely ill-posed deconvolution problem with a logarithmic minimax rate (\citeNP{carroll1988optimal}).\footnote{ Throughout the paper, we use “prior distribution” and “mixing distribution” interchangeably for $G$. In empirical Bayes language, $G$ is the prior distribution for $\theta_i$, while the term “mixing distribution” reflects that the marginal law of $Y_i$ is obtained by mixing the normal distribution over $G$, see Equation ((ref)). } Although posterior mean and posterior quantiles are both nonlinear functionals of $G$, the posterior mean can circumvent the ill-posedness and be estimated at a nearly parametric rate (\citeNP{zhang1997empirical,zhang2009generalized}).\footnote{In the nonparametric empirical Bayes literature, faster rate results for the posterior mean are typically stated in terms of average mean squared errors for estimating $\theta_i$; see e.g. \citeA{jiang2009general,chen2026empirical}. } We show that posterior quantiles cannot: because they are defined through a non-smooth criterion function, their estimation inherits the severe ill-posedness of nonparametric deconvolution.
We next investigate the conditional coverage properties of the feasible NP-EBCI, which relies on estimated posterior quantiles. By conditional coverage, we mean the probability that the reported interval contains \(\theta_i\) after \(Y_i\) has been observed. This conditional perspective has a long history in statistics, dating back at least to \citeA{cox1958some},\footnote{Related classic discussions of the conditional perspective include \citeA{buehler1959some,wallace1959conditional,birnbaum1962foundations,robinson1979conditional}; see also \citeA{rubin1984bayesianly,robins2000conditioning} for later discussions.} who argued that conditioning is needed so that the uncertainty statement reflects what can be learned from the realized data. Our second contribution is to characterize the minimax optimal rate at which the error in conditional coverage probability vanishes, and to show that the feasible NP-EBCI attains this rate. This optimal rate is again logarithmic, reflecting the fundamental limits of nonparametric empirical Bayes inference when targeting exact conditional coverage. Ultimately, the feasible NP-EBCI achieves asymptotically exact conditional coverage, but its coverage probability converges to the nominal level at this minimax optimal logarithmic rate.
As our third contribution, we show that the feasible NP-EBCI has asymptotically exact marginal coverage, with coverage error that again vanishes at a logarithmic rate. By marginal coverage, we mean the ex-ante probability that the reported interval contains \(\theta_i\) under repeated sampling of both $\theta_i$ and the observed data. This notion is commonly referred to as “empirical Bayes coverage” (\citeNP{morris1983parametric}). Our result differs from the parametric EBCI literature, where the error in marginal coverage can vanish at a faster parametric rate under a correctly specified Gaussian prior.\footnote{See e.g. \citeA{datta2002asymptotic,chatterjee2008parametric,yoshimori2014second}.} However, their theoretical guarantees are vulnerable to misspecification when the true prior is non-Gaussian. Our result also differs from \citeA[ Theorem 4.1]{armstrong2022robust}, who show that the empirical Bayes coverage of their robust EBCI is asymptotically at least the nominal level, thus their actual coverage need not be exact even asymptotically. The feasible NP-EBCI, by contrast, achieves asymptotically exact empirical Bayes coverage, but with coverage error of logarithmic order.
This theoretical contrast motivates a comparison of the finite-sample performance of these EBCIs. In our simulations, all three EBCI procedures deliver average length reductions relative to the naive $z$-interval. AKP robust EBCIs with the second moment often exceed nominal coverage and yield the least length reductions. Cox-Morris parametric EBCIs are typically the shortest, but this efficiency can come at the cost of substantial undercoverage when the true prior is non-Gaussian, especially in low signal-to-noise ratio regimes. The feasible NP-EBCI, despite its severe ill-posedness, remains much closer to the nominal coverage level, even in small samples, while still delivering substantial length reductions.
\paragraph{Related literature.} This paper contributes to the nonparametric empirical Bayes literature on the normal means problem. The central idea in this literature is to treat many related decision problems jointly, learn prior information from related units, and then use it to construct Bayes posterior summaries or decision rules for each unit (\citeNP{robbins1956empirical,robbins1964empirical,efron2014two,efron2019bayes}, see also \citeNP{koenker2024empirical,walters2024empirical} and references therein). Much of this literature has focused on estimation of individual effects,\footnote{See \citeA{chen2026empirical} and references therein. For other settings beyond normal mean models, see e.g. \citeA{gu2017unobserved,liu2020forecasting,gilraine2020new,gaillac2024predicting,kwon2026optimal,cheng2025optimal}. } while more recent work, including \citeA{jiang2019comment} and \citeA{koenker2020empirical}, has considered confidence intervals for individual effects. Both papers provide constructions and numerical evidence, but do not develop theoretical guarantees for their procedures. Our contribution is to establish a general theory for posterior quantile estimation and the coverage properties of posterior-quantile-based EBCIs. More broadly, because posterior quantiles are Bayes rules under asymmetric loss, these results are also relevant for other empirical Bayes decision problems that involve posterior quantiles (e.g. \citeNP{gu2023invidious}, Section 2.5 in \citeNP{walters2024empirical}).
This paper is broadly related to, but different from \citeA{ignatiadis2022confidence} and \citeA{ignatiadis2025empirical}. Specifically, \citeA{ignatiadis2022confidence} study confidence intervals for empirical Bayes estimands, such as the posterior mean, rather than confidence intervals for the individual effect $\theta_i$ itself. \citeA{ignatiadis2025empirical} study empirical partially Bayes multiple testing. As they note in Footnote 2, their $p$-values can be inverted to construct individual intervals for $\theta_i$. A key distinction is that they place a prior on the variances rather than individual effects, pooling information across units through the variance distribution. Consequently, in the homoskedastic model, their interval reduces to the naive $z$-interval (see their Example 14), while our procedure continues to exploit the common prior on the individual effects.
Our objective differs from that of the deconvolution literature. Deconvolution methods aim to recover the latent distribution of unobserved heterogeneity from noisy measurements (see e.g. \citeNP{horowitz2014ill,schennach2020mismeasured} for reviews). Our goal, by contrast, is individualized inference for $\theta_i$ itself. The role of deconvolution in our analysis is therefore technical rather than substantive: the posterior quantiles underlying our NP-EBCIs are functionals of the latent distribution, so the problem is built on the same Gaussian deconvolution structure.
Our results on minimax optimal rates for posterior quantiles connect to the literature on functional estimation in ill-posed inverse problems. For linear functionals of the mixing distribution $G$, the minimax rate can depend sharply on the smoothness of the weighting function that defines the functional (\citeNP{butucea2009adaptive,pensky2017minimax}). Our setting is different, because both posterior mean and posterior quantiles are nonlinear functionals of the mixing distribution. But a similar insight emerges after local linearization: the two functionals differ sharply in their local dependence on the mixing distribution. The posterior mean involves a comparatively smooth weighting function and therefore can circumvent the severe ill-posedness. By contrast, the posterior quantile is locally governed by a discontinuous weighting function and therefore inherits the logarithmic difficulty of the underlying inverse problem.
\paragraph{Outline.} (ref) presents the set-up, defines the oracle NP-EBCI, and contrasts it with the Cox-Morris parametric EBCI and AKP robust EBCIs. Assuming a homoskedastic normal mean model, (ref) derives the minimax optimal rate of estimating posterior quantiles, and establishes the conditional and marginal coverage properties of the feasible NP-EBCI. (ref) returns to the heteroskedastic setting and presents a feasible implementation of the NP-EBCI, including bandwidth selection. (ref) uses Monte Carlo simulations to compare the finite-sample performance of NP-EBCI against existing methods. (ref) concludes. Proofs and extended discussions are collected in the Appendix.
We consider the following normal means model:
where $Y_i$ is a noisy estimate of the unobserved outcome of interest $\theta_i$ for each unit $i$. For example, $\theta_i$ may represent the value-added of teacher $i$, and $Y_i$ is an estimate constructed from the associated students' test scores. The unit-specific parameter $\theta_i$ is thought to be random, and drawn from the population distribution $G(\cdot)$. The variances $\sigma_i^2$ are known, but $G$ and ${\bm \theta}=\{\theta_1,...,\theta_n\}$ are unknown.
The model ((ref)) arises naturally from a hierarchical data structure. For instance, let $\tilde{Y}_{i,j}$ denote the test outcome of student $j=1,\ldots,J_{i}$ taught by teacher $i$, and suppose that these test scores are normally distributed with a teacher-specific mean and variance: $\tilde{Y}_{i,j}|\theta_i\sim N(\theta_i,\sigma_i^2)$ for $j=1,...,J_{i}$ where $\theta_i$ represents the value-added of teacher $i$. A natural estimator for $\theta_i$ is the average test score of students taught by teacher $i$: $Y_i=J_{i}^{-1}\sum_{j=1}^{J_i}\tilde{Y}_{ij}$. Then we have that $Y_i|\theta_i\sim N(\theta_i,\sigma_i^2/J_i)$ which falls into the model ((ref)).
The model ((ref)) fits into the classical measurement errors model, with $Y_i=\theta_i+\varepsilon_i$ and $\varepsilon_i\sim N(0,\sigma_i^2)$. The normal error assumption is restrictive, but standard in the empirical Bayes literature with cross-sectional data (\citeNP{efron2012large,koenker2024empirical,walters2024empirical}). It can be replaced by any other known error distribution $F_{\varepsilon_i}(\cdot)$, or even fully relaxed to an unknown distribution when panel data or other repeated measurements are available.\footnote{See, for example, \citeA{horowitz1996semiparametric} and \citeA{li1998nonparametric}, as well as the surveys by \citeA{schennach2020mismeasured} and \citeA[Section 7]{bonhomme2024estimating}, for further details. }
The mixing distribution $G(\cdot)$ is nonparametrically identified from the marginal law of $Y_i$. Assume that the mixing distribution $G(\cdot)$ has a continuous density $g(\cdot)$ with respect to Lebesgue measure. Given the known variance $\sigma_i^2$, the marginal density of $Y_i$ is given by
where $\phi(\cdot)$ is the standard normal density. For an integrable function $f:\mathbb{R}\rightarrow\mathbb{C}$, write $f^{\star}(t)=\int_{\mathbb{R}}\exp(itx)f(x)\mathrm{d}x$ for its Fourier transform. Let $f_{Y_i}^{\star}(t), g^{\star}(t)$ be the Fourier transforms of the densities $f_{Y_i}(y), g(\theta)$ respectively. The convolution theorem (Section 3.3.2, \citeNP{schennach2020mismeasured}) gives that $f_{Y_i}^{\star}(t)=g^{\star}(t)f_{\varepsilon_i}^{\star}(t)$ where $f_{\varepsilon_i}^{\star}(t)=\exp\left(-\frac{1}{2}\sigma_i^2t^2\right)$ is the Fourier transform of the normal error distribution. Since $f_{\varepsilon_i}^{\star}(t)\neq 0$ for every $t\in\mathbb{R}$, we recover $g^{\star}(t)=f_{Y_i}^{\star}(t)/f_{\varepsilon_i}^{\star}(t)$. Because characteristic functions uniquely determine probability laws, this establishes the identification of $G$.\footnote{This simplified exposition is intended to align with the theoretical analysis in (ref). However, we note that the identification argument applies to any Borel probability measure $G$ (\citeNP[Proposition 8.50]{folland1999real}). }
We aim to construct an individual interval $\mathrm{CI}_i$ for each unobservable individual effect $\theta_i$ with a pre-specified coverage probability. Following \citeA{cox1975prediction,morris1983parametric}, we take an empirical Bayesian approach: assume that $\theta_i\sim G$ for some prior distribution $G(\cdot)$ before any observation is taken. After observing $Y_i=y$, uncertainty about $\theta_i$ is summarized by the posterior distribution $\theta_i\mid Y_i$.
If the true prior $G$ were known, we could construct an oracle Bayesian credible interval by taking the lower and upper $\alpha/2$ posterior quantiles.\footnote{Our interval is the “equal-tailed credible set” for $\theta_i$. One could alternatively construct a highest posterior density (HPD) credible set, defined as the set where the posterior density exceeds some threshold (\citeNP[Definition 5]{Berger1985}). However, a practical limitation of the HPD credible set is that it may be disjoint, making it difficult to report. } Let $\mathcal{D}_i=(Y_i,\sigma_i)$, we denote this infeasible oracle interval by
where the posterior quantile is defined as
The subscript $G$ is used to emphasize that these posterior quantiles are functionals of the underlying prior distribution.
Because $G$ is unknown in practice, we must distinguish between the oracle target ((ref)) and its feasible analogue. Our proposed nonparametric empirical Bayes confidence interval (NP-EBCI) replaces the oracle posterior quantiles with nonparametric estimates:
The theoretical justification of the feasible procedure ((ref)) relies on large-sample asymptotics as $n\rightarrow\infty$, in the spirit of \citeA{robbins1964empirical}. Before turning to that analysis in (ref), we discuss several oracle properties of the interval ((ref)). First, the oracle interval ((ref)) admits a simple decision-theoretic characterization.
(ref) gives a direct decision-theoretic justification for using posterior quantiles as interval endpoints. The interval ((ref)) is not only an equal-tailed posterior credible interval, but also the Bayes rule for an interval-valued decision problem that trades off interval length against linear penalties for undercoverage and overcoverage. In this sense, the lower and upper posterior quantiles arise as optimal endpoints under a particular loss function, in line with the decision-theoretic treatment of Bayesian interval estimation studied by \citeA{winkler1972decision,Berger1985}. The criterion in (ref) also coincides with the interval score used to evaluate prediction intervals (\citeNP{gneiting2007strictly}).
The oracle interval ((ref)) has the usual conditional coverage property of an empirical Bayes confidence interval (EBCI) for $\theta_i$ (\citeNP[Definition 3.2]{carlin2000bayes}):
that summarizes the uncertainty about $\theta_i$ after the data are observed. It is interpreted as a data-specific measure of uncertainty: it evaluates inference procedure based on the realized sample rather than averaging over hypothetical, unobserved outcomes. This conditional perspective goes back at least to \citeA{cox1958some} who wrote that inference should reflect what can be learned from the data that we have, and it is central to Bayesian literature which treats the observed data as known and concerns the remaining uncertainty through the conditional distribution of unknowns given knowns (\citeNP{rubin1984bayesianly}, see also \citeNP[Section 1.6]{Berger1985}).
The oracle interval ((ref)) also satisfies the marginal coverage property, often called empirical Bayes coverage (\citeNP{morris1983parametric}, \citeNP[Definition 3.1]{carlin2000bayes}). To see this, take averages over the observed data distribution:
the probability in $\operatorname{P}\left( \theta_i\in\mathrm{CI}_i \right)$ is taken over both ${\bm \theta}=(\theta_1,...,\theta_n)$ and the data $\mathbf{Y}=(Y_1,...,Y_n)$.\footnote{This notion of marginal coverage also appears in the conformal prediction literature, see e.g. \citeA{angelopoulos2023gentle}. } Thus marginal coverage is an ex ante statement: it characterizes the uncertainty prior to observing any data. It can also be interpreted as the population average of frequentist coverage $\operatorname{P}(\theta_i\in\mathrm{CI}_i^{\mathrm{NP}*}| \bm \theta)$ that considers the uncertainty from the data. This marginal coverage statement is weaker than the conditional coverage property ((ref)) because it need only hold on average and therefore may be poorly calibrated for the sample data at hand. Even so, the notion of marginal coverage is standard in the empirical Bayes literature.
We now compare the oracle NP-EBCI with the existing alternatives in the literature. The natural starting point is the naive $z$-interval $Y_i\pm z_{1-\alpha/2}\sigma_i$, which treats each unit in isolation and has exact frequentist coverage for $\theta_i$ without imposing any prior structure. The limitation is that when $\sigma_i^2$ is large, the interval can be unnecessarily wide. This concern leads naturally to the empirical Bayes confidence interval (EBCI) literature, which aims to improve inference for each unit by combining the noisy observation $Y_i$ with prior information learned from the distribution of individual effects.
\citeA{cox1975prediction} and \citeA{morris1983parametric} initiated the parametric EBCI literature by constructing individual intervals for $\theta_i$ from the posterior distribution under a Gaussian prior. In the homoskedastic case $\sigma_i^2\equiv \sigma^2$, they assume $\theta_i\stackrel{\rm i.i.d.}{\sim} N(0,A)$. Under this normal-normal specification, the posterior distribution takes a simple form: $\theta_i|Y_i=y\sim N(\frac{A}{A+\sigma^2}y,\frac{A\sigma^2}{A+\sigma^2})$ which yields the oracle interval
A feasible version of this interval replaces the unknown variance component $A$ with an estimator. This construction can produce substantially shorter intervals than the naive $z$-interval. This gain, however, comes from the parametric prior restriction.
Like \citeA{cox1975prediction} and \citeA{morris1983parametric}, we construct intervals from the posterior distribution of $\theta_i|Y_i$. The difference is that we make no parametric assumptions about the prior distribution. Instead we treat the prior $G$ as fully nonparametric and take the oracle posterior interval in ((ref)) as the target. The interval proposed in practice is the feasible analog in ((ref)), obtained by replacing the oracle posterior quantiles with nonparametric estimates.
A different alternative is the robust EBCI of \citeA{armstrong2022robust}, which targets worst-case marginal coverage over a moment class rather than the posterior interval under a point-identified prior. In the homoskedastic case $\sigma_i^2\equiv \sigma^2$, their oracle interval takes the form
where $A=E[Y_i^2]-\sigma^2$. The critical value $\chi(\alpha)$ is chosen such that the optimal value of the following optimization problem is at least $1-\alpha$:
By construction, the marginal coverage probability of their oracle robust EBCI is at least $1-\alpha$ for all prior distributions that satisfy the second moment constraint. However, the actual coverage probability may exceed the nominal level when the true prior is not least-favorable distribution that solves ((ref)). The oracle NP-EBCI, by contrast, relies on nonparametric identification of the prior distribution and fully exploits the entire distribution of $\theta_i$. Consequently, the interval ((ref)) attains exact marginal coverage at the nominal level $1-\alpha$. Moreover, the AKP robust EBCI does not target the exact posterior conditional coverage property in ((ref)) which conditions on the observed data. By contrast, the oracle NP-EBCI ((ref)) satisfies this property by construction.
We make two additional remarks to clarify how certain terminology in empirical Bayes inference are understood in our setting.
In this section, we analyze the theoretical properties of the feasible NP-EBCI. First we consider nonparametric estimation of the posterior quantile: we establish its optimal rates of convergence in (ref) and propose a rate-optimal estimator in (ref). Second, in (ref), we establish the optimal rates of the errors in conditional coverage probability and show that the feasible NP-EBCI based on plug-in estimate of posterior quantile achieves this rate. Finally, we examine its marginal coverage probability and connects it to the empirical Bayes literature.
\paragraph{Notation.} For \( 1 \leq p < \infty \), the space \( L^p(\mathbb{R}) \) is defined as \( L^p(\mathbb{R}) = \{ f : \mathbb{R} \to \mathbb{C} \mid \int_{\mathbb{R}} |f(x)|^p \, \mathrm{d}x < \infty \} \), with associated norm \( \|f\|_p = \left( \int_{\mathbb{R}} |f(x)|^p \, \mathrm{d}x \right)^{1/p} \). For \( p = \infty \), we define the norm as \( \|f\|_\infty = \operatorname{ess\,sup}_{x \in \mathbb{R}} |f(x)| \). For any function \( f : \mathbb{R} \to \mathbb{C} \) such that \( f \in L^1(\mathbb{R}) \), we define its Fourier transform by $f^{\star}(t) = \int_{\mathbb{R}} \exp(itx) f(x) \, \mathrm{d}x$ for \( t \in \mathbb{R} \).
Throughout (ref), we assume that $\sigma_i^2\equiv \sigma^2$ for technical convenience, which implies that $(Y_i,\theta_i)\stackrel{\rm i.i.d.}{\sim} \mathcal{F}$. In (ref), we consider the posterior quantile associated with a hypothetical random variable $(Y_{n+1},\theta_{n+1})\sim \mathcal{F}$, evaluated at $Y_{n+1}=y$, and estimate it using the data $\{Y_i\}_{i=1}^n$. In (ref), since our goal is to construct individual intervals $\mathrm{CI}_i$ for each $\theta_i$, the relevant posterior quantile corresponds to $(Y_i,\theta_i)$ and is estimated using the leave-one-out data $\{Y_j\}_{j\neq i}$. We focus on a fixed quantile level $\tau\in(0,1)$, tailored to the objective of the application.
Both posterior mean\footnote{The posterior mean, defined as $m(y)=\operatorname{E}_G[\theta|Y=y]$, is one of the most widely studied objects in empirical Bayes literature. It is the optimal predictor of the individual parameter $\theta_i$ under the mean-square error, conditional on the observed data $Y=y$. } and posterior quantiles are nonlinear integral functionals of the unknown mixing distribution $G$. In the normal means model, the observed distribution of $Y$ is the convolution of $G$ with the normal error distribution, so recovering $G$ from the data is a deconvolution problem. With normal errors, this inverse problem is severely ill-posed: the optimal rate for estimating the mixing density $g$ is only logarithmic (\citeNP{carroll1988optimal}, see also \citeA{horowitz2014ill} for a review).\footnote{ The normal density is a canonical example of a “super-smooth” density, which leads to the worst-case rates of convergence in nonparametric deconvolution. When the error density is “ordinarily smooth”, the optimal rate is polynomial (\citeNP{fan1991optimal}); see \citeA[page 511]{schennach2020mismeasured} for further discussion. } In contrast to this severe ill-posedness, the posterior mean can be estimated at a nearly parametric rate (\citeNP{zhang1997empirical,zhang2009generalized}). Since the posterior quantile is also an integral functional of the unknown mixing density, it is natural to ask whether such “integration” can likewise lead to faster rates of convergence.
We will show that the answer is negative. Unlike the posterior mean, the posterior quantile is defined through a non-smooth criterion function. This non-smoothness prevents the associated Fourier weighting term from offsetting the exponential decay induced by the Gaussian error density. Thus, posterior quantile estimation inherits the severe ill-posedness of Gaussian deconvolution. To show that this difficulty is intrinsic, rather than a deficiency of a particular estimator, we derive a minimax lower bound on the risk of posterior quantile estimation over a functional class. This minimax result characterizes the fundamental limits of the estimation problem and provides a benchmark for evaluating the quality of an estimation procedure.
We begin by rewriting the posterior quantile in a form that exposes its non-smoothness. Although defined in ((ref)) as the solution to a moment condition weighted by the normal error density, it is more useful to reframe it as the minimizer of a convex optimization problem.
We transform the objective function ((ref)) into the frequency domain using Fourier basis functions. Define the weighting function $M(q,\theta;y)=\rho_{\tau}\left( \theta-q \right)\phi(\frac{y-\theta}{\sigma})$ that appears in the integral of the objective function ((ref)). By the Parseval-Plancherel theorem (Theorem 9.13 in \citeA{rudin1987}), the objective function $W(q)$ in (ref) can be expressed as
where $\overline{M^{\star}(q,t;y)}=\int \exp(-it\theta)M(q,\theta;y)\mathrm{d}\theta$ is the conjugate Fourier transform of $M(q,\theta;y)$.
The Fourier representation ((ref)) shows that the problem has two ingredients. First, Gaussian deconvolution is severely ill-posed because $f_{\varepsilon}^{\star}(t)$ decays exponentially as $|t|\rightarrow \infty$, so division by $f_{\varepsilon}^{\star}(t)$ amplifies high-frequency noise. Second, whether this difficulty can be mitigated depends on the decay of the Fourier weighting term. For smooth functionals, such as the posterior mean, that weighting term decays fast enough to downweigh the exponential decay at the tail, leading to a nearly parametric rate.\footnote{ A similar insight appears in the semiparametric models with measurement errors that achieve the parametric $\sqrt{n}$ rate, see \citeA[Page 512]{schennach2020mismeasured} for related discussion. } For posterior quantiles, however, the check function in ((ref)) is non-smooth at zero, so $\overline{M^{\star}(q,t;y)}$ decays too slowly to offset it. The posterior quantile therefore inherits the severe ill-posedness of Gaussian deconvolution.
To state the formal results, we introduce the following regularity assumptions.
We denote the posterior density evaluated at $\theta = u$ by
(ref) is a regularity condition that facilitates the application of Fourier-based methods. (ref) is equivalent to requiring that the sum of the squared $L^2(\mathbb{R})$-norms of the first $s$ derivatives of the mixing density $g(\cdot)$ be bounded, i.e. $\sum_{l=0}^{s}\|g^{(l)}\|_2^2<L$, see e.g. \citeA[Appendix A.2]{meister2009deconvolution}. This Sobolev condition provides global control of the smoothness and is well suited for analyzing posterior quantile, a nonlinear integral functional of mixing density, while H\"{o}lder condition imposes smoothness only locally around some specific point.
(ref) is the identification condition for the posterior quantile by requiring that the criterion function ((ref)) is non-flat at its solution, ensuring a unique minimizer. (ref) together with (ref), ensures that the inverse Fourier transform of $g^{\star}(t)$ is well-defined and yields a continuous mixing density $g(\cdot)$. Consequently, the posterior distribution has no atoms at the target quantile. Finally, (ref) ensures that the posterior quantile is well-defined and lies in a compact subset of $\mathbb{R}$.
(ref) establishes a minimax lower bound for estimating posterior quantiles uniformly over the Sobolev space $\mathcal{G}(s,L)$. In particular, it shows that no estimator can attain a faster rate of convergence than logarithmic uniformly over this function class. This lower bound reflects that posterior quantile estimation, due to the non-smoothness of its criterion function, inherits the intrinsic difficulty of nonparametric deconvolution. The next subsection shows that this lower bound is attained by a feasible estimator. Hence the logarithmic rate is the minimax optimal rate for estimating posterior quantiles in the homoskedastic normal-means model.
(ref), together with the nearly parametric rate result for posterior mean estimation (\citeNP{zhang1997empirical,zhang2009generalized}), highlights a broader point in functional estimation for ill-posed inverse problems: the minimax optimal rate is determined jointly by the degree of ill-posedness of the inverse operator and the regularity of the target functional. Posterior quantiles are too irregular to circumvent the difficulty of the underlying severely ill-posed deconvolution problem, and their minimax rate therefore remains logarithmic.
In this subsection, we propose a one-step estimator for the posterior quantile that achieves the optimal rate established in (ref).
Observe that the representation ((ref)) expresses the posterior quantile directly as a function of the data, and naturally suggests an estimator for the posterior quantile, defined as the minimizer of the empirical counterpart of the objective function. In this representation, the only unknown component is $f_{Y}^{\star}(t)$. Therefore, to construct the empirical analog $W_n(q)$, we replace $f_Y^{\star}(t)$ with $\hat{f}^{\star}_Y(t)K^{\star}(h_nt)$, where $\hat{f}^{\star}_Y(t)=n^{-1}\sum_{i=1}^n \exp(it Y_i)$ is the empirical characteristic function, and $K^{\star}(t)$ denotes the Fourier transformation of a kernel function $K(\cdot)$ with bandwidth $h_n$. The kernel serves as a regularization device to ensure the existence of the integral when $f_Y^{\star}(t)$ is replaced by its empirical analog. The resulting kernel estimator of the posterior quantile is then
where the sample objective function is denoted by $W_n(q)$.
Our proposed estimator ((ref)) is a one-step procedure that bypasses the need to estimate the mixing distribution $G(\cdot)$ as an intermediate step. This estimation strategy is in the spirit of $f$-modeling for posterior mean estimation, which involves estimating only the marginal density $f_Y(y)$ and its derivative (\citeNP{robbins1956empirical,efron2011tweedie,efron2014two}). Another advantage of this direct estimation is that the regularization parameter can be tuned directly for the target functional, rather than for an intermediate estimate of the mixing distribution.
Our estimator for the posterior quantile is, by construction, similar to the classical deconvolution kernel density estimator (\citeNP{carroll1988optimal,stefanski1990deconvolving}). However, as discussed in (ref), the posterior quantile is a nonlinear integral functional of the deconvolution density. This leads to two key distinctions and challenges of the theoretical analysis. First, our estimator is defined as the solution to a Fourier integral, rather than an integral itself. Second, the integrand involves a weighting function $\overline{M^{\star}(q,t;y)}$ whose (non-)smoothness crucially affects the convergence rate of the estimator.
We then formally establish the theoretical properties of the proposed estimator $\hat{q}_{G}(\tau;y)$.
(ref) (i) requires that the Fourier transform of the kernel function be bounded, integrable and have compact support. (ref) (ii) requires $K$ to be higher-order kernel, that is a device designed to increase the rate of convergence of the asymptotic bias. Here we use a kernel of the order $s+1$, which is one degree higher than the smoothness of the mixing density.
(ref) establishes that our proposed kernel estimator for the posterior quantile is consistent and achieves the optimal rate specified in (ref). The slow rate of convergence reflects the difficulty of the estimation problem, and is not an indicator that our estimator is deficient.
The key step in the proof of (ref) is to derive the following asymptotic expansion
where $r_n$ denotes the rate of convergence, $\hat{q}=\hat{q}_G(\tau;y)$ is the estimated posterior quantile, and $q_0=q_{G}(\tau;y)$ its population counterpart. This expansion is analogous to the classical Bahadur-Kiefer representation of the sample quantiles (\citeNP{bahadur1966note,kiefer1967bahadur}), and more generally to the linearization of M-estimators (see e.g. Section 3.2.4 in \citeNP{van1996weak}).
However, establishing the asymptotic expansion in (ref) is delicate when the functional of interest is possibly irregular, as is often the case in ill-posed inverse problems. One can invoke stochastic equicontinuity and empirical process arguments to control the nonlinear remainder term,\footnote{See Assumption 3.5 in \citeA{chen2015sieve} for a formal statement of the required high-level conditions. } but this approach has a direct implication on the rates of the functional estimation.\footnote{The issue arises for regular functionals in severely ill-posed inverse problems. Consider a two-step procedure for estimating the posterior mean, where a deconvolution kernel density estimator is used as a plug-in. In this context, stochastic equicontinuity arguments typically bound the nonlinear remainder term by requiring uniform control over the first-stage estimate, i.e. $\|\hat{g}-g\|_{\infty}$, which results in a suboptimal convergence rate for the posterior mean itself. } We instead exploit the convexity of the sample criterion in ((ref)) and apply the convexity lemma of \citeA{pollard1991asymptotics,hjort2011asymptotics}, thereby avoiding a separate stochastic equicontinuity argument. The key step is to establish a uniform quadratic approximation of the localized criterion on compact sets; the convexity lemma then implies that the sample argmin is close to the argmin of this approximation. Thus, the proof relies on the local uniform control already established for the deconvolution kernel density estimator, rather than additional high-level empirical process conditions.
Having established the asymptotic expansion ((ref)), the rate of convergence is governed by the leading term $\partial W_n(q_0)/\partial q$. The same Fourier-space heuristic described after ((ref)) applies here: once the criterion function is linearized, this term behaves like a functional with a non-smooth weighting function. Its Fourier coefficient therefore decays too slowly to offset the exponential decay of $f_{\varepsilon}^{\star}(t)$ as $|t|\rightarrow\infty$, resulting in a logarithmic convergence rate. For comparison, Appendix (ref) studies the analogous kernel-based estimator for the posterior mean. There the leading term is associated with a smooth weighting function whose Fourier coefficient does offset this decay, so with a suitable bandwidth choice the estimator achieves a nearly parametric rate.
We now examine the conditional coverage of the feasible NP-EBCI. Because the target is the conditional coverage probability $\operatorname{P}_{G}(\theta_i\in\mathrm{CI}_i \mid Y_i=y)$, we evaluate $\mathrm{CI}_i$ at the realized observation $Y_i$, while estimating its endpoints from the remaining sample $\mathbf{Y}_{-i}:=\{Y_j\}_{j\neq i}$. This leave-one-out construction separates the observation being conditioned on from the data used to learn the prior.
Let $\hat{q}_G(\tau;Y_i)$ denote the leave-one-out analogue of the posterior quantile estimator in (ref). We define the feasible NP-EBCI as
To evaluate accuracy, we consider the error in conditional coverage probability, that is the difference between the actual and nominal coverage probabilities. The next theorem characterizes the best possible rate at which this error can vanish uniformly over the function class $\mathcal{G}(s,L)$.
(ref) provides a minimax lower bound on the errors in conditional coverage probabilities. The infimum is taken over all two-sided intervals constructed from the data, and the supremum is over all mixing densities satisfying (ref). Thus, if conditional coverage is the target criterion, no leave-one-out two-sided interval can have a coverage error that vanishes faster than a logarithmic rate. In this sense, (ref) reveals the fundamental limits on the accuracy of nonparametric empirical Bayes inference if the conditional coverage is the target criterion.
The proof of (ref) relies on Fano's inequality (\citeNP{scarlett2019introductory}). Unlike the lower bound construction for mean square error, where alternative data distributions are separated by the induced target parameters, the separation here is based on the conditional coverage probabilities. We construct many alternative mixing densities around a baseline such that the induced observed laws of $\mathbf{Y}_{-i}$ are statistically indistinguishable, while the conditional coverage errors are bounded away from zero at the desired rate. Hence, if an interval had conditional coverage error uniformly smaller than this rate, it could decode which alternative generated the data. But Fano's inequality precludes this possibility because the alternatives are statistically close.
We next show that the feasible NP-EBCI ((ref)) achieves the rate in the lower bound established in (ref). Because we construct the feasible interval by replacing infeasible posterior quantiles with their estimators, our task is to translate the theoretical results of posterior quantile estimation ((ref)) into a statement about conditional coverage. To do so, we impose one additional regularity condition on the posterior density at the relevant quantiles.
(ref) shows that the conditional coverage probability of the feasible NP-EBCI ((ref)) asymptotically fully adapts to the nominal $1-\alpha$, but differs at a logarithmic rate. The logarithmic conditional coverage error comes from that the endpoints themselves can only be estimated at a logarithmic rate. Nevertheless, (ref) together establish the optimal rate at which the worst-case conditional coverage probability converges to the nominal level, and the feasible NP-EBCI ((ref)) attains this rate.
We turn to the analysis of the marginal coverage probability of the feasible NP-EBCI ((ref)). As this coverage statement is defined with respect to the joint distribution of both $\theta$ and $\{Y_i\}_{i=1}^n$, our analysis requires additional regularity conditions that hold uniformly for $y\in\mathcal{Y}\subset\mathbb{R}$.
(ref) is a standard condition used to establish $L_2$ or uniform convergence in the quantile regression literature (see, e.g., Assumption 3 in \citeA{horowitz2005nonparametric}, Condition D1 in \citeA{belloni2011}). (ref) is weaker than the uniform version of (ref), requiring only that the expectation be finite.
For the oracle NP-EBCI, marginal coverage follows immediately from conditional coverage by iterated expectations. (ref), however, concerns the feasible NP-EBCI whose endpoints are estimated. To obtain the marginal coverage result for this feasible interval, we need to control the errors in conditional coverage uniformly over $y\in\mathcal{Y}$ so they can be averaged over the distribution of $Y_i$. The averaging does not improve the convergence rate, since the dominant error still comes from estimating posterior quantiles, which suffers from the severe ill-posedness. Hence the errors in the coverage probability of the feasible NP-EBCI ((ref)) vanish at the same logarithmic rate under both marginal and conditional coverage.
(ref) gives the accuracy assessment about the feasible NP-EBCI under a fully nonparametric prior. The resulting slow logarithmic rate again arises from the intrinsic difficulty of nonparametric estimation in ill-posed inverse problems. It is useful to compare this result with the two alternatives discussed above: Cox-Morris parametric EBCIs and robust EBCIs of \citeA{armstrong2022robust}.
We first compare (ref) with the marginal coverage result for Cox-Morris parametric EBCIs. Under a correctly specified Gaussian prior, Cox-Morris EBCIs can achieve a parametric rate for errors in marginal coverage probabilities; see e.g. Theorem 1 in \citeA{yoshimori2014second}. This advantage, however, relies on correct specification of the normal-normal model. When the true prior is non-Gaussian, especially in low signal-to-noise ratio regimes, Cox-Morris EBCIs can substantially undercover. The feasible NP-EBCI avoids this source of parametric misspecification by leaving the prior fully nonparametric. The price is that its marginal coverage error vanishes only at a logarithmic rate.
We then compare our result with AKP robust EBCIs. The oracle NP-EBCI targets exact marginal coverage under the true, point-identified prior. AKP instead calibrate their robust EBCI to guarantee marginal coverage over a moment class of priors, and their Theorem 4.1 implies that, under exchangeability, the feasible interval satisfies $\lim\inf_{n\rightarrow\infty} \operatorname{P}(\theta_i\in \mathrm{CI}_i^{\mathrm{AKP}})\geq 1-\alpha $ (\citeNP{armstrong2022robust}). Because this calibration relies only on finite-dimensional moments, rather than on the full prior distribution, the statistical problem is simpler than posterior quantile estimation. Furthermore, their Theorem 4.1 requires only the convergence of the relevant moments, while our (ref) imposes a stronger continuity condition on the prior distribution. The gain of AKP is robustness over the moment class; the cost is conservatism, since their EBCIs need not be asymptotically exact. The NP-EBCI instead targets exact marginal coverage under the true point-identified prior, but pays the logarithmic price of the severely ill-posed inverse problem.
Taken together, (ref) clarifies the trade-offs behind the feasible NP-EBCI. Relative to the Cox-Morris parametric EBCI, it sacrifices the faster coverage error rates achievable under a correctly specified Gaussian prior in exchange for robustness to prior misspecification. Relative to the AKP robust EBCI, it gives up worst-case moment-class robustness in exchange for asymptotically exact marginal coverage under a point-identified prior. (ref) examines how these trade-offs play out in finite samples.
This section returns to the baseline heteroskedastic normal mean model in ((ref)) and describes how to implement the proposed NP-EBCI when units have known but unequal variances $\sigma_i^2$. (ref) imposed the homoskedastic restriction $\sigma_i^2\equiv \sigma^2$ to focus on the optimal rate results. In empirical applications, however, heteroskedasticity is the natural default. Accordingly, in this section, we outline the practical implementation of our approach and discuss bandwidth selection.
We now describe the baseline implementation of the proposed NP-EBCI under heteroskedastic Gaussian noise. Let $\mathcal{D}_i=(Y_i,\sigma_i)$. The oracle posterior quantile can be written as
where $\overline{M^{\star}(q,t;\mathcal{D}_i)}=\int \exp(-it\theta)\rho_{\tau}(\theta-q)\phi\left( (Y_i-\theta)/\sigma_i\right)\mathrm{d}\theta.$ The Fourier transform of the mixing density $g^{\star}(t)$ can be identified from the following population identity
Here, the leave-one-out form is used for estimation rather than identification: it allows us to estimate the posterior quantile from the other units, while $\mathcal{D}_i$ enters only through the evaluation step.
The representation ((ref)) leads to a natural plug-in estimator. Define
For a given bandwidth $h_n$, we estimate the posterior quantile by
As in the homoskedastic case, the kernel $K^{\star}(h_nt)$ regularizes the Fourier integral. In practice, we use the compactly supported flat-top spectral kernel (\citeNP{politis1999multivariate}) but incorporate a cosine taper (\citeNP{harris1978use}) to smooth the transition from the flat-top region to zero and avoid a hard spectral cutoff. We discuss the choice of the bandwidth parameter $h_n$ in (ref).
Using the estimated posterior quantiles, we construct the NP-EBCI
So far, we have presented the baseline implementation for heteroskedastic variances treated as fixed and known. The same construction also covers the case when $\sigma_i^2$ is random across units and independent with $\theta_i$. To see this, under $\theta_i\perp \!\!\! \perp \sigma_i$, we have $\operatorname{E}[\exp(it Y_i)|\sigma_i]=g^{\star}(t)\exp\left(-\frac{1}{2}\sigma_i^2t^2\right)$, so averaging over the realized units gives the same population identity in ((ref)). In this sense, the baseline implementation can also be viewed as conditioning on observed heteroskedastic variances. Appendix (ref) discusses the extension to the case in which $\theta_i$ and $\sigma_i^2$ may be correlated.
We choose a common bandwidth $h_n $ to minimize average interval length subject to an average coverage requirement.\footnote{One could, in principle, choose a separate bandwidth for each $\mathrm{CI}_i$ by minimizing the error of conditional coverage probability of that interval. This selection method is natural but computationally impractical because it requires solving $n$ distinct bandwidth selection problems.} Since this criterion depends on the unknown mixing distribution $G$, we implement it by a plug-in rule that replaces $G$ with its nonparametric maximum likelihood estimator (NPMLE) $\hat{G}(\cdot)$ (\citeNP{soloff2025multivariate}). In this respect, our bandwidth selection method is close in spirit to plug-in or rule-of-thumb bandwidth choice in nonparametric kernel density estimation (see e.g. \citeNP[Section 2.7]{li2007nonparametric}). However, unlike these methods relying on preliminary pilot values, our procedure uses the plug-in NPMLE $\hat{G}$, which itself is a consistent estimator of the mixing distribution $G$ (\citeNP{kiefer1956consistency,pfanzagl1988consistency}) and requires no additional smoothing tuning parameter (\citeNP{polyanskiy2020self}).
It is important to distinguish the role of the NPMLE in our implementation from the two-step plug-in approach discussed in (ref). Our feasible NP-EBCI is based on the one-step kernel estimator for posterior quantiles. The NPMLE is only used for bandwidth selection. We do not suggest estimating posterior quantiles with the NPMLE plug-in. In simulations not reported here, the NP-EBCI based on these plug-in estimates can substantially undercover, which aligns with the simulation evidence in \citeA{koenker2020empirical}. One plausible explanation is that the NPMLE has discrete support (\citeNP{lindsay1983geometry}), so plugging it directly into the posterior quantile can yield an overly concentrated posterior distribution and hence intervals that are too short.
The bandwidth selection procedure is detailed as follows:
\noindentStep 1. Partition the units $\{1,\ldots,n\}$ into $V$ disjoint folds $\mathcal{I}_1,\dots,\mathcal{I}_V$, and let $\mathcal{I}_v^{c}=\{1,\ldots,n\}\setminus \mathcal I_v$. For each fold $v=1,...,V$, we compute the leave-$v$-out nonparametric maximum likelihood estimator for the mixing distribution $\hat{G}^{(-v)}(\cdot)$.
\noindentStep 2. Let $\mathcal{H}_n$ be a grid of candidate bandwidths. For each $h\in\mathcal{H}_n$, each fold $v=1,\dots,V$ and each held-out unit $i\in\mathcal{I}_v$, we estimate the lower and upper endpoints $\hat{q}_h^{(-v)}(\alpha/2;\mathcal{D}_i)$ and $\hat{q}_h^{(-v)}(1-\alpha/2;\mathcal{D}_i)$ on the fold $\mathcal{I}_v^c$, but evaluate at the held-out unit $\mathcal{D}_i$.
\noindentStep 3. For each $h\in\mathcal{H}_n$, we compute the average conditional coverage and average length of intervals under the out-of-fold prior estimate $\hat{G}^{(-v)}(\cdot)$:
\noindentStep 4. We set the bandwidth that yields the shortest average-length interval among those that satisfy the nominal coverage constraint: \[ \hat h = \underset{h\in\mathcal H_n}{\arg\min}\ \widehat L(h) \qquad\text{s.t.}\qquad \widehat C(h)\ge 1-\alpha. \] Appendix (ref) provides further strategies that can substantially reduce the computational cost of this bandwidth selection procedure.
Both NP-EBCI and the AKP robust EBCI move beyond the parametric Gaussian prior, but use the information about the prior distribution in different ways. NP-EBCI targets the oracle posterior interval under a point-identified nonparametric prior and thus exploits the entire distribution, but suffers from the severe ill-posedness. The AKP robust EBCI, by contrast, calibrate intervals to guarantee worst-case coverage uniform over a moment class of priors, so the resulting intervals may be conservative when the true prior is not least favorable. This section uses Monte Carlo simulations to examine how these differences appear in finite samples. We compare the feasible NP-EBCI with AKP (using the second moment) and the Cox-Morris parametric EBCI, and summarize performance by average marginal coverage and average length reduction relative to the naive $z$-interval.
We generate simulated data from the heteroskedastic normal means model ((ref)) with known variances:
We consider sample sizes $n\in\{100,200,500,1000\}$ and target signal-to-noise ratios
For each design cell $(G,n,\mathrm{SNR})$, we first draw $\sigma_i^2 \stackrel{\rm i.i.d.}{\sim} \mathrm{Lognormal}(0,1)$, $i=1,...,n$, and rescale it to attain the target SNR, and then keep the resulting variance vector fixed across all Monte Carlo replications. We then generate 500 independent replications by redrawing only $\{Y_i,\theta_i\}_{i=1}^n$.
We consider eight designs for the mixing distribution $G$:
\noindentDesign 1. Gaussian: $\theta_i\sim N(0,1)$
\noindentDesign 2. Laplace: $\theta_i\sim \text{Laplace}(0,1/\sqrt{2})$
\noindentDesign 3. Student-$t$: $\theta_i\sim \sqrt{0.2}t_{2.5}$
\noindentDesign 4. Bimodal Gaussian mixture: $\theta_i\sim 0.5N(-5/\sqrt{29},4/29)+0.5N(5/\sqrt{29},4/29)$
\noindentDesign 5. Spike-and-slab: $\theta_i\sim 0.9 N(0,0.05^2)+0.1N(0,3.159^2)$
\noindentDesign 6. 3-point distribution that places masses $(0.25,0.5,0.25)$ on $(-\sqrt{2},0,\sqrt{2})$
\noindentDesign 7. The least favorable distribution for AKP robust EBCI
\noindentDesign 8. The least favorable distribution for Cox-Morris parametric EBCI
All designs for the mixing distribution $G$ are normalized to have unit variance, so differences across designs are driven by the shape of $G$. Design 1 is the Gaussian benchmark under which the Cox-Morris normal-normal model is correctly specified. Designs 2 and 3 introduce heavier tails and hence represent progressively larger departures from the Gaussian benchmark. Designs 1 and 3 have supersmooth mixing densities in the Fourier sense, so they are relatively favorable cases for nonparametric deconvolution (\citeNP{lacour2006rates}), and thus for feasible NP-EBCI. Design 4 is included to examine how prior multimodality affects the posterior quantiles that form NP-EBCI. Design 5 is a nearly sparse prior with most mass concentrated near zero and a small fraction of large effects.
Designs 6-8 are discrete priors, thus violate the absolute continuity condition on $G$ ((ref)), which underlies our theoretical analysis of NP-EBCI. Design 7 is the least favorable prior for the AKP robust EBCI (Design 5, \citeNP{armstrong2022robust}). Among priors with the same second moment, it attains the smallest marginal coverage and therefore comes closest to the worst-case coverage bound. Design 8 is the least favorable prior for Cox-Morris parametric EBCI (Design 6, \citeNP{armstrong2022robust}), under which the interval attains its lowest marginal coverage.
For each EBCI method, we focus on the average marginal coverage and the average length reduction relative to the naive $z$-interval, both averaged across units $i=1,\dots,n$ and across all Monte Carlo replications. The nominal coverage level is $1-\alpha=95\%$ throughout.
(ref) summarizes the full simulation grid; each point corresponds to one $\text{DGP}\times \text{sample size}\times \text{SNR}$ cell.\footnote{Monte Carlo standard errors are small: across the three EBCI procedures, they are at most 0.17 percentage points for average coverage and 0.45 percentage points for average length reduction.} First, all three EBCI procedures deliver substantial average length reductions relative to the naive $z$-interval, reflecting the classical empirical Bayes gain from borrowing strength across related units (\citeNP{morris1983parametric}). Second, (ref) reveals a clear ordering among the three EBCI procedures in terms of average coverage and length. AKP (second moment) is conservative: many of its points lie above the nominal $95\%$ coverage line, but its length reductions are comparatively modest. Cox-Morris parametric EBCI often delivers the largest length reductions, but many of these gains come with substantial undercoverage. NP-EBCI achieves large length reductions while keeping empirical marginal coverage much closer to the nominal level.
(ref) further shows that the average coverage and length of EBCIs depend strongly on the signal-to-noise ratio (SNR). Lower-SNR cells tend to exhibit larger average length reductions for all three EBCIs, which is consistent with the fact that empirical Bayes pooling is most valuable when noise is large relative to signal (\citeNP{morris1983parametric}). These same cells also generate the sharpest differences in coverage across methods. At SNR=$0.1$, AKP tends to lie above the nominal $95\%$ line, Cox-Morris falls below it in many cells, and NP-EBCI remains much closer to the nominal $95\%$ line. This pattern motivates a closer examination of the low-SNR regime in (ref).
(ref) focuses on the regime $\text{SNR}=0.1$. Panel A reports the results for small samples $n=100$. NP-EBCI maintains average coverage near the nominal $95\%$ level across all designs (0.934 to 0.960), while reducing the average length by $53\%$ to $58.6\%$ relative to the naive $z$-interval. Although Cox-Morris EBCIs achieve larger average length reductions ($70.2\%$ to $74.0\%$), these gains come with substantial undercoverage in most cells, even in the Gaussian design, coverage falls to 0.848. One possible explanation is that when the prior variance $A$ is small relative to $\sigma_i^2$, the linear shrinkage factor $w_i=A/(A+\sigma_i^2)$ is also correspondingly small. Consequently, if the estimate $\hat{A}$ is biased downward in small samples, both the estimated posterior means and variances shrink aggressively toward zero, producing artificially narrow intervals (see e.g. \citeA{li2010adjusted}).
Panel B of (ref) reports the results for large samples $n=1000$. NP-EBCI continues to maintain average coverage close to the nominal $95\%$ level (0.941 to 0.972), while reducing average length by $57.3\%$ to $67.9\%$. By contrast, AKP (second moment) substantially overcovers in nearly every design (0.973 to 0.998) except its own least-favorable design, and yields modest length reductions ($44.0\%$-$48.3\%$). While Cox-Morris EBCI produces the shortest intervals (reductions of $67.8\%$ to $70.8\%$) and improves coverage compared to Panel A, it still undercovers markedly in some specific least-favorable designs. Furthermore, as \citeA{armstrong2022robust} emphasize, the Cox-Morris EBCI's empirical Bayes coverage can drop to 0.74 for the nominal $95\%$ level in the low-SNR limit.
The strong finite-sample performance of feasible NP-EBCI in average marginal coverage may seem at odds with the logarithmic ill-posedness established in (ref). It is useful to notice that the bandwidth choice in (ref) is designed to minimize average interval length subject to an average coverage requirement, rather than to minimize the mean squared errors of estimating posterior quantiles. Such coverage-oriented regularization may keep average marginal coverage close to the nominal level, with the finite-sample cost of ill-posedness showing up primarily as excess length relative to oracle NP-EBCI rather than as severe undercoverage. For this reason, we compare the average length of feasible NP-EBCI with that of its oracle counterpart in (ref).
In (ref), we define “relative feasible-oracle length gap” as
where $\overline{\text{Length}}(\mathrm{CI}_i^{\mathrm{NP}})$ and $\overline{\text{Length}}(\mathrm{CI}_i^{\mathrm{NP}^{*}} )$ denote, respectively, the average lengths of the feasible NP-EBCI in (ref) and its oracle counterpart in (ref).
(ref) first shows that the relative feasible-oracle length gap of NP-EBCI changes only marginally with sample size $n$. This aligns with the logarithmic ill-posedness of posterior quantile estimation being driven primarily by the structure of the inverse problem rather than by the sample size. Notably, this feasible-oracle gap varies significantly more across prior designs than across SNRs. In every SNR panel, the gap is small for the smooth priors but largest for the discrete and spike-and-slab priors. Although the gap varies slightly with SNR within a given design, the ordering across designs is quite stable. This pattern suggests that the relevant ill-posedness is governed mainly by the shape of $G$: smooth priors are comparatively favorable, while priors with sharp local structure, such as point masses or narrow spikes, make posterior quantiles harder to estimate.
This paper proposes empirical Bayes confidence intervals based on posterior quantiles for unobservable individual effects under a fully nonparametric prior, and establishes their theoretical properties. On the one hand, the feasible NP-EBCI achieves asymptotically exact conditional and marginal coverage. On the other hand, this flexibility comes at an unavoidable statistical cost: because posterior quantiles inherit the severe ill-posedness of nonparametric deconvolution, the relevant errors in quantile estimation and in coverage probabilities vanish at logarithmic rates. Thus the accuracy of nonparametric empirical Bayes inference is fundamentally limited by the difficulty of underlying inverse problem. Despite these slow asymptotic rates, the simulations show that the feasible NP-EBCI remains much closer to nominal coverage, even in small samples, is less conservative than the AKP robust EBCI, and delivers substantial length reductions relative to the naive $z$-interval.
While this paper focuses on empirical Bayes confidence intervals, an important direction for future work is to study other empirical Bayes decision problems that involve posterior quantiles. As \citeA{walters2024empirical} emphasizes, posterior quantiles can be more appropriate when the objective is to identify units in the tails of the distribution while controlling the risk of costly selection errors. This concern arises naturally in empirical Bayes compound decision problems such as ranking and selection; see e.g. \citeA{gu2023invidious}. An interesting question for future research is whether our minimax results can be extended to characterize the statistical difficulty of such tail-oriented empirical Bayes decision problems.
A useful feature of the NP-EBCI construction is that it relies on point identification of the distribution of individual effect $G$, rather than on the normal means structure alone. This suggests another direction for future work: extending posterior-quantile-based empirical Bayes inference to a broader class of latent-variable models in which $G$ remains point identified (\citeNP{schennach2020mismeasured}). Such extensions would help clarify which parts of our theory are specific to Gaussian deconvolution and which reflect a broader difficulty of posterior-quantile-based empirical Bayes inference.