Ismaël Castillo
email: ismael.castillo@upmc.fr

Étienne Roquain
email: etienne.roquain@upmc.fr

On spike and slab empirical Bayes multiple testing

Keywords: AMS 2000 subject classifications: Primary 62C12, 62G10 Frequentist properties of Bayesian procedures, False discovery rate, sparsity, multiple testing

This paper explores a connection between empirical Bayes posterior distributions and false discovery rate (FDR) control. In the Gaussian sequence model, this work shows that empirical Bayes-calibrated spike and slab posterior distributions allow a correct FDR control under sparsity. Doing so, it offers a frequentist theoretical validation of empirical Bayes methods in the context of multiple testing. Our theoretical results are illustrated with numerical experiments.

Introduction

Context

In modern high dimensional statistical models, several aims are typically pursued, often at the same time: testing of hypotheses on the parameters of interest, estimation and uncertainty quantification, among others. Due to their flexibility, in particular in the choice of the prior, Bayesian posterior distributions are routinely used to provide solutions to a variety of such inference problems. However, although practitioners may often directly read off quantities such as the posterior mean or credible sets once they have simulated posterior draws, the question of mathematical justification of the use of such quantities, in particular from a frequentist perspective, has recently attracted a lot of attention. While the seminal papers [START_REF] Ghosal | Convergence rates of posterior distributions[END_REF], [START_REF] Shen | Rates of convergence of posterior distributions[END_REF] set the stage for the study of posterior estimation rates in general models, the case of estimation in high dimensional models has been considered only recently from the point of view of estimation, see [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF], [START_REF] Castillo | Needles and straw in a haystack: posterior concentration for possibly sparse sequences[END_REF], [START_REF] Van Der Pas | The horseshoe estimator: posterior concentration around nearly black vectors[END_REF] among others, while results on frequentist coverage of credible sets are just starting to emerge, see e.g. [START_REF] Belitser | Needles and straw in a haystack: robust empirical Bayes confidence for possibly sparse sequences[END_REF], [START_REF] Van Der Pas | Uncertainty quantification for the horseshoe (with discussion)[END_REF]. Some of the previous approaches rely on automatic data-driven calibration of the prior parameters, following the so-called empirical Bayes approach, notably [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF], estimating the proportion of significant parameters, and [START_REF] Jiang | General maximum likelihood empirical Bayes estimation of normal means[END_REF], where the full distribution function of the unknowns is estimated.

Our interest here is on the issue of multiple testing of hypotheses. Typically, the problem is to identify the active variables among a large number of candidates. This task appears in a wide variety of applied fields as genomics, neuro-imaging, astrophysics, among others. Such data typically involve more than thousands of variables with only a small part of them being significant (sparsity).

In this context, a typical aim is to control the false discovery rate (FDR), see [START_REF] Bogdan | Beyond parametrics in interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen[END_REF] below, that is, to find a selection rule that ensures that the averaged proportion of errors among the selected variables is smaller than some prescribed level α. This multiple testing type I error rate, introduced in [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF], became quickly popular with the development of high-throughput technologies because it is "scalable" with respect to the dimension: the more rejections are possible, the more false positives are allowed. A common way to achieve this goal is to compute the p-values (probability under the null that the test statistic is larger than the observed value) and to run the Benjamini-Hochberg (BH) procedure [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF], which is often considered as a benchmark procedure. In the last decades, an extensive literature aimed at studying the BH method, by showing that it (or versions of it) controls the FDR in various frameworks, see [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF][START_REF] Benjamini | Adaptive linear step-up procedures that control the false discovery rate[END_REF][START_REF] Sarkar | Stepup procedures controlling generalized FWER and generalized FDR[END_REF][START_REF] Finner | Dependency and false discovery rate: asymptotics[END_REF], among others.

In a fundamental work [START_REF] Abramovich | Adapting to unknown sparsity by controlling the false discovery rate[END_REF], Abramovich, Benjamini, Donoho and Johnstone proved that a certain hard thresholding rule deduced from the BH procedure -keeping only observations with significant p-values -satisfies remarkable risk properties: it is minimax adaptive simultaneously for a range of losses and sparsity classes over a broad range of sparsity parameters. In addition, similar results hold true for the misclassification risks, see [START_REF] Bogdan | Asymptotic bayes-optimality under sparsity of some multiple testing procedures[END_REF][START_REF] Neuvial | On false discovery rate thresholding for classification under sparsity[END_REF]. These results in particular suggest a link between FDR controlling procedures and adaptation to sparsity. Further results in this direction include [START_REF] Bogdan | SLOPE-adaptive variable selection via convex optimization[END_REF], [START_REF] Su | SLOPE is adaptive to unknown sparsity and asymptotically minimax[END_REF] for the SLOPE estimate. Here, we shall follow a questioning that can be seen as 'dual' to the former one: starting from a commonly used Bayesian procedure that is known to optimally adapt to the sparsity in terms of risk over a broad range of sparsity classes (and even, under appropriate self-similarity type conditions, to produce adaptive confidence sets), we ask whether a uniform FDR control can be guaranteed.

Setting

In this paper, we consider the Gaussian sequence model. One observes, for 1 ≤ i ≤ n,

X i = θ 0,i + ε i , (1)
for an unknown n-dimensional vector θ 0 = (θ 0,i ) 1≤i≤n ∈ R n and ε i i.i.d. N (0, 1). This model can be seen as a stylized version of an high-dimensional model. The problem is to test H 0,i : "θ 0,i = 0" against H 1,i : "θ 0,i = 0", simultaneously over i ∈ {1, . . . , n}. We also introduce the assumption that the vector θ 0 is s n -sparse, that is, is supposed to belong to the set for some sequence s n ∈ {0, 1, . . . , n}, typically much smaller than n, measuring the sparsity of the vector.

Bayesian multiple testing methodology

From the point of view of posterior distributions, one natural approach for testing is simply based on comparing posterior probabilities of the hypotheses under consideration. Yet, to do so, a choice of prior needs to be made, and for this reason it is important to carefully design a prior that is flexible enough to adapt to the unknown underlying structure (and, here, sparsity) of the model. This is one of the reasons behind the use of empirical Bayes approaches, that aim at calibrating the prior in a fully automatic, data-driven, way. Empirical Bayes methods for multiple testing have been in particular advocated by Efron (see e.g. [START_REF] Efron | Microarrays, empirical Bayes and the two-groups model[END_REF] and references therein) in a series of works over the last 10-15 years, reporting excellent behaviour of such procedures -we describe two of them in more detail in the next paragraphs -in practice. Fully Bayes methods, that bring added flexibility by putting prior on sensible hyperparameters, are another alternative. In the sequel Bayesian multiple testing procedures will be referred to as BMT for brevity.

Several popular BMT procedures rely on two quantities that can be seen as possible Bayesian counterparts of standard p-values:

• the -value: the probability that the null is true conditionally on the fact that the test statistics is equal to the observed value, see e.g. [START_REF] Efron | Empirical Bayes analysis of a microarray experiment[END_REF]; • the q-value: the probability that the null is true conditionally on the fact that the test statistics is larger than the observed value, introduced in [START_REF] Storey | The positive false discovery rate: a Bayesian interpretation and the q-value[END_REF].

(Note that the -value is usually called "local FDR". Here, we used another terminology to avoid any confusion between the procedure and the FDR.) Obviously, these quantities are well defined only if the trueness/falseness of a null hypothesis is random, which is obtained by introducing an appropriate prior distribution.

Once the prior is calibrated (in a data-driven way or not), the q-values (resp. -values) can be computed and combined to produce BMT procedures. For instance, existing strategies reject null hypotheses with:

• a -value smaller than a fixed cutoff t = 0.2 [START_REF] Efron | Size, power and false discovery rates[END_REF];

• a q-value smaller than the nominal level α [START_REF] Efron | Microarrays, empirical Bayes and the two-groups model[END_REF];

• averaged -values smaller than the nominal level α [START_REF] Müller | Optimal sample size for multiple testing: the case of gene expression microarrays[END_REF][START_REF] Sun | Oracle and adaptive compound decision rules for false discovery rate control[END_REF][START_REF] Sun | Large-scale multiple testing under dependence[END_REF].

For alternatives see, e.g., [START_REF] Sarkar | A general decision theoretic formulation of procedures controlling FDR and FNR from a Bayesian perspective[END_REF]. In particular, one popular fact is that the use of Bayesian quantities "automatically corrects for the multiplicity of the tests", see, e.g., [START_REF] Storey | The positive false discovery rate: a Bayesian interpretation and the q-value[END_REF]; while using p-values requires to use a cutoff t that decreases with the dimension n, using -values/q-values can be used with a cutoff t close to the nominal level α, without any further correction. This is well known to be valid from a decision theoretic perspective for the Bayes FDR, that is, for the FDR integrated w.r.t. the prior distribution, as we recall in Proposition 1 below.

When the hyper-parameters are estimated from the data within the BMT, the Bayes FDR is still controlled to some extent, as proved in [START_REF] Sun | Oracle and adaptive compound decision rules for false discovery rate control[END_REF][START_REF] Sun | Large-scale multiple testing under dependence[END_REF]. However, controlling the Bayes FDR does not give theoretical guarantees for the usual frequentist FDR, that is, for the FDR at the true value of the parameter, as the pointwise FDR may deviate from an integrated version thereof.

Frequentist control of BMT

In this paper, our main aim is to study whether BMT procedures have valid frequentist multiple testing properties.

A first hint has already been given in [START_REF] Storey | The positive false discovery rate: a Bayesian interpretation and the q-value[END_REF][START_REF] Efron | Microarrays, empirical Bayes and the two-groups model[END_REF]: it turns out that the BH procedure can loosely be seen as a "plug-in version" of the procedure rejecting the q-values smaller than α (namely, the theoretical c.d.f. of the p-values is estimated by its empirical counterpart). Since the BH procedure controls the (frequentist) FDR, this might suggest a possible connection between BMT and successful frequentist multiple testing procedures.

In regard to the rapidly increasing literature on frequentist validity of Bayesian procedures from the estimation perspective, the multiple testing question for BMT procedures has been less studied so far from the theoretical, frequentist, point of view. This is despite a number of very encouraging simulation performance results, see e.g. [START_REF] Müller | Optimal sample size for multiple testing: the case of gene expression microarrays[END_REF][START_REF] Cao | Bayesian optimal discovery procedure for simultaneous significance testing[END_REF][START_REF] Guindani | A Bayesian discovery procedure[END_REF][START_REF] Martin | A nonparametric Empirical Bayes framework for Large-scale significance testing[END_REF]. A recent exception is the interesting preprint [START_REF] Salomond | Risk quantification for the thresholding rule for multiple testing using gaussian scale mixtures[END_REF] that shows a frequentist FDR control for a BMT based on a continuous shrinkage prior; yet, this control holds under a certain signal-strength assumption only. One main question we ask in the present work is whether a fully uniform control (over sparse vectors) of the frequentist FDR is possible for some posterior-based BMT procedures. Also, we would like to clarify whether the final FDR control is made at, or close to, the required level α.

Spike and slab prior distributions and sparse priors

Let w ∈ (0, 1) be a fixed hyper-parameter. Let us define the prior distribution

Π = Π w,γ on R n as Π w,γ = ((1 -w)δ 0 + wG) ⊗n , (3)
where G is a distribution with a symmetric density γ on R. Such a prior is a tensor product of a mixture of a Dirac mass at 0 (spike), that reflects the sparsity assumption, and of an absolutely continuous distribution (slab), that models nonzero coefficients. This is arguably one of the most natural priors on sparse vectors and has been considered in many key contributions on Bayesian sparse estimation and model selection, see, e.g., [START_REF] Mitchell | Bayesian variable selection in linear regression[END_REF], [START_REF] George | Calibration and empirical Bayes variable selection[END_REF]. Of course, an important question is that of the choice of w and γ. A popular choice of w is data-driven and based on a marginal maximum likelihood empirical Bayes method (to be described in more details below). The idea is to make the procedure learn the intrinsic sparsity while also incorporating some automatic multiplicity correction, as discussed e.g. in [START_REF] Scott | Bayes and empirical-Bayes multiplicity adjustment in the variable-selection problem[END_REF][START_REF] Bogdan | Beyond parametrics in interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen[END_REF]. Following such imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018 an approach in a fundamental paper, Johnstone and Silverman [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF] show that, provided γ has tails at least as heavy as Laplace, the posterior median of the empirical Bayes posterior is rate adaptive for a wide range of sparsity parameters and classes, is fast to compute and enjoys excellent behaviour in simulations (the corresponding R-package EBayesThresh [START_REF] Johnstone | EbayesThresh: R Programs for Empirical Bayes Thresholding[END_REF] is widely used). Namely, if • denotes the euclidian norm and θ = θ(X) is the coordinate-wise median of the empirical Bayes posterior distribution, there exists c 1 > 0 such that sup

θ0∈ 0 [sn] E θ0 θ -θ 0 2 ≤ c 1 s n log(n/s n ). (4)
Thus, asymptotically (in the regime s n , n → ∞, s n /n → 0), it matches up to a constant the minimax risk for this problem ( [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF]). In the recent work [START_REF] Castillo | Empirical Bayes analysis of Spike and Slab posterior distributions[END_REF], the convergence of the empirical Bayes full posterior distribution (not only aspects such as median or mean) is considered, and similar results can be obtained, under stronger conditions on the tails of γ (for instance γ Cauchy works). More precisely, for Π(

• | X) = Π(• | X) the empirical Bayes posterior, one can find a constant C 1 > 0 such that sup θ0∈ 0[sn ] E θ0 θ -θ 0 2 dΠ(θ | X) ≤ C 1 s n log(n/s n ). (5)
Further, under some conditions, one can show that certain credible sets from the posterior distributions are also adaptive confidence sets in the frequentist sense [START_REF] Castillo | Spike and Slab Empirical Bayes sparse credible sets[END_REF]. Alternatively, one can also follow a hierarchical approach and put a prior on w. The paper [START_REF] Castillo | Needles and straw in a haystack: posterior concentration for possibly sparse sequences[END_REF] obtains adaptive rates for such a fully Bayes procedure over a variety of sparsity classes, and presents a polynomial time algorithm to compute certain aspects of the posterior. Empirical Bayes approaches have also been successfully applied to a variety of different sparse priors such as empirically recentered Gaussian slabs as in [START_REF] Belitser | Needles and straw in a haystack: robust empirical Bayes confidence for possibly sparse sequences[END_REF][START_REF] Belitser | Empirical bayes oracle uncertainty quantification for regression[END_REF], or the horseshoe [START_REF] Van Der Pas | Adaptive posterior contraction rates for the horseshoe[END_REF][START_REF] Van Der Pas | Uncertainty quantification for the horseshoe (with discussion)[END_REF], both studied in terms of estimation and the possibility to construct adaptive confidence sets. In [START_REF] Jiang | General maximum likelihood empirical Bayes estimation of normal means[END_REF], an empirical Bayes approach based on the 'empirical' cdf of the θs is shown to allow for optimal adaptive estimation over various sparsity classes. For an overview on the rapidly growing literature on sparse priors, we refer to the discussion paper [START_REF] Van Der Pas | Uncertainty quantification for the horseshoe (with discussion)[END_REF].

Yet, most of the previous results are concerned with estimation or confidence sets, although a few of them report empirical false discoveries, e.g. [START_REF] Van Der Pas | Uncertainty quantification for the horseshoe (with discussion)[END_REF], Figure 7, without theoretical analysis though.

Aim and results of the paper

Here we wish to find -if this is at all possible -a posterior-based procedure using a prior Π (possibly an empirical Bayes one i.e. Π = Π), that can perform simultaneous inference in that a) it behaves optimally up to constants in terms of the quadratic risk in the sense of (4) (or (5)), b) its frequentist FDR at any sparse vector is bounded from above by (a constant times) a given nominal level. More precisely, given a nominal level t ∈ (0, 1) and ϕ t a multiple testing procedure deduced from Π ( -values or q-values procedure, as listed in Section 1.3) we want to validate its use in terms of a uniform control of its false discovery rate FDR(θ 0 , ϕ t ), see [START_REF] Bogdan | Beyond parametrics in interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen[END_REF] below, over the whole parameter space. That is, we ask whether we can find C 2 > 0 independent of t such that, for n large enough, sup

θ0∈ 0[sn] FDR(θ 0 , ϕ t ) ≤ C 2 t. (6)
Our main results are as follows: for a sparsity s n = O(n υ ) with υ ∈ (0, 1),

• Theorem 1 shows that (6) holds with C 2 arbitrary small for the BMT procedure rejecting the nulls whenever the corresponding -value is smaller than t. • Theorem 2 shows that (6) holds for some C 2 > 0 for the BMT procedure rejecting the nulls whenever the corresponding q-value is smaller than t (with a slight modification if only few signals are detected).

These results hold for spike and slab priors, for γ being Laplace or Cauchy, or even for slightly more general heavy-tailed distributions. The hyperparameter ŵ is chosen according to a certain empirical Bayes approach to be specified below (with minor modifications with respect to the choice of [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF]).

In addition, from a pure multiple testing perspective, it is important to evaluate the amplitude of C 2 > 0 in [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]. Our numerical experiments support the fact that, roughly, C 2 = 1. Furthermore, Theorem 3 shows that for some subset L 0 [s n ] ⊂ 0 [s n ] (containing strong signals), we have for the q-value BMT, for any (sequence)

θ 0 ∈ L 0 [s n ], lim n FDR(θ 0 , ϕ t ) = t , (7)
so the FDR control is exactly achieved asymptotically in that case. It follows from these results (combined with previous results of [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF][START_REF] Castillo | Empirical Bayes analysis of Spike and Slab posterior distributions[END_REF]) that the posterior distribution associated to a spike and slab prior, with γ Cauchy and a suitably empirical Bayes-calibrated w, is appropriate to perform several tasks: (6) (multiple testing), ( 5)-(4) (posterior concentration in L 2 -distance). The posterior can also be used to build honest adaptive confidence sets ( [START_REF] Castillo | Spike and Slab Empirical Bayes sparse credible sets[END_REF]). The present work, focusing on the multiple testing aspect, then completes the inference picture for spike and slab empirical Bayes posteriors, confirming their excellent behaviour in simulations.

Organisation of the paper

In Section 2, we introduce Bayesian multiple testing procedures associated to spike and slab posterior distributions as well as the considered empirical Bayes choice of w. In Section 3, our main results are stated, while Section 4 contains numerical experiments and Section 5 a short discussion. Preliminaries for the proofs are given in Section 6, while the proof of the main results can be found in Section 7. The supplementary file [START_REF] Castillo | Supplement to "On spike and slab empirical Bayes multiple testing[END_REF] gathers a number of useful lemmas used in the proofs, as well as the proofs of Propositions 1-2 and Theorem 3. The sections and equations of this supplement are referred to with an additional symbol "S-" in the numbering.

Notation

In this paper, we use the following notation:

• for F a cdf, we set F = 1 -F • φ(x) = (2π) -1/2 e -x 2 /2 and Φ(x) = x -∞ φ(u)du • u n v n means that there exists constants C, C > 0 such that |v n |c ≤ |u n | ≤ C|v n | for n large enough; • u n v n means that there exists constants C > 0 such that |u n | ≤ C|v n | for n large enough; • f (y) g(y), for y ∈ A means that there exists constants C, C > 0 such that for all y ∈ A, c|g(y)| ≤ |f (y)| ≤ C|g(y)|; • f (y) g(y), as y → ∞ means that there exists constants C, C > 0 such that c|g(y)| ≤ |f (y)| ≤ C|g(y)| for y large enough; • u n ∼ v n means u n -v n = o(u n ).
Also, for τ ∈ R n , the symbols E τ (resp. P τ ) denotes the expectation (resp. probability) under θ 0 = τ in the model [START_REF] Abramovich | Adapting to unknown sparsity by controlling the false discovery rate[END_REF]. The support of θ 0 ∈ R n is denoted by S θ0 = {i : θ 0,i = 0} or sometimes S 0 for simplicity. The cardinality of the support S θ0 is denoted by σ 0 = |S 0 |.

Preliminaries

Procedure and FDR

A multiple testing procedure is a measurable function of the form ϕ(X) = (ϕ i (X)) 1≤i≤n ∈ {0, 1} n , where each ϕ i (X) = 0 (resp. ϕ i (X) = 1) codes for accepting H 0,i (resp. rejecting H 0,i ). For any such procedure ϕ, we let

FDR(θ 0 , ϕ) = E θ0 n i=1 1{θ 0,i = 0}ϕ i (X) 1 ∨ n i=1 ϕ i (X) . (8)
A procedure ϕ is said to control the FDR at level α if FDR(θ 0 , ϕ) ≤ α for any θ 0 in R n . Note that under θ 0 = 0, we have FDR(θ 0 , ϕ) = P θ0=0 (∃i : ϕ i (X) = 1), which means that an α-FDR controlling procedure provides in particular a (single) test of level α of the full null "θ 0,i = 0 for all i". As already mentioned, in the framework of this paper, our goal is a control of the FDR around the pre-specified target level, as in ( 6) or ( 7) (where t = α).

Prior, posterior, -values and q-values

Recall the definition of the prior distribution Π = Π w,γ from (3) and let

g(x) = γ(x -u)φ(u)du. (9)
imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018

The posterior distribution Π

[• | X] of θ is then explicitly given by θ | X ∼ n i=1 i (X) δ 0 + (1 -i (X)) G Xi ( 10
)
where G x is the distribution with density γ x (u) := φ(x -u)γ(u)/g(x) and

i (X) = (X i ; w, g); (11)
(x; w, g) = Π(θ 1 = 0 | X 1 = x) = (1 -w)φ(x) (1 -w)φ(x) + wg(x) . (12)
The quantities i (X), 1 ≤ i ≤ n, given by ( 11) are called the -values. Note that, although we do not emphasize it in the notation for short, the -values depend also on w and g. The -value measures locally, for a given observation X i , the probability that the latter comes from pure noise. This is why it is sometimes called 'local-FDR', see [START_REF] Efron | Empirical Bayes analysis of a microarray experiment[END_REF]. If one has in mind a range of values -i.e. those that exceed a given amplitude-, a different measure is given by the q-values defined by

q i (X) = q(X i ; w, g); (13)
q(x; w, g) = Π(θ 1 = 0 | |X 1 | ≥ |x|) = (1 -w)Φ(|x|) (1 -w)Φ(|x|) + w G(|x|) ; (14)
G(s) = +∞ s g(x)dx. (15)
The identity ( 14) relating the q-value to Φ, G is proved in Section S-3.

Assumptions

We follow throughout the paper assumptions similar to those of [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF]. The prior γ is assumed to be unimodal, symmetric and so that

| log γ(x) -log γ(y)| ≤ Λ|x -y|, x, y ∈ R; (16)
γ(y) -1 ∞ y γ(u)du y κ-1 , as y → ∞, κ ∈ [1, 2]; (17)
y ∈ R → y 2 γ(y) is bounded. ( 18
)
Conditions ( 16), ( 17) and ( 18) above are for instance true when γ is Cauchy (κ = 2, Λ = 1) or Laplace (κ = 1, Λ is the scaling parameter). As we show in Remark S-11, explicit expressions exist for g, see [START_REF] Bogdan | SLOPE-adaptive variable selection via convex optimization[END_REF], in the Laplace case. In the Cauchy case, the integral is not explicit, but in practice (to avoid approximating the integral) one can work with the quasi-Cauchy prior, see [START_REF] Johnstone | EbayesThresh: R Programs for Empirical Bayes Thresholding[END_REF], that satisfies the above conditions and corresponds to

γ(x) = (2π) -1/2 (1 -|x|Φ(x)/φ(x)); (19)
g(x) = (2π) -1/2 x -2 (1 -e -x 2 /2 ). ( 20
)
imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018

Bayesian Multiple Testing procedures (BMT)

We define the multiple procedures defined from the -values/q-values in the following way:

ϕ -val i (t; w, g) = 1 { i(X )≤t} , 1 ≤ i ≤ n; (21)
ϕ q-val i (t; w, g) = 1 {qi(X)≤t} , 1 ≤ i ≤ n, (22)
where t ∈ (0, 1) is some threshold, that possibly depends on X. As we will see in Section 6.2, these two procedures, denoted ϕ -val (t), ϕ q-val (t) for brevity, simply correspond to (hard) thresholding procedures that select the |X i |'s larger than some (random) threshold. The value of the threshold is driven by the posterior distribution in a very specific way: it depends on γ, t, and on the whole data vector X through the empirical Bayes choice of the hyper-parameter w, that automatically "scales" the procedure according to the sparsity of the data.

Controlling the Bayes FDR

If the aim is to control the FDR at some level α, a first result indicates that choosing t = α in ϕ -val (t) and ϕ q-val (t) may be appropriate, because the corresponding procedures control the Bayes FDR, that is, the FDR where the parameter θ has been integrated with respect to the prior distribution (see, e.g., [START_REF] Sarkar | A general decision theoretic formulation of procedures controlling FDR and FNR from a Bayesian perspective[END_REF]). More formally, for any multiple testing procedure ϕ, and hyper-parameters w and γ, define

BFDR(ϕ; w, γ) = R n FDR(θ, ϕ)dΠ w,γ (θ). ( 23
)
Then the following result holds.

Proposition 1. Let α ∈ (0, 1) and w ∈ (0, 1) and consider any density γ satisfying the assumptions of Section 2.3. Let ϕ = ϕ -val (α; w, g) as defined in (21) and ϕ q = ϕ q-val (α; w, g) as defined in [START_REF] Guindani | A Bayesian discovery procedure[END_REF]. Then we have

BFDR(ϕ ; w, γ) ≤ α P (∃i : i (X) ≤ α) (24)
≤ α P (∃i :

q i (X) ≤ α) = BFDR(ϕ q ; w, γ) ≤ α. (25)
This result can be certainly considered as well known, as (24) (resp. ( 25)) is similar in essence to Theorem 4 of [START_REF] Sun | Large-scale multiple testing under dependence[END_REF] (resp. Theorem 1 of [START_REF] Storey | The positive false discovery rate: a Bayesian interpretation and the q-value[END_REF]). It is essentially a consequence of Fubini's theorem, see Section S-2.1 for a proof. While Proposition 1 justifies the use of /q-values from the purely Bayesian perspective, it does not bring any information about FDR(θ 0 , ϕ ) and FDR(θ 0 , ϕ q ) at an arbitrary sparse vector θ 0 ∈ R n .

Marginal maximum likelihood

In order to choose the hyper-parameter w, we explore now the choice made in [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF], following the popular marginal maximum likelihood method. Let us imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018 introduce the auxiliary functions

β(x) = g φ (x) -1; β(x, w) = β(x) 1 + wβ(x) . (26)
A useful property is that β is increasing on [0, ∞) from β(0) ∈ (-1, 0) to infinity, see Section 6.1. The marginal likelihood for w is by definition the marginal density of X, given w, in the Bayesian setting. Its logarithm is equal to

L(w) = n i=1 log φ(X i ) + n i=1 log (1 + wβ(X i )) ,
which is a differentiable function on [0, 1]. The derivative S of L, the score function, can be written as

S(w) = n i=1 β(X i ) 1 + wβ(X i ) = n i=1 β(X i , w). ( 27
)
The function w ∈ [0, 1] → S(w) is (a.s.) decreasing and thus w ∈ [0, 1] → L(w) is (a.s.) strictly concave. Hence, almost surely, the maximum of the function L on a compact interval exists, is unique, and we can define the marginal maximum likelihood estimator ŵ by

ŵ = argmax w∈[ 1 n ,1] L(w) (a.s.). (28)
This choice of ŵ is close to the one in [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF]. The only difference is in the lower bound, here 1/n, of the maximisation interval, which differs from the choice in [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF] by a slowly varying term. This difference is important for multiple testing in case of weak or zero signal (in contrast to the estimation task, for which this different choice does not modify the results). Another slightly different choice of interval, still close to [1/n, 1], will also be of interest below.

In addition, if ŵ ∈ (1/n, 1), it solves the equation S(w) = 0 in w. However, note that the maximiser ŵ can be at the boundary and thus may not be a zero of S.

Main results

Let us first describe the -value algorithm.

Algorithm EBayesL

Input: X 1 , . . . , X n , slab prior γ, target confidence t Output: BMT procedure ϕ -val 1. Find the maximiser ŵ given by (28); 2. Compute ˆ i (X) = (X i ; ŵ, g) given by (12); 3. Return, for 1

≤ i ≤ n, ϕ -val i = 1{ ˆ i (X) ≤ t}. ( 29
)
imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018

Theorem 1. Consider the parameter space 0 [s n ] given by (2) with sparsity s n ≤ n υ for some υ ∈ (0, 1). Let γ be a unimodal symmetric slab density that satisfies (16)-( 18) with κ as in [START_REF] Efron | Microarrays, empirical Bayes and the two-groups model[END_REF]. Then the algorithm EBayesL produces as output the BMT ϕ -val defined in (29) that satisfies the following: there exists a constant C = C(γ, υ) such that for any t ≤ 3/4, there exists an integer N 0 = N 0 (γ, υ, t) such that, for any n ≥ N 0 ,

sup θ0∈ 0 [sn] FDR(θ 0 , ϕ -val ) ≤ C log log n (log n) κ/2 . ( 30
)
Theorem 1 is proved in Section 7. The proof relies mainly on two different arguments: first, a careful analysis of the concentration of ŵ, which requires to distinguish between two regimes (signal weak/moderate or strong, basically). Second, study the FDR behavior of the -value procedure taken at some sparsity parameter w (not random but depending on n) in each of these two regimes. This requires to analyse the mathematical behavior of a number of functions of w, θ 0 , uniformly over a wide range of possible sparsities, which is one main technical difficulty of our results. In particular, the concentration of ŵ is obtained uniformly over all sparse vectors with polynomial sparsity, without any strong-signal or self-similarity-type assumption, as would typically be the case for obtaining adaptive confidence sets. Such assumptions would of course simplify the analysis significantly, but the point here is precisely that a uniform FDR control is possible for rate-adaptive procedures without any assumption on the true sparse signal. The uniform concentration of ŵ is expressed implicitly and requires sharp estimates, contrary to rate results for which a concentration in a range of values is typically sufficient. In particular, some of our lemmas in the supplementary file [START_REF] Castillo | Supplement to "On spike and slab empirical Bayes multiple testing[END_REF] are refined versions of lemmas in [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF].

As a corollary, (30) entails

lim n sup θ0∈ 0[sn]
FDR(θ 0 , ϕ -val ) = 0, and this for any chosen threshold t ∈ (0, 1) in ϕ -val . From a pure α-FDR controlling point of view, while making a vanishing small proportion of errors is obviously desirable, it implies that ϕ -val is, as far as the FDR is concerned, somewhat conservative, in the sense that it does not spend all the allowed type I errors (0 instead of α) and thus will make too few (true) discoveries at the end. It turns out that in the present setting -values are not quite on the "exact" scale for FDR control. An alternative is to consider the q-value scale, as we now describe.

Algorithm EBayesq Input: X 1 , . . . , X n , slab prior γ, target confidence t Output: BMT procedure ϕ q-val 1. Find the maximiser ŵ given by (28). 2. Compute qi (X) = q(X i ; ŵ, g)

3. Return, for 1 ≤ i ≤ n, ϕ q-val i = 1{q i (X) ≤ t} (31)
We also consider the following variant of the procedure EBayesq, which is mostly the same, except that it does not allow for too small estimated weight ŵ. Set, for L n tending slowly to infinity,

ω n = L n nG( √ 2.1 log n) . ( 32
)
For instance, for γ Cauchy or quasi-Cauchy, we have

ω n (L n /n) √ log n while for γ Laplace(1) we have ω n (L n /n) exp{C √ log n}.
Algorithm EBayesq.0

Input: X 1 , . . . , X n , slab prior γ, target confidence t, sequence L n Output: BMT procedure ϕ q-val.0

1.-2. Same as for EBayesq, returning qi (X) 3. Return, for 1 ≤ i ≤ n, and ω n as in [START_REF] Sarkar | A general decision theoretic formulation of procedures controlling FDR and FNR from a Bayesian perspective[END_REF],

ϕ q-val.0 i = 1{q i (X) ≤ t}1{ ŵ > ω n } (33)
Theorem 2. Consider the same setting as Theorem 1. Then the algorithm EBayesq produces the BMT procedure ϕ q-val in (31) that satisfies the following: there exists a constant C = C(γ, υ) such that for any t ≤ 3/4, there exists an integer N 0 = N 0 (γ, υ, t) such that, for any n ≥ N 0 ,

sup θ0∈ 0[sn ] FDR(θ 0 , ϕ q-val ) ≤ Ct log(1/t). ( 34
)
In addition, the algorithm EBayesq.0 produces the BMT procedure ϕ q-val.0 in (33) that satisfies, for ω n as in [START_REF] Sarkar | A general decision theoretic formulation of procedures controlling FDR and FNR from a Bayesian perspective[END_REF] with L n → ∞, L n ≤ log n, t ≤ 3/4 and C, N 0 as before (but with possibly different numerical values), for any n ≥ N 0 ,

sup θ0∈ 0[sn ]
FDR(θ 0 , ϕ q-val.0 ) ≤ Ct.

The proof of Theorem 2 is technically close to that of Theorem 1 and is given in Section 7, see also Section 7.2 for an informal heuristic that serves as imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018 guidelines for the proof. The statements of Theorem 2 are however of different nature, because the q-value threshold t appears explicitly in the bounds (34)- [START_REF] Storey | The positive false discovery rate: a Bayesian interpretation and the q-value[END_REF], that do not vanish as n → ∞.

The two bounds [START_REF] Shen | Rates of convergence of posterior distributions[END_REF] and ( 35) differ from a log(1/t) term, which may become significant for small t. This term appears in the case where the signal is weak (only few rejected nulls), for which the calibration ŵ is slightly too large. This may not be the case using a different type of sparsity-adaptation, or a different estimate ŵ. Indeed, this phenomenon disappears when using EBayesq.0 , since ŵ is then set to 0 when it is not large enough, in which case the FDR control is shown to be guaranteed, and we retrieve a dependence in terms of a constant times the target level t.

A consequence of Theorem 2 is that an α-FDR control can be achieved with EBayesq/EBayesq.0 procedures by taking t = t(α) sufficiently small (although not tending to zero). Again, it is important to known how small the constant C > 0 can be taken in [START_REF] Shen | Rates of convergence of posterior distributions[END_REF] and [START_REF] Storey | The positive false discovery rate: a Bayesian interpretation and the q-value[END_REF]. When the signal is strong enough, the following result shows that C = 1 and the log(1/t) factor can be removed in [START_REF] Shen | Rates of convergence of posterior distributions[END_REF].

Let us first introduce a set L 0 [s n ] of 'large' signals, for arbitrary a > 1, [START_REF] Su | SLOPE is adaptive to unknown sparsity and asymptotically minimax[END_REF] with an arbitrary a > 1, for s n → ∞ and s n ≤ n υ for some υ ∈ (0, 1). Assume that γ is a unimodal symmetric slab density that satisfies (16)-( 18) with κ as in [START_REF] Efron | Microarrays, empirical Bayes and the two-groups model[END_REF]. Then, for any pre-specified level t ∈ (0, 1), EBayesq produces the BMT procedure ϕ q-val in (31) such that

L 0 [s n ] = θ ∈ 0 [s n ] : |θ i | ≥ a 2 log(n/s n ) for i ∈ S θ , |S θ | = s n . ( 36
) Theorem 3. Consider L 0 [s n ] = L 0 [s n ; a] defined by
lim n sup θ0∈L0[sn]
FDR(θ 0 , ϕ q-val ) = lim

n inf θ0∈L0[sn]
FDR(θ 0 , ϕ q-val ) = t .

In addition, EBayesq.0 with L n → ∞, satisfies the same property whenever s n /n ≥ 2ω n , for ω n as in [START_REF] Sarkar | A general decision theoretic formulation of procedures controlling FDR and FNR from a Bayesian perspective[END_REF], which is in particular the case if s n grows faster than a given power of n and L n ≤ log n.

Theorem 3, although focused on a specific regime, shows that empirical Bayes procedures are able to produce an asymptotically exact FDR control. Again, this may look surprising at first, as the prior slab density γ is not particularly linked to the true value of the parameter θ 0 ∈ L 0 [s n ] in [START_REF] Sun | Oracle and adaptive compound decision rules for false discovery rate control[END_REF]. This puts forward a strong adaptive property of the spike and slab prior for multiple testing. Remark 4. Our three results can be extended to the case where g is not of the form (9) (that is, not necessarily of the form of a convolution with the standard gaussian), but satisfies some weaker properties, see Section 6.1. This extended setting corresponds to a "quasi-Bayesian" approach where the -values (resp. q-values) are directly given by the formulas (11) (resp. (13)), without specifying a slab prior γ.

Numerical experiments

In this section, our theoretical findings are illustrated via numerical experiments. A motivation here is also to evaluate how the parameters s n , θ 0 ∈ 0 [s n ], and the hyper-parameter γ (or g) affect the FDR control, in particular the value of the constant in the bound of Theorem 2.

For this we consider n = 10 4 , s n ∈ {10, 10 2 , 10 3 } and the following two possible scenarios for θ 0 ∈ 0 [s n ]:

• constant alternatives: θ 0,i = µ if 1 ≤ i ≤ s n and 0 otherwise; • randomized alternatives: θ 0,i i.i.d. uniformly distributed on (0, 2µ) if 1 ≤ i ≤ s n and 0 otherwise.

The parameter range for µ is taken equal to {0.01, 0.5, 1, 2, . . . , 10}. The marginal likelihood estimator ŵ given by ( 28) is computed by using a modification of the function wfromx of the package EbayesThresh [START_REF] Johnstone | EbayesThresh: R Programs for Empirical Bayes Thresholding[END_REF], that accommodates the lower bound 1/n in our definition (instead of

w n = ζ -1 ( √ 2 log n), see (49)
, in the original version). The parameter γ is either given by the quasi-Cauchy prior ( 19)- [START_REF] George | Calibration and empirical Bayes variable selection[END_REF] or by the Laplace prior of scaling parameter a = 1/2 (see Remark S-11 for more details). For any of the above parameter combinations, the FDR of the procedures EBayesL, EBayesq (defined in Section 3) is evaluated empirically via 500 replications.

Figure 1 displays the FDR of the procedures EBayesL ( -values) and EBayesq (q-values). Concerning EBayesL, in all situations, the FDR is small while not exactly equal to the value 0, which seems to indicate that the bound found in Theorem 1 is not too conservative. Moreover, the quasi-Cauchy version seems more conservative than the Laplace version, which corroborates our theoretical findings (in our bound [START_REF] Salomond | Risk quantification for the thresholding rule for multiple testing using gaussian scale mixtures[END_REF], we have the factor (log n) -1 for quasi-Cauchy and (log n) -1/2 for Laplace). As for EBayesq, when the signal is large, the FDR curves are markedly close to the threshold value t when s n /n is small, which is in line with Theorem 3. However, for a weak sparsity s n /n = 0.1, the FDR values are slightly inflated (above the threshold t), which seems to indicate that the asymptotical regime is not yet reached for this value. Looking now at the whole range of signal strengths, one notices the presence of a 'bump' in the regime of intermediate values of µ, especially for the Laplace prior. However, this bump seems to disappear when s n /n decreases. We do not known presently whether this bump is vanishing with n or if this corresponds to a necessary additional constant C = C(γ, υ) > 1 (or log(1/t)) in the achieved FDR level, but we suspect that this is related to the fact that the intermediate regime was the most challenging part of our proofs. Overall, the Cauchy slab prior seems to have a particularly suitable behavior. This was not totally surprising for us as it already showed more stability than the Laplace prior in the context of estimation with the full empirical Bayes posterior distribution, as seen in [START_REF] Castillo | Empirical Bayes analysis of Spike and Slab posterior distributions[END_REF].

Finally, we provide additional experiments in the supplement, see Section S-7. The findings can be summarized as follows:

• the curves behave qualitatively similarly for randomized alternatives (second scenario);

• at least in the considered framework, the classical BMT procedure based on averaged -values [START_REF] Sun | Oracle and adaptive compound decision rules for false discovery rate control[END_REF] seems to fail to control the frequentist FDR for large signals, which is markedly different from EBayesq; • the procedure EBayesq.0 (with L n = log log n) has a global behavior similar to EBayesq, with more conservativeness for weak signal (as expected). • it is possible to uniformly improve EBayesq.0 by considering the following modification (named EBayesq.hybrid below): if w ≤ ω n , instead of rejecting no null, EBayesq.hybrid performs a standard Bonferroni correction, that is, rejects the H 0,i 's such that p i (X) ≤ t/n. Note that a careful inspection of the proof of Theorem 2 (EBayesq.0 part) shows that the bound ( 35) is still valid for EBayesq.hybrid.

Discussion

Our results show that spike and slab priors produce posterior distributions with particularly suitable multiple testing properties. One main challenge in deriving the results was to build bounds that are uniform over sparse vectors. We demonstrate that such a uniform control is possible up to a constant term away from the target control level. This constant is very close to 1 in simulations, and can even be shown to be 1 asymptotically for some subclass of sparse vectors. The results of the paper are meant as a theoretical validation of the common practical use of posterior-based quantities for (frequentist) FDR control. While the main purpose here was validation, and as such the proposed procedures were not particularly meant to improve upon the classical FDR-controlling procedures (which specifically target the FDR control, while here the starting point is the posterior, which is rather a global inference object over the whole vector θ), it is remarkable that a uniform control of the FDR very close to the target level can be obtained for the spike and slab BMT procedure in the present unstructured sparse high-dimensional model.

While many studies focused on controlling the Bayes FDR with Bayesian multiple testing procedures, this work paves the way for a frequentist FDR analysis of Bayesian multiple testing procedures in different settings. In our study, the perhaps most surprising fact is how well marginal maximum likelihood estimation combines with FDR control under sparsity: as shown in our proof (and summarized in our heuristic) the score function is linked to a peculiar equation that makes perfectly the link between the numerator and the denominator in the FDR of the q-value-based multiple testing procedure. This phenomenon has not been noticed before to the best of our knowledge. We suspect that this link is only part of a more general picture, in which the concentration of the score process in general sparse high dimensional models plays a central role. While this exceeds the scope of this paper, generalizing our results to such settings is a very interesting direction for future work. FDR of EBayesL and EBayesq procedures with threshold t ∈ {0.05, 0.1, 0.2}. α = 0.2; n = 10, 000; 500 replications; alternative all equal to µ (on the X-axis).

imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018 6. Preliminaries for the proofs

Working with general g

As noted in Remark 4, the results of Theorems 1, 2 and 3 are also true under slightly more general assumptions, that do not impose that g is coming from a γ by a convolution product. Namely, let us assume g is a positive, symmetric, differentiable density that decreases on a vicinity of +∞ [START_REF] Sun | Large-scale multiple testing under dependence[END_REF] (g decreasing on a vicinity of +∞ means that x → g(x) is decreasing for x > M , for a suitably large constant M = M (g)). Assume moreover that

|(log g) (y)| ≤ Λ, for all y ∈ R, Λ > 0; (39)
G(y) g(y) y κ-1 , as y → ∞, for some κ ∈ [1, 2]; (40)
y ∈ R → (1 + y 2 )g(y) is bounded; (41) g/φ is increasing on [0, ∞) from (g/φ)(0) < 1 to ∞; (42)
By Lemma S-9, it is worth to note that (42) implies

G/Φ is increasing on [0, ∞) from 1 to ∞. (43)
In the case where g is of the form of a convolution with γ, see (

, ( 40) and ( 41) are easy consequences of the fact g(y) γ(y) when y → ∞ and condition (42) follows from the fact that for all fixed u > 0, the function x ∈ [0, ∞) → (φ(x + u) + φ(x -u))/φ(x) is increasing, see Lemma 1 of [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF] for a detailed derivation.

A consequence of ( 39) is that g and G have at least Laplace tails g(y) ≥ g(0)e -Λy , y ≥ 0;

G(y) ≥ g(0)Λ -1 e -Λy , y ≥ 0. (45)

BMT as thresholding-based procedures

Recall the definitions ( 21) and [START_REF] Guindani | A Bayesian discovery procedure[END_REF]. Let, for any w and t in [0, 1),

r(w, t) = wt (1 -w)(1 -t) . (46)
The following quantity plays the role of threshold for -values,

ξ = (φ/g) -1 : (0, (φ/g)(0)] → [0, ∞), (47)
i.e. ξ is the decreasing continuous invert of φ/g (that exists thanks to (42)).

Simple algebra shows that for w, t ∈ [0, 1) with r(w, t) ≤ φ(0)/g(0),

i (X) ≤ t ⇔ |X i | ≥ ξ(r(w, t)). (48)
imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018

7. Proof of of the main results

Notation

The following moments are useful when studying the score function S. Let us set

m(w) = -E 0 β(X, w) = ∞ -∞ β(t, w)φ(t)dt (56)
and further denote

m 1 (τ, w) = E τ [β(X, w)] = ∞ -∞ β(t, w)φ(t -τ )dt. ( 57
)
m 2 (τ, w) = E τ [β(X, w) 2 ] = ∞ -∞ (β(t, w)) 2 φ(t -τ )dt. ( 58
)
These expectations are well defined and studied in detail in Appendix S-5, refining previous results established in [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF].

In order to study the FDR of a procedure ϕ, we introduce the notation

V (ϕ) = i: θ0,i=0 ϕ i , S(ϕ) = i: θ0,i =0 ϕ i , (59)
counting for ϕ the number of false and true discoveries, respectively.

Heuristic

Why should the marginal empirical Bayes choice of w lead to a correct control of the FDR? Here is an informal argument that will give a direction for our proofs. We consider the case of ϕ q-val here as it is expected to reject more nulls than ϕ -val and thus to have a larger FDR. First, let us note that, when there is enough signal, one can expect ŵ to be approximately equal to the solution w of the score equation in expectation E θ0 (S(w )) = 0, that is, by using [START_REF] Mitchell | Bayesian variable selection in linear regression[END_REF]

, i:θ0,i =0 m 1 (θ 0,i , w ) = (n -s n ) m(w ),
where m and m 1 are defined by ( 56) and (57), respectively, if there θ 0 has exactly s n nonzero coordinates. As seen in Section S-5, up to log-terms,

i:θ0,i =0 m 1 (θ 0,i , w ) ≈ i:θ0,i =0 Φ(ζ(w ) -θ 0,i ) + Φ(ζ(w ) + θ 0,i ) w ; m(w ) ≈ 2G(ζ(w )).
imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018

Now consider the FDR and assume that all quantities are well concentrated (in particular, take the expectation both in the numerator and denominator in ( 8)). Then, by using (55), we have FDR(θ 0 , ϕ q-val (α; ŵ, g)) ≈ FDR(θ 0 , ϕ q-val (α; w , g))

≈ i:θ0,i=0 P θ0,i (q i (X) ≤ α) i:θ0,i=0 P θ0,i (q i (X) ≤ α) + i:θ0,i =0 P θ0,i (q i (X) ≤ α) ≈ (n -s n )r(w , α) 2G (ζ(w )) (n -s n )r(w , α) 2G (ζ(w )) + i:θ0,i =0 P θ0,i (q i (X) ≤ α) ,
where we denoted q i (X) = q(X i ; w , g) and we used that χ(r(w , t)) is close to ζ(w ), as seen in Section S-4. Now, by using the definition of q i (X),

i:θ0,i =0 P θ0,i (q i (X) ≤ α) = i:θ0,i =0 Φ (χ(r(w , α)) -θ 0,i ) + Φ (χ(r(w , α)) + θ 0,i ) ≈ i:θ0,i =0 Φ (ζ(w ) -θ 0,i ) + Φ (ζ(w ) + θ 0,i ) ,
where we used again χ(r(w , t)) ≈ ζ(w ). Now using the above properties of w , the latter is

≈ w i:θ0,i =0 m 1 (θ 0,i , w ) = (n -s n )w m(w ) ≈ (n -s n )w 2G(ζ(w )).
Putting the previous estimates together yields FDR(θ 0 , ϕ q-val (α; ŵ, g))

≈ (n -s n )r(w , α) 2G (ζ(w )) (n -s n )r(w , α) 2G (ζ(w )) + (n -s n )w 2G(ζ(w )) = r(w , α) r(w , α) + w = w 1-w α 1-α w 1-w α 1-α + w ≈ α 1-α α 1-α + 1 = α.
We will see that this heuristic holds, up to some constant terms that may come in factor of the target level α.

Proof of Theorems 1 and 2

We prove results for -and q-values together. The proof for EBayesq.0 is given at the end of this section. First, let w 0 be the solution of the equation,

nw 0 m(w 0 ) = M, (60)
for M to be chosen below in the range [1, log n] (more precisely, equal to either C log(1/t) or Ct -1 log log n for a constant C independent of t and large enough; both bounds belong to the previous interval for n large enough). For any M ∈ [1, log n], this equation has always a unique solution, as m is continuous increasing (see Lemma S-22) so the map w → w m(w) increases from 0 at w = 0 to a constant at w = 1, and in particular has a continuous inverse. This implies that w 0 goes to 0 with n, which we use freely in the sequel. Also, we note that w 0 is larger than 1/n for C in the choice of M large enough. Indeed, w 0 ≥ m(1) -1 M/n by monotonicity of m. But m(1) is at most a constant, so, provided M is large enough, w 0 ≥ 1/n. Thus w 0 is always inside the interval [n -1 , 1] over which the maximiser ŵ is defined. Let ν ∈ (0, 1) and θ 0 ∈ 0 [s n ]. Recall that S 0 denotes the support of θ 0 and that σ 0 = |S 0 | denotes the exact number of nonzero coefficients of θ 0 , so that 0 ≤ σ 0 ≤ s n . The next equation, depending on the configuration θ 0 , and on the just defined w 0 , plays a key role in the proof:

i∈S0 m 1 (θ 0,i , w) = (1 -ν)(n -σ 0 ) m(w), w ∈ [w 0 , 1). ( 61
)
This equation may or may not have a solution, depending on the true θ 0 and the values of n and ν. We will now assume n ≥ N 0 for some universal constant N 0 to be determined below.

Case 1: (61) has no solution

For a given value of n, let us consider the case where (61) has no solution in w ∈ [w 0 , 1). First, the maps w ∈ [0, 1] → m(w) and w ∈ [0, 1] → m 1 (µ, w) (µ ∈ R) are continuous, see Lemmas S-22 and S-24 and, for any µ ∈ R,

|m 1 (µ, 1)| ≤ β(x) 1 + β(x) φ(x -µ)dx ≤ max x∈R β(x) 1 + β(x)
, so thatD i∈S0 m 1 (θ 0,i , 1) ≤ Cσ 0 < (1 -ν)(n -σ 0 ) m(1) for n ≥ N 0 , where we use σ 0 ≤ s n ≤ dn υ and m(1) > 0 and N 0 = N 0 (g, υ). This means

i∈S0 m 1 (θ 0,i , w) < (1 -ν)(n -σ 0 ) m(w), for w ∈ [w 0 , 1), (62)
as otherwise by the intermediate value theorem the graphs of the functions on the two sides of the previous inequality would have to cross on [w 0 , 1) and (61) would have a solution. Lemma S-3 shows that, under (62), we have

P θ0 ( ŵ > w 0 ) ≤ e -C0ν 2 M , (63)
for some constant C 0 = C 0 (g, υ). Now consider ϕ being either ϕ -val or ϕ q-val and upper-bound the FDR by the so-called family-wise error rate by distinguishing the two cases ŵ ≤ w 0 and ŵ > w 0 as follows:

FDR(θ 0 , ϕ(t; ŵ, g)) ≤ P θ0 (∃i : θ 0,i = 0, ϕ i (t; ŵ, g) = 1)

≤ P θ0 (∃i : θ 0,i = 0, ϕ i (t; w 0 , g) = 1) + P θ0 ( ŵ > w 0 ) ≤ (n -σ 0 )P θ0,i=0 (ϕ i (t; w 0 , g) = 1) + e -C0ν where we use that w → ϕ i (t; w, g) is nondecreasing, see Lemma S-7, together with a union bound.

-value part Let ξ 0 = ξ(r(w 0 , t)) and ζ 0 = ζ(w 0 ), hen (53) leads to (provided r(w 0 , t) ≤ (φ/g)(0), which holds for e.g. t ≤ 3/4 and w 0 ≤ 1/4)

FDR(θ 0 , ϕ -val (t; ŵ, g)) ≤ 2 nw 0 1 -w 0 t 1 -t g (ξ 0 ) ξ 0 + e -C0ν 2 M .
Combining the definition of w 0 and Lemma S-24, taking n large enough so that w 0 is appropriately small, with t ≤ 3/4,

FDR(θ 0 , ϕ -val (t; ŵ, g)) ≤ 5M ξ 0 g (ξ 0 ) G(ζ 0 ) t + e -C0ν 2 M . Noting that |ξ 0 -ζ 0 | 1, g(ξ 0 ) ≤ Dg(ζ 0 ) and G(ζ 0 ) ζ κ-1 0 g(ζ 0
) by Lemma S-17 and S-24, one obtains

FDR(θ 0 , ϕ -val (t; ŵ, g)) ≤ C(g)M ζ κ 0 t + e -C0ν 2 M . (65)
q-value part For the q-value case, we come back to (64) and use (55) instead of (53) to get, setting χ 0 = χ(r(w 0 , t)),

FDR(θ 0 , ϕ q-val (t; ŵ, g)) ≤ 2 nw 0 1 -w 0 t 1 -t G (χ 0 ) + e -C0ν 2 M .
As a result, by (60) and Lemma S-24, one gets for n large enough, t ≤ 3/4, FDR(θ 0 , ϕ q-val (t; ŵ, g))

≤ 5M t G(χ 0 ) G (ζ 0 ) + e -C0ν 2 M .
Now, by the last assertion of Lemma S-17, the ratio in the last display is bounded by 2 (say) provided n is large enough, which gives FDR(θ 0 , ϕ q-val (t; ŵ, g)) ≤ 10M t + e -C0ν 2 M . (66)

Case 2: (61) has a solution

In this case we denote the solution by w 1 ∈ [w 0 , 1), so that one can write

i∈S0 m 1 (θ 0,i , w 1 ) = (1 -ν)(n -σ 0 ) m(w 1 ). ( 67
)
Now consider the slightly different equation in

w i∈S0 m 1 (θ 0,i , w) = (1 + ν)(n -σ 0 ) m(w), w ∈ [0, 1). ( 68
)
imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018

Equation (68) always has a (unique) solution w 2 ∈ [0, w 1 ). To see this, first note that the case θ 0 = 0 is excluded from (67), as m 1 (0, w) = -m(w) < 0 if w = 0. By Lemma S-22, w → m 1 (µ, w) and w → m(w) are continuous and respectively decreasing and increasing (both strictly), and m(0) = 0, while it can be seen that m 1 (µ, 0) > 0 if µ = 0, see Lemma S-22. On the other hand, the value at w = 1 of the left hand side of (68) is at most σ 0 C/w σ 0 , and so is of smaller order than (1 + ν)(n -σ 0 ) m(1) n.

The purpose of w 1 , w 2 is to provide (implicit) deterministic upper and lower bounds for the random ŵ: this is the content of Lemma S-4. Additionally, the key Lemma S-5 shows that, in case where the solution w 1 of (67) exists, we have w 1 w 2 ; that is, the bounds are of the same order. q-value part Recall the notation (59). We focus on the case of q-values first. We come back to the case of -values at the end, its proof being similar. For simplicity, we write V q (w) = V (ϕ q-val (t; w, g)) and S q (w) = S(ϕ q-val (t; w, g)).

By definition of the FDR, FDR(θ 0 , ϕ q-val (t; ŵ, g))

= E θ0 V q ( ŵ) (V q ( ŵ) + S q ( ŵ)) ∨ 1 ≤ E θ0 V q ( ŵ) (V q ( ŵ) + S q ( ŵ)) ∨ 1 1{w 2 ≤ ŵ ≤ w 1 } + P θ0 [ ŵ / ∈ [w 2 , w 1 ]] .
The last expectation in the previous display is now bounded by, using first the monotonicity of the maps w → V q (w), w → S q (w), x → x/(1 + x) and x → 1/(1 + x), then bounding the indicator variable by 1, and finally combining with Lemma S-37 applied to the independent variables U = V q (w 1 ) and T = S q (w 2 ),

E θ0 V q ( ŵ) (V q ( ŵ) + S q ( ŵ)) ∨ 1 1{w 2 ≤ ŵ ≤ w 1 } ≤ E θ0 V q (w 1 ) (V q (w 1 ) + S q (w 2 )) ∨ 1 ; ≤ exp{-E θ0 S q (w 2 )} + 12 E θ0 V q (w 1 ) E θ0 S q (w 2 ) .
Next, by using the definition of V q , one writes

E θ0 V q (w 1 ) = i: θ0,i=0 2Φ(χ(r(w 1 , t))) = 2(n -σ 0 )Φ(χ(r(w 1 , t))).
Using the definition of χ, we have Φ(χ(u)) = G(χ(u))u for u ∈ (0, 1), so

Φ(χ(r(w 1 , t))) = r(w 1 , t)G(χ(r(w 1 , t))).
Then (S-21) in Lemma S-17 implies, for small enough w 1 ,

G(χ(r(w 1 , t))) ≤ 2G(ζ(w 1 )).
Combining (S-3) in Lemma S-5, that is w 1 /C ≤ w 2 ≤ w 1 , for a constant C = C(ν, υ, g) > 0 and Lemma S-19, we have (with, say,

c 1 = 1/2), (1/2)G(ζ(w 1 )) ≤ G(ζ(w 1 /C)) ≤ G(ζ(w 2 )).
imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018

Next using Lemma S-24, one obtains G(χ(r(w 1 , t))) ≤ 3 m(w 2 ), so that

E θ0 V q (w 1 ) ≤ 3(n -σ 0 ) w 1 1 -w 1 m(w 2 ) t 1 -t ≤ 3C(n -σ 0 ) w 2 1 -Cw 2 m(w 2 ) t 1 -t ≤ C * (n -σ 0 )w 2 m(w 2 )t,
because t ≤ 3/4 for some constant C * = C * (ν, υ, g) > 0. On the other hand, by definition of S q , one can write

E θ0 S q (w 2 ) = i:θ0,i =0 Φ (χ(r(w 2 , t)) -θ 0,i ) + Φ (χ(r(w 2 , t)) + θ 0,i ) . Let us introduce the set of indices, for K 1 = 2/(1 -υ), C 0 (w, K 1 ) = 1 ≤ i ≤ n : |θ 0,i | ≥ ζ(w) K 1 . (69)
Moreover, χ(r(w 2 , t)) ≤ ζ(w 2 ) by Lemma S-16. Hence,

E θ0 S q (w 2 ) ≥ i∈C0(w2,K1) Φ (ζ(w 2 ) -θ 0,i ) + Φ (ζ(w 2 ) + θ 0,i ) ≥ i∈C0(w2,K1) Φ (ζ(w 2 ) -|θ 0,i |) . (70)
Φ (ζ(w 2 ) -|θ 0,i |) ≥ Cw 2 i∈C0(w2,K1) m 1 (θ 0,i , w 2 ) ≥ Cw 2 i∈S0 m 1 (θ 0,i , w 2 ) -Cn 1-D m(w 2 ) = Cw 2 (1 + ν)(n -σ 0 ) m(w 2 ) -C n 1-D m(w 2 ) ,
where the last equality comes from (68). As a consequence, for n large enough, for a positive constant C * = C * (υ, g) > 0, we have

E θ0 S q (w 2 ) ≥ C * (n -σ 0 )w 2 m(w 2 ).
Combining the previous bounds leads to

E θ0 V q ( ŵ) V q ( ŵ) + S q ( ŵ) ∨ 1 1{w 2 ≤ ŵ ≤ w 1 } ≤ e -C * (n-σ0)w2 m(w2) + 12 C * C * t.
As w → w m(w) is increasing, and w 1 /C ≤ w 2 by Lemma S-5, we have w 2 m(w 2 ) ≥ (w 1 /C) m(w 1 /C). Recall that w 1 ≥ w 0 by definition, so Lemma S-24 together with (S-24) of Lemma S-19 imply m(w 1 /C) ≥ (1/2) m(w 1 ) ≥ (1/2) m(w 0 ).

Combining the obtained inequalities leads to

(n -σ 0 )w 2 m(w 2 ) ≥ C (n -σ 0 )w 0 m(w 0 ) ≥ C M, (71)
where the last inequality follows from the definition of w 0 . Now turning to a bound on the FDR, Lemma S-4 and the above inequality imply, with ν = 1,

P θ0 [ ŵ / ∈ [w 1 , w 2 ]] ≤ 2e -C1ν 2 nw2 m(w2) ≤ 2e -CM , (72)
for some C = C(υ, g) > 0. Conclude that in the considered case, for some constants c 1 = c 1 (υ, g), c 2 = c 2 (υ, g) > 0, FDR(θ 0 , ϕ q-val (t; ŵ, g)) ≤ c 2 t + 3e -c1M .

-value part In the case of -values, one can follow a similar argument. We write V (w) = V (ϕ -val (t; w, g)) and S (w) = S(ϕ -val (t; w, g)). Again, the maps w → V (w) and w → S (w) are monotone. As above for q-values, one deduces FDR(θ 0 , ϕ -val (t; ŵ, g)) ≤ exp{-E θ0 S (w 2 )} + 12 E θ0 V (w 1 ) E θ0 S (w 2 ) .

By definition of V and ξ, one can write Φ (ξ(r(w 2 , t)) -θ 0,i ) + Φ (ξ(r(w 2 , t)) + θ 0,i ) . Appart from the Ct term in factor, it is the same bound as for q-values, see (70). Hence, using the bound obtained above, for n large enough and c * = c * (υ, g) > 0, E θ0 S (w 2 ) ≥ c * t(n -σ 0 )w 2 m(w 2 ).

Combining the previous bounds leads to For q-values, for ν = 1 and t ≤ 3/4, we get by combining (66) and ( 73) FDR(θ 0 , ϕ q-val (t; ŵ, g)) ≤ max 10M t + e -C0M , c 2 t + 3e -c1M

E θ0 V ( ŵ) V ( ŵ) + S ( ŵ) ∨ 1 1{w 2 ≤ ŵ ≤ w 1 } ≤ e
Taking M = (C 0 ∧ c 1 ) -1 log(1/t) gives the upper bound FDR(θ 0 , ϕ q-val (t; ŵ, g)) ≤ max{C t log(1/t) + e -log(1/t) , c 2 t + 3e -log 1/t }, which is smaller than Ct log(1/t), giving the result for q-values.

In the -values case, with ζ 1 ≤ ζ 0 and setting ν = 1, we get by combining (65) and ( 74 First notice that FDR(θ 0 , ϕ q-val.0 (t; ŵ, g)) = E θ0 n i=1 1{θ 0,i = 0}ϕ q-val.0 (t; ŵ, g) 1 ∨ n i=1 ϕ q-val.0 (t; ŵ, g) = E θ0 n i=1 1{θ 0,i = 0}ϕ q-val (t; ŵ, g) 1 ∨ n i=1 ϕ q-val (t; ŵ, g)

1{ ŵ > ω n } , (75)
by definition of algorithm EBayesq.0 . The strategy of proof is similar to the q-value case. Let us take M in the definition (60) of w 0 equal to L n from the statement of Theorem 2, see [START_REF] Sarkar | A general decision theoretic formulation of procedures controlling FDR and FNR from a Bayesian perspective[END_REF] = w 0 , leading to (76). Next, on the one hand, in Case 1, the FDR is bounded by FDR(θ 0 , ϕ q-val.0 (t; ŵ, g)) ≤ P θ0 ( ŵ > ω n ) ≤ P θ0 ( ŵ > w 0 ).

By using (63), the last display is at most e -C0ν 2 Ln . On the other hand, in Case 2, we simply use that by (75), FDR(θ 0 , ϕ q-val.0 (t; ŵ, g)) ≤ FDR(θ 0 , ϕ q-val (t; ŵ, g)) ≤ c 2 t + 3e -c1Ln , which concludes the proof.

Fig 1 .

 1 Fig 1.FDR of EBayesL and EBayesq procedures with threshold t ∈ {0.05, 0.1, 0.2}. α = 0.2; n = 10, 000; 500 replications; alternative all equal to µ (on the X-axis).

First, we

  apply Corollary S-26 with K = K 1 , w = w 2 to bound each term in the sum in terms of m 1 , noting that |θ 0,i | ≥ ζ(w 2 )/K 1 by definition of the set C 0 (w 2 , K 1 ). Next, one uses Lemma S-31 restricting the suprema to w = w 2 (which is in the prescribed interval by Lemmas S-1, S-2 and S-5) and K = K 1 , to get for n large enough and constants C = C(υ, g) > 0, C = C (υ, g) > 0, D = D(υ, g) ∈ (0, 1), i∈C0(w2,K1)

E

  θ0 V (w 1 ) = 2(n -σ 0 )Φ(ξ(r(w 1 , t))).The bound Φ(u) ≤ φ(u)/u for u > 0 (see Lemma S-33), combined with the definition of ξ and that |ξ(r(w 1 , t)) -ζ(w 1 )| 1 by Lemma S-17 leads toE θ0 V (w 1 ) ≤ 3(n -σ 0 )ζ(w 1 ) -1 r(w 1 , t)g(ξ(r(w 1 , t))).Lemma S-17 then implies g(ξ(r(w 1 , t))) ≤ 2g(ζ(w 1 )) (say), for n large enough. Using w 1 /C ≤ w 2 ≤ w 1 , and (S-24) in Lemma S-19, we have(1/2)g(ζ(w 1 )) ≤ g(ζ(w 1 /C)) ≤ g(ζ(w 2 )).Next using the relation ζκ-1 g(ζ)m(w) from Lemma S-24, one obtains g(ξ(r(w1 , t))) ζ(w 2 ) 1-κ m(w 2 ) ζ(w 1 ) 1-κ m(w 2 ), so that E θ0 V (w 1 ) ≤ Ct(n -σ 0 )w 1 m(w 2 )ζ(w 1 ) -κ ≤ c * t(n -σ 0 )w 2 m(w 2 )ζ(w 1 ) -κ ,for a constant c * = c * (υ, g) > 0. On the other hand, by definition of S , one can write E θ0 S (w 2 ) = i:θ0,i =0

  imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date:August 14, 2018 Lemma S-18 now enables to bound from below the two terms in the previous display in terms of ζ(w 2 ), and further restricting the sum to the set of indices C 0 (w 2 , K 1 ) defined by (69) with the same choice of K 1 leads toE θ0 S (w 2 ) ≥ Ct i∈C0(w2,K1) Φ (ζ(w 2 ) -|θ 0,i |) .

) 7 . 3 . 3 .

 733 -c * t(n-σ0)w2 m(w2) + 12c * c * 1 ζ(w 1 ) κ .As in (71), we have (n-σ 0 )w 2 m(w 2 ) ≥ C (n -σ 0 )w 0 m(w 0 ) ≥ C M . One concludes that, in Case 2, for some constants d 1 = d 1 (υ, g), d 2 = d 2 (υ, g) > 0 and taking ν = 1, setting ζ(w 1 ) = ζ 1 , FDR(θ 0 , ϕ -val (t; ŵ, g)) ≤ d 2 ζ -κ 1 t + e -C M c * t + 2e -CM ≤ d 2 ζ -κ 1 t + 3e -d1M t .(74Combining cases 1 and 2

2 1log( 1 /w 1 ) log n and ζ 2 1

 2112 )FDR(θ 0 , ϕ -val (t; ŵ, g)) ≤ max CM ζ -κ 0 t + e -C0M , d 2 ζ -κ 1 t + 3e -d1tM ≤ d 3 {(M + 1)ζ -κ 1 t + e -d4tM }.The announced bound is obtained upon setting M = t -1 d -1 4 log(ζ κ 1 ) and noting that ζ log(1/w 1 ) log n by using Lemmas S-1, S-2 to bound w 1 and Lemma S-15 to bound ζ(w 1 ). This concludes the proof of Theorem 1 for -values and Theorem 2 for q-values. imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date: August 14, 2018 7.3.4. Proof for EBayesq.0

  , and suppose L n ∈ [1, log n]. Let us show for n large enough, ω n ≥ w 0 . (76) As ζ(w 0 ) ≤ ζ(1/n) ≤ √ 2.1 log n for n large enough by Lemmas S-1, S-15,

	ω n =	nG( √	L n 2.1 log n)	≥	L n nG(ζ(1/n))	≥	L n nG(ζ(w 0 ))	.
	Now, by using Lemma S-24, for n large enough,
	L n nG(ζ(w 0 ))	≥ 0.9	2L n n m(w 0 )	≥	L n n m(w 0 )

imsart-generic ver. 2014/10/16 file: CR2018_arXiv.tex date:August 14, 2018

Acknowledgments

This work has been supported by ANR-16-CE40-0019 (SansSouci) and ANR-17-CE40-0001 (BASICS).

* Work partly supported by grants of the ANR, projects ANR-16-CE40-0019 (SansSouci) and ANR-17-CE40-0001 (BASICS) 1 imsart-generic ver. θ i = 0} ≤ s n },

When u becomes small, the order magnitude of ξ(u) is given in Lemma S-13: ξ(u) slightly exceeds (-2 log u)

1/2 but not by much, which comes from the fact that g has heavy tails.

Another quantity close to ξ we shall use in the sequel is the threshold ζ introduced in [START_REF] Johnstone | Needles and straw in haystacks: empirical Bayes estimates of possibly sparse sequences[END_REF] and defined as, for any w ∈ (0, 1],

Combining the definitions leads, see (S-8) for details, to ζ(w) = ξ(w/(1 + w)) and ξ(w) ≤ ζ(w). Similarly, let us introduce a threshold for q-values as

which is the decreasing continuous invert of Φ/G (that exists thanks to (43)). For all w ∈ [0, 1) and t ∈ [0, 1) with r(w, t) ≤ 1,

Lemma S-14 shows that, for small u, the order of magnitude of χ(u) is slightly more than Φ -1 (u) but not by much, which comes from the fact that G has heavy tails. Also, Lemma S-10 together with (48)-(51) imply

(52)

Single type I error rates

The single type I error rates of our procedures are evaluated by the following result (proved in Section S-2.2).

Proposition 2. Consider any function g satisfying the assumptions of Section 6.1. Consider r(•, •) as in (46), ξ as in (47) and χ as in (50). Then the following bounds hold. For all t, w such that r(w, t) ≤ (φ/g)(0),

For q-values, we have, for all t, w such that r(w, t) ≤ 1,

As a result, for a fixed w, we see that the more g is heavily tailed, the more the type I error rate is large. This is well-expected, as the heavier the tails of g, the more mass the prior puts on large values.


Samir M Perlaza

Ali Tajer

H Vincent Poor

Simultaneous Information and Energy Transmission: A Finite Block-Length Analysis

Keywords: Information and Energy Transmission, Information and Power Transfer, Finite Block-Length Regime

In this paper, a non-asymptotic analysis of the fundamental limits of simultaneous information and energy transmission (SIET) is presented. The notion of the informationenergy capacity region, i.e., the largest set of simultaneously achievable information and energy rates, is revisited in a context in which transmissions occur within a finite number of channel uses and strictly positive decoding error probability (DEP) and energy shortage probability (ESP) are tolerated. The focus is on the case of one transmitter, one information receiver and one energy harvester communicating through binary symmetric memoryless channels. In this case, some outer bounds on the information transmission rate and the energy transmission rate are presented. More specifically, given a finite block-length, a DEP, and an ESP, four scenarios arise depending on whether an average or maximal probability constraint is imposed on the DEP and the ESP. For each scenario, the limits on the information rate and energy rate beyond which a transmission is no longer possible are presented (impossibility results). These results reveal the competitive interaction between the information transmission and energy transmission tasks identifying a certain regime in which increasing the information rate necessarily implies decreasing the energy rate and vice versa.

I. INTRODUCTION

Simultaneous information and energy transmission (SIET) refers to a communication system in which a set of transmitters aim to simultaneously carry on two tasks: information transmission to a set of information receivers (IRs); and energy transmission to a set of energy harvesters (EHs). The performance of SIET is often measured by the information and energy transmission rates that can be simultaneously achieved under certain reliability constraints. Reliability can be measured by two metrics: decoding error probability (DEP); and energy shortage probability (ESP). The fundamental limits of SIET consist of the largest set of information and energy rate tuples that can be simultaneously achieved. This set is Samir M. Perlaza (samir.perlaza@inria.fr) is with the CITI Laboratory, a joint laboratory between the Institut National de Recherche en Informatique et en Automatique (INRIA), the Université de Lyon and the Institut National de Sciences Apliquées (INSA) de Lyon. 6 Av. des Arts 69621 Villeurbanne, France. He is also with the Department of Electrical Engineering at Princeton University, Princeton, NJ 08544 USA.

Ali Tajer is with the Department of Electrical, Computer, and Systems Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180.

H. Vincent Poor (poor@princeton.edu) is with the Department of Electrical Engineering at Princeton University, Princeton, NJ 08544 USA.

This research was supported in part by the European Commission under Marie Skłodowska-Curie Individual Fellowship No. 659316 and Euro-Mediterranean Cooperation ERA-NET project COM-MED, and in part by the U.S. National Science Foundation under Grants CNS-1702808 and ECCS-1647198.

often referred to as the information-energy capacity region [START_REF] Belhadj Amor | Fundamental limits of simultaneous energy and information transmission[END_REF]. Traditionally, information-energy capacity regions are calculated subject to the fact that both DEP and ESP must be arbitrarily close to zero. This strong reliability constraint leads to the unavoidable use of infinitely long communication blocks and thus, these fundamental limits are meaningful only under the assumption that the communication might last a very long time. Within this context, the case of one transmitter and one co-located IR and EH is studied in [START_REF] Varshney | Transporting information and energy simultaneously[END_REF], [START_REF] Grover | Shannon meets Tesla: Wireless information and power transfer[END_REF], and [START_REF] Varasteh | Wireless information and power transfer over an AWGN channel: Nonlinearity and asymmetric Gaussian signaling[END_REF]. The multiuser case is studied in [START_REF] Belhadj Amor | Feedback enhances simultaneous wireless information and energy transmission in multiple access channels[END_REF], [START_REF] Park | Joint wireless information and energy transfer in a two-user MIMO interference channel[END_REF], [START_REF] Fouladgar | On the transfer of information and energy in multi-user systems[END_REF], [START_REF] Ulukus | Energy harvesting wireless communications: A review of recent advances[END_REF], [START_REF] Popovski | Interactive joint transfer of energy and information[END_REF], [START_REF] Zhang | MIMO broadcasting for simultaneous wireless information and power transfer[END_REF], [START_REF] Gastpar | Gaussian multiple-access channels under received-power constraints[END_REF], [START_REF] Krikidis | Simultaneous wireless information and power transfer in modern communication systems[END_REF], [START_REF] Ding | Power allocation strategies in energy harvesting wireless cooperative networks[END_REF], [START_REF] Khalfet | Simultaneous information and energy transmission in Gaussian interference channels[END_REF], [START_REF]Simultaneous information and energy transmission in Gaussian interference channels with feedback[END_REF], and references therein. Despite the existing literature, these approaches quickly lose relevance in scenarios in which communications must occur within a short period, e.g., the Internet of things.

More relevant fundamental limits, from an engineering perspective, are those that take into account that SIET occurs within a finite number of channel uses and the system tolerates strictly positive DEP and ESP. This paper introduces this novel approach building upon the existing results on the fundamental limits of information transmission in the non-asymptotic block-length regime (c.f. [START_REF] Polyanskiy | Channel coding: Non-asymptotic fundamental limits[END_REF] and [START_REF] Polyanskiy | Channel coding rate in the finite blocklength regime[END_REF]). The main result is the characterization of the information-energy capacity region of an SIET system with one transmitter, one IR, and one EH communicating over binary symmetric memoryless channels. More specifically, for a given number of channel uses and the biggest DEP and ESP that can be tolerated, the largest set of simultaneously achievable information and energy rates is characterized. The trade-off between information rate and energy rate, first reported in [START_REF] Varshney | Transporting information and energy simultaneously[END_REF] and [START_REF] Grover | Shannon meets Tesla: Wireless information and power transfer[END_REF] in the asymptotic regime, is studied in the non-asymptotic regime. The system parameters for which information and energy transmission are conflicting tasks are also characterized. These results reveal the central role of the ESP in both the information transmission rate and the energy transmission rate.

The paper is organized as follows. Sec. II formulates the problem and introduces the notion of the information-energy capacity region in the non-asymptotic block-length regime. Sec. III focuses on the case of binary symmetric memoryless channels and presents the main results. Sec. IV concludes this work.

II. SYSTEM MODEL

Consider a three-party communication system in which a transmitter aims at simultaneously sending information to an IR and energy to an EH through a noisy communication medium. Such a system can be modeled by a random transformation (X , Y × Z, P Y Z|X ) from an input alphabet X to an output alphabet Y × Z consisting of a transition probability kernel P Y Z|X . That is, given an input x ∈ X , the output (y, z) ∈ Y ×Z is observed with probability P Y Z|X (y, z|x). In the following, the noisy communication medium is represented by a memoryless channel. A memoryless channel is a random transformation

(X n , Y n × Z n , P Y Z|X ), (1)
where n ∈ N is the block length and

Y = (Y 1 , Y 2 , . . . , Y n ) ∈ Y n , Z = (Z 1 , Z 2 , . . . , Z n ) ∈ Z n
, and X = (X 1 , X 2 , . . . , X n ) ∈ X n are n-dimensional vectors of random variables, such that given an input x = (x 1 , x 2 , . . . , x n ), the output (y 1 , y 2 , . . . , y n , z 1 , z 2 , . . . , z n ) is observed with probability

P Y Z|X (y, z|x) = n t=1 P Y Z|X (y t , z t |x t ). (2)
Within this context, two tasks are carried out by the transmitter: (a) the information transmission task; and (b) the energy transmission task.

A. Information Transmission Task

The message to be sent from the transmitter to the IR is a realization of a random variable that is uniformly distributed in {1, 2, . . . , M }, with 1 M < ∞. To carry out this task within n channel uses, the transmitter uses an (n, M )-code.

Definition 1 ((n, M )-code): An (n, M )-code for the random transformation in (1) is a system

{(u(1), D 1 ), (u(2), D 2 ), . . . , (u(M ), D M )} , (3)
where for all (i, j) ∈ {1, 2, . . . , M } 2 , with i = j,

u(i) (u 1 (i), u 2 (i), . . . , u n (i)) ∈ X n , ( 4a
) D i ∩ D j = ∅; and (4b) D 1 ∪ D 2 ∪ . . . ∪ D M ⊆ Y n . ( 4c
)
Given the system in (3), to transmit the message index i, the transmitter inputs the symbol u t (i) to the channel at time t ∈ {1, 2, . . . , n}. The IR observes at the end of channel use t, the output y t . At the end of n channel uses, the IR states that the symbol i was transmitted if it satisfies the rule

(y 1 , y 2 , . . . , y n ) ∈ D i . (5)
The probability of error associated with the transmission of message index i, denoted by

λ i ∈ [0, 1], is λ i Pr [Y ∈ D c i | X = u(i)] , (6)
where the probability is with respect to the marginal P Y |X of the joint probability distribution in (2). The average probability of error, denoted by λ, is

λ 1 M M m=1 λ m . ( 7
)
Given a parameter ∈ [0, 1], information transmission is said to be reliably performed in two different senses:

(a) the average decoding error probability is smaller than , i.e., λ < , or (b) the maximum decoding error probability is smaller than , i.e., for all i ∈ {1, 2, . . . , M }, λ i < . An (n, M )-code that exclusively satisfies the first criterion is said to be an (n, M, )-code with average decoding error probability. Alternatively, an (n, M )-code that satisfies both criteria is said to be an (n, M, )-code with maximal decoding error probability.

B. Energy Transmission Task

Let g : Z → R + be a positive-valued function. The amount of energy delivered to the EH by the channel outputs z = (z 1 , z 2 , . . . , z n ), given by the function B n :

Z n → R + , is B n (z) = n t=1 g(z t ). ( 8
)
The objective of the transmitter is to ensure that a minimum amount of energy b is harvested at the EH at the end of n channel uses. An energy shortage is said to occur when the energy harvested at the EH is less than the minimum required at the end of the transmission. The probability of an energy shortage when transmitting the message index i ∈ {1, 2, . . . , M }, denoted by θ i , is

θ i Pr [B n (Z) < b | X = u(i)] , (9)
where the probability is with respect to the marginal P Z|X of the joint probability distribution in [START_REF] Varshney | Transporting information and energy simultaneously[END_REF]. The average probability of energy shortage, denoted by θ, is

θ 1 M M i=1 θ i . ( 10
)
Given a parameter δ ∈ [0, 1], energy transmission is said to be reliably performed in two different senses: (c) the average energy shortage probability is smaller than δ, i.e., θ < δ, or (d) the maximum energy shortage probability is smaller than δ, i.e., for all i ∈ {1, 2, . . . , M }, θ i < δ. An (n, M, )-code that exclusively satisfies the first criterion is said to be an (n, M, , δ, b)-code with average energy shortage probability. Alternatively, an (n, M, )-code that satisfies both criteria is said to be an (n, M, , δ, b)-code with maximal energy shortage probability.

C. Fundamental Limits

The non-asymptotic fundamental limits of SIET are described by the notion of the information-energy capacity region.

↵2 ↵2 0 0 1 1 0 0 1 1 ↵1 1 ↵1 Source Transmitter W Ut(W ) Yt Z t c W Receiver EH ↵1 1 ↵1 1 ↵2 ↵2
Fig. 1. A three-party communication system in which a transmitter aims at simultaneously sending information to an information receiver and energy to an energy harvester through binary symmetric channels.

The information-energy capacity region C in Definition 2 is a subset of R 5 that is difficult to characterize. In this work, given fixed parameters (n, , δ), the focus is on the subset C(n, , δ) that contains the pairs

(M, b) ∈ N × R + such that (n, M, , δ, b) ∈ C. Therefore, if (M, b) ∈ C(n, , δ),
then the information rate log 2 (M ) n bits per channel use, and the energy rate b n energy units per channel use, are jointly achievable within n channel uses with a (maximal and/or average) decoding error probability and (maximal and/or average) energy shortage probability δ.

The subset C(n, , δ) is a subset of R 2 . In order to emphasize the reliability criterion, when needed, the informationenergy capacity region C(n, , δ) is often written as

C (a,c) (n, , δ), C (a,d) (n, , δ), C (b,c) (n, , δ) or C (b,d) (n, , δ), respectively.
A similar notation is used for the larger region C.

III. MEMORYLESS BINARY SYMMETRIC CHANNELS

This section focuses on the case in which X = Y = Z = {0, 1} and the random transformation in (1) is such that for all x ∈ {0, 1} n , the outputs y ∈ {0, 1} n and z ∈ {0, 1} n are observed with probability

P Y Z|X (y, z|x) = n t=1 P Y |X (y t |x t )P Z|X (z t |x t ), (11)
where for all (x, y, z) ∈ {0, 1} 3 ,

P Y |X (y|x)=α 1 1 {x =y} + (1 -α 1 )1 {x=y} , (12)
P Z|X (z|x)=α 2 1 {x =z} + (1 -α 2 )1 {x=z} , (13)
and

α 1 ∈ [0, 1
2 ) and α 2 ∈ (0, 1 2 ) are the crossover probabilities of the channel. This scenario is depicted by Figure 1. Let

b 0 = g(0) 0 and b 1 = g(1) 0, (14)
be the energy harvested when the output of the channel at the EH is 0 and 1, respectively. The case in which b 0 = b 1 is trivial, since the harvested energy is always nb 1 = nb 0 energy units, independently of the codebook. This implies that the information transmission task can be carried out without taking into account the energy transmission task. Hence, for avoiding the trivial cases, the following assumption is adopted without loss of generality:

b 1 < b 0 . (15)
For all z ∈ {0, 1} n , it follows that

B n (z)=b 0 N (0|z) + b 1 N (1|z) (16) =(b 0 -b 1 )N (0|z) + nb 1 , (17)
where N (0|z) and N (1|z) are the numbers of zeros and ones in the vector z, respectively. Note that B n (z) is bounded for all z ∈ Z n , i.e.,

nb 1 B n (z) nb 0 . (18)
The inequalities in [START_REF] Feller | An Introduction to Probability Theory and Its Application[END_REF] imply that there exists a case in which energy transmission might occur with zero (maximal or average) energy shortage probability at a given energy rate b n b 1 . This is because the event B n (Z) < nb 1 is observed with zero probability. Basically, transmitting any symbol either zero or one is indifferent from the energy transmission perspective. In this case, the information transmission task can be carried out independently of the energy transmission task given that nb 1 energy units can always be reliably transmitted in n channel uses. Alternatively, any energy transmission rate b n > b 0 cannot be achieved with an average or maximal energy shortage probability strictly smaller than one.

Given an (n, M )-code described by the system in (3), let the empirical probability distribution of the channel input symbols induced by the codeword u(i), denoted by P (i) X , be such that

P (i) X (0) = 1 - P (i) X (1) = 1 n n t=1 1 {ut(i)=0} , (19)
for all i ∈ {1, 2, . . . , M }. Let also the empirical distribution of the channel input symbols jointly induced by all codewords, denoted by PX (0), be such that

PX (0) = 1 -PX (1) = 1 nM M i=1 n t=1 1 {ut(i)=0} . ( 20
)
Using these empirical distributions some upper bounds can be obtained on both the energy and information transmission rates.

A. An Outer Bound on the Energy Rate of (n, M )-Codes Define Q : R → [0, 1] as the complementary cumulative distribution function of the standard normal distribution, and define Q -1 : (0, 1) → R as the functional inverse of Q. Using this notation, the following proposition provides an outer bound on the number of energy units b that can be reliably delivered by any given (n, M )-code with an average or maximal energy shortage probability δ ∈ (0, 1).

Proposition 1: Consider an (n, M, , δ, b)-code described by the system in (3) for the random transformation in [START_REF] Gastpar | Gaussian multiple-access channels under received-power constraints[END_REF] satisfying [START_REF]Simultaneous information and energy transmission in Gaussian interference channels with feedback[END_REF]. Then, subject to a maximal energy shortage probability constraint, it holds that for all i ∈ {1, 2, . . . , M }, b<n

Å (b 0 -b 1 ) (1 -2α 2 ) P (i) X (0) + α 2 + b 1 ã - » n(b 0 -b 1 ) 2 α 2 (1 -α 2 )Q -1 (δ) + O(1), (21)
and subject to an average energy shortage probability constraint, it holds that b<n

Å (b 0 -b 1 ) (1 -2α 2 ) PX (0) + α 2 + b 1 ã - » n(b 0 -b 1 ) 2 α 2 (1 -α 2 )Q -1 (δ) + O(1). ( 22
)
Proof of Proposition 1: Note that an (n, M, , δ, b)-code with maximal energy shortage probability satisfies for all i ∈ {1, 2, . . . , M },

Pr [B n (Z) < b | X = u(i)] < δ, ( 23
)
where the probability operator in (23) applies with respect to the marginal P Z|X of the joint distribution in [START_REF] Gastpar | Gaussian multiple-access channels under received-power constraints[END_REF]. Note also that the random variable B n (Z) in ( 23) is the sum of n independent binary random variables 1 {Zt=0} , where Z t follows the distribution P Z|X=ut(i) . That is,

Pr [B n (Z)<b | X = u(i)] = Pr n t=1 1 {Zt=0} < Å b -nb 1 b 0 -b 1 ã .
Hence, from the Berry-Essen theorem [START_REF] Feller | An Introduction to Probability Theory and Its Application[END_REF], it follows that for all i ∈ {1, 2, . . . , M },

Pr [B n (Z) < b | X = u(i)] - (1 -α 2 ) 2 + α 2 2 2 nα 2 (1 -α 2 ) +Q Ñ n Ä (1 -2α 2 ) P (i) X (0) + α 2 ä -b-nb1 b0-b1 nα 2 (1 -α 2 ) é . (24)
Thus, from ( 23) and ( 24), it follows that b<n

Å (b 0 -b 1 ) (1 -2α 2 ) PX (0) + α 2 + b 1 ã - » n(b 0 -b 1 ) 2 α 2 (1 -α 2 )Q -1 δ + (1 -α 2 ) 2 + α 2 2 2 nα 2 (1 -α 2 ) <n Å (b 0 -b 1 ) (1 -2α 2 ) PX (0) + α 2 + b 1 ã - » n(b 0 -b 1 ) 2 α 2 (1 -α 2 )Q -1 δ + O(1), (25)
where the last inequality in (25) follows from a Taylor approximation of the function Q -1 . A similar procedure can be performed to prove (22). This completes the proof. For fixed parameters (n,

α 2 , δ, b 0 , b 1 ), let ρ * : R + → R + be defined as ρ * (b) = min 1, ρ + (b) , (26)
with

ρ + (b) = b -n (1 -α 2 )b 1 + α 2 b 0 n(b 0 -b 1 )(1 -2α 2 ) + α 2 (1 -α 2 ) √ n(1 -2α 2 ) Q -1 (δ)+O Å 1 n ã .
Using this notation, the following corollary follows immediately from Proposition 1.

Corollary 1: Consider an (n, M, , δ, b)-code described by the system in (3) for the random transformation in [START_REF] Gastpar | Gaussian multiple-access channels under received-power constraints[END_REF] satisfying [START_REF]Simultaneous information and energy transmission in Gaussian interference channels with feedback[END_REF]. Then, subject to a maximal energy shortage probability constraint, for all i ∈ {1, 2, . . . , M }, the empirical input distribution

P (i) X satisfies P (i) X (0)>ρ * (b), (27)
and subject to an average energy shortage probability constraint, it follows that

PX (0)>ρ * (b), (28)
where ρ * (b) is defined in (26). Corollary 1 leads to interesting conclusions by noticing that ρ * (b) is a lower bound on the fraction of zeros in each codeword (maximal energy shortage probability) or the fraction of zeros among all codewords (average energy shortage probability) when energy is transmitted at an average energy rate b n . This is natural from the perspective of the assumption in [START_REF]Simultaneous information and energy transmission in Gaussian interference channels with feedback[END_REF], which implies that the symbol zero carries more energy that the symbol one.

Note that the input distribution that achieves the largest information transmission rate is the uniform distribution [START_REF] Polyanskiy | Channel coding: Non-asymptotic fundamental limits[END_REF]. That is, PX (0) = 1 -PX (1) = 1 2 . Hence, Proposition 1 provides an outer bound on the energy rate that can be transmitted by an (n, M, , δ, b)-code that possesses an empirical input distribution that is information-rate optimal. The following corollary describes this observation.

Corollary 2: Consider an (n, M, , δ, b)-code for the random transformation in [START_REF] Gastpar | Gaussian multiple-access channels under received-power constraints[END_REF] satisfying [START_REF]Simultaneous information and energy transmission in Gaussian interference channels with feedback[END_REF]. Assume that such a code exhibits an information-rate optimal empirical distribution. Then, it follows that b < b(n, δ), with b :

N × [0, 1] → R, such that b(n, δ)=n Å b 0 + b 1 2 ã - » n (b 0 -b 1 ) 2 α 2 (1 -α 2 )Q -1 (δ)
+O [START_REF] Belhadj Amor | Fundamental limits of simultaneous energy and information transmission[END_REF].

Essentially, Corollary 2 determines a threshold on the number of energy units b beyond which an (n, M, , δ, b)-code, if it exists, exhibits a conflict between the energy transmission task and the information transmission task. More specifically, if there exists an (n, M, , δ, b)-code whose energy transmission rate b is beyond the threshold b, it exhibits an empirical input distribution for which PX (0) > PX [START_REF] Belhadj Amor | Fundamental limits of simultaneous energy and information transmission[END_REF]. This implies than a zero is transmitted more often than a one, which is not information-rate optimal.

Proposition 1 also provides an outer bound on the largest energy rate that can be transmitted by any (n, M, )-code. Note that the largest energy-transmission rate is achieved by a zero information-rate code whose codewords contain only zeros, i.e., PX (0) = 1 -PX (1) = 1. The following corollary describes this observation.

Corollary 3: Consider an (n, M, , δ, b)-code for the random transformation in [START_REF] Gastpar | Gaussian multiple-access channels under received-power constraints[END_REF] satisfying [START_REF]Simultaneous information and energy transmission in Gaussian interference channels with feedback[END_REF]. Then, it follows that:

b < b(n, δ), with b : N × [0, 1] → R, such that b(n, δ)=n (1 -α 2 )b 0 + α 2 b 1 - » n (b 0 -b 1 ) 2 α 2 (1 -α 2 )Q -1 (δ) + O(1).(30)
Finally, from Corollary 2 and Corollary 3, it might be expected that the energy transmission task and the information transmission task exhibit a conflicting interaction. The following section explores this particular interaction.

B. An Outer Bound on the Information Rate of (n, M )-Codes Given an (n, M, , δ, b)-code, the following proposition describes a bound on M . This bound does not take into account the decoding error probability and thus, it might be loose in some cases. However, it plays an important role when b(n, δ) < b < b(n, δ).

Proposition 2: Consider an (n, M, , δ, b)-code for the random transformation in [START_REF] Gastpar | Gaussian multiple-access channels under received-power constraints[END_REF] satisfying [START_REF]Simultaneous information and energy transmission in Gaussian interference channels with feedback[END_REF]. Then, subject to a maximal energy shortage probability constraint, it holds that

M Ç n nρ * (b) å 2 (n-nρ * (b) ) , (31)
and subject to an average energy shortage probability constraint, it follows that

M Ç M n M nρ * (b) å 2 (M n-M nρ * (b) ) , (32)
where ρ * (b) is defined by (26).

Proof of Proposition 2: Corollary 1 provides an approximation to the minimum number of zeros in each codeword in any given (n, M, , δ, b)-code with maximal probability of energy shortage. That is, for all

i ∈ {1, 2, . . . , M } if follows that N (0|u(i)) nρ * (b) . (33)
This immediately provides the upper-bound on M given that all codewords must contain at least nρ * (b) zeros. Hence, the right-hand side of (31) is the maximum number of codewords of length n for which at least nρ * (b) symbols are zeros. The same argument can be used to prove the fixed-point equation in (32). This completes the proof. Note that ρ * (b) is monotonically increasing with the energy rate b. Interestingly, when ρ * (b) ∈ ( 1 2 , 1], the right-hand sides of (31) and (32) are monotonically decreasing with b. This highlights the existing trade-off between the information transmission task and the energy transmission task. That is, in the regime in which ρ * (b) ∈ ( 1 2 , 1], increasing the energy rate would necessarily imply decreasing the information rate.

IV. CONCLUSIONS

In this paper, the fundamental limits of SIET have been studied under the assumption that the transmission occurs during a finite number of channel uses at the expense of strictly positive DEP and ESP. From this perspective, a non-asymptotic fundamental limit has been introduced: the information-energy capacity region, that is, the largest set of jointly achievable energy and information rates. The focus has been on the case of one transmitter, one IR and one EH communicating via binary symmetric memoryless channels. In this case, given a finite block-length, a DEP, and an ESP, four scenarios have been observed depending on whether an average or maximal probability constraint is imposed on the DEP and the ESP. For each scenario, the limits on the information rate and energy rate beyond which a transmission is no longer possible have been characterized. These results have revealed the competitive interaction between the information transmission task and energy transmission task. In particular, a certain regime in which increasing the information rate necessarily implies decreasing the energy rate and vice versa has been identified.

Definition 2 (

 2 Information-Energy Capacity Region): The information-energy capacity region C ∈ N 2 × [0, 1] 2 × R of the random transformation in (1) is the set of all tuples (n, M, , δ, b) for which there exists an (n, M, , δ, b)-code satisfying at least one of the reliability criterion pairs: (a, c); (a, d); (b, c); or (b, d).

1

 1