Gilles Blanchard

Thorsten Dickhaus

Etienne Roquain

Fanny Villers

ON LEAST FAVORABLE CONFIGURATIONS FOR STEP-UP-DOWN TESTS

Keywords: False discovery rate, least favorable configuration, multiple testing, Steck's recursions, step-up-down

This paper investigates an open issue related to false discovery rate (FDR) control of step-up-down (SUD) multiple testing procedures. It has been established

in earlier literature that for this type of procedure, under some broad conditions, and in an asymptotical sense, the FDR is maximum when the signal strength under the alternative is maximum. In other words, so-called "Dirac uniform configurations" are asymptotically least favorable in this setting. It is known that this property also holds in a non-asymptotical sense (for any finite number of hypotheses), for the two extreme versions of SUD procedures, namely step-up and step-down (with extra conditions for the step-down case). It is therefore very natural to conjecture that this non-asymptotical least favorable configuration property could more generally be true for all "intermediate" forms of SUD procedures. We prove that this is, somewhat surprisingly, not the case. The argument is based on the exact calculations proposed earlier by Roquain and Villers (2011a), that we extend here by generalizing Steck's recursion to the case of two populations. Secondly, we quantify the magnitude of this phenomenon by providing a nonasymptotic upper-bound and explicit vanishing rates as a function of the total number of hypotheses.

Introduction

In mathematical statistics, so-called least favorable parameter configurations (LFCs) play a pivotal role. For a given statistical decision problem over a parameter space Θ and a given decision rule δ, we define an LFC θ * (δ) as any element of Θ that maximizes the risk (expected loss) of δ under this parameter, i.e., ∀θ ∈ Θ, R(θ, δ) ≤ R(θ * (δ), δ), where R(θ, δ) denotes the risk of rule δ under θ. When available, the knowledge of an LFCs allows one to obtain a bound on the risk over a possibly very large parameter space, including non-or semi-parametric cases where Θ has infinite dimensionality. Theoretical investigations of minimax properties can rely on the computation of an LFC. Such knowledge is also relevant for practice, because a user of the procedure can be provided with a performance guarantee if an LFC is known. In this case, even if the risk under the LFC cannot be computed in closed form, it can be approximated by a Monte-Carlo method simulating the distribution corresponding to the LFC. Finally, if the parameter space Θ is partitioned into disjoint, restricted submodels, it can be of interest to gain knowledge of the LFC of a decision rule δ separately over each submodel, thus providing finer-grained information.

LFC considerations naturally occur in hypothesis testing problems. A classical example is that of one-sided tests over a one-dimensional parameter space admitting an isotonic likelihood ratio: it is well-known that the LFC for the type I error probability is located at the boundary of the null hypothesis. This fact is used to derive critical values for uniformly more powerful tests in that setting.

The LFC problem is particularly delicate for multiple hypotheses testing and the latter has been investigated by many authors in previous literature [START_REF] Finner | On the false discovery rate and expected type I errors[END_REF][START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF][START_REF] Lehmann | Generalizations of the familywise error rate[END_REF][START_REF] Finner | Dependency and false discovery rate: asymptotics[END_REF][START_REF] Romano | Control of generalized error rates in multiple testing[END_REF][START_REF] Guo | On control of the false discovery rate under no assumption of dependency[END_REF][START_REF] Somerville | Step-up and step-down procedures controlling the number and proportion of false positives[END_REF][START_REF] Finner | On the false discovery rate and an asymptotically optimal rejection curve[END_REF][START_REF] Finner | Controlling the familywise error rate with plugin estimator for the proportion of true null hypotheses[END_REF][START_REF] Gontscharuk | Asymptotic and Exact Results on FWER and FDR in Multiple Hypothesis Testing[END_REF]. In that setting, a family of m ≥ 2 null hypotheses H 1 , . . . , H m is to be tested simultaneously under the scope of a common statistical model with parameter space Θ, and some type I error criterion is used that accounts for multiplicity. For theoretical as well as practical applications, it is relevant to determine LFCs over the restricted parameter spaces Θ m,m 0 where exactly m 0 out of m of the null hypotheses are true. In this setting, LFC results can be derived straightforwardly only in special situations.

In the present work, we restrict our attention to multiple testing procedures that depend on the observed data only through a collection of marginal p-values, each associated to an individual null hypothesis. This is a commonly used setting for multiple testing problems in high dimension. Moreover, we consider procedures that reject exactly those null hypotheses having their p-value less than a certain common threshold t * , which can possibly be data-dependent. That is, t * may depend in a complex way of the entire family of p-values. We call such procedures threshold-based for short.

In this setting, LFCs crucially depend on the type I error criterion considered.

One frequently encountered family of such criteria is given through loss functions that only depend on the number of type I errors, denoted In other words, the risk takes the form R(θ, δ) = E θ [φ(V m )]. Natural assumptions are that t * is a nonincreasing function in each p-value and that φ is a nondecreasing function. Then, by additionally assuming that the p-values are jointly independent, it is known that the LFC over Θ m,m 0 is a Dirac-uniform (DU) distribution (introduced by [START_REF] Finner | On the false discovery rate and expected type I errors[END_REF], i. e. , such that p-values corresponding to true nulls are independent uniform variables, while p-values under alternatives follow a Dirac distribution with point mass 1 in zero. This result is formally restated in Appendix C. For example, this LFC property holds true under the above assumptions for the k-family-wise error rate (k-FWER). For a given θ ∈ Θ, the k-FWER under θ is defined by FWER k,θ := P θ (V m ≥ k). Strong control of the (1-)family-wise error rate, i. e., ensuring that sup θ FWER 1,θ ≤ α for a pre-defined level α ∈ (0, 1), is the usual type I error concept in traditional multiple hypotheses testing theory.

V m = V m (θ, δ) := |{1 ≤ i ≤ m : H i is
However, over the last two decades, progress in application fields such as genomics, proteomics, neuroimaging, and astronomy has lead to massive multiple testing problems with very large systems of hypotheses [START_REF] Dudoit | Multiple testing procedures with applications to genomics[END_REF][START_REF] Pantazis | A comparison of random field theory and permutation methods for statistical analysis of meg data[END_REF][START_REF] Miller | Controlling the false-discovery rate in astrophysical data analysis[END_REF]. In this type of applications, (k-)FWER control is typically too strict a requirement, and a less stringent notion of type I error control is needed in order to ensure reasonable power of corresponding multiple tests. In particular, the false discovery rate (FDR) introduced by [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] has become a standard criterion for type I error control in large-scale multiple testing problems. The FDR is defined as the expected proportion of type I errors among all rejections.

Unfortunately, it does not fall into the class of type I error measures considered in the previous paragraph, so that the above result does not apply. Furthermore, because the average ratio of two dependent random variables is not necessarily increasing in the value of the numerator, the LFC problem for the FDR criterion turns out to be a challenging issue -even for simple classes of multiple tests and under independence assumptions.

In this work, we contribute to the theory of LFCs under the FDR criterion for so-called step-up-down multiple tests (SUD procedures, for short). These procedures constitute an important and wide subclass of threshold-based multiple testing procedures, wherein the threshold t * is obtained by comparing the reordered p-values to a fixed set of critical values (see [START_REF] Tamhane | A generalized step-up-down multiple test procedure[END_REF][START_REF] Sarkar | Some results on false discovery rate in stepwise multiple testing procedures[END_REF]. Furthermore, recent research has reinforced the interest of this type of procedures. For instance, Finner et al. (2009) have shown that step-up-down tests can be used is association with the so-called asymptotically optimal rejection curve (AORC) to provide an asymptotically (as m → ∞) valid FDR control which is additionally optimal in some specific sense.

Namely, the contributions of the paper are as follows:

• a survey of known LFC results for SUD procedures in specific model classes is provided in Section 3;

• new results and surprising counterexamples for LFCs of SUD procedures are derived in Section 4.

• in Section 5, Steck's recursion is extended to the case of two populations and we provide a summary of the exact formulas for computing the FDR proposed by Roquain and Villers (2011a,b); these formulas are used to derive the counterexamples previously mentioned. In what follows, we denote by F the set containing c.d.f.'s from [0, 1] into [0, 1] that are continuous.

Mathematical setting

Models

• The p-value family p follows the (two group) fixed mixture model with parameters m ≥ 2, 1 ≤ m 0 ≤ m and F ∈ F, for which the corresponding distribution is denoted by FM(m, m 0 , F ), if p = (p i , i ∈ {1, ..., m}) is a family of mutually independent variables and for all i,

p i ∼ U (0, 1) if 1 ≤ i ≤ m 0 , F if m 0 + 1 ≤ i ≤ m,
where U (0, 1) denotes the uniform distribution on (0, 1).

• The p-value family p follows the (two group) random mixture model with parameters m ≥ 2, π 0 ∈ [0, 1] and F ∈ F, for which the corresponding distribution is denoted by RM(m, π 0 , F ), if there is an (unobserved) binomial random variable m 0 ∼ B(m, π 0 ) such that p follows the FM(m, m 0 , F ) model conditionally on m 0 . In that case, the p-values are i.i.d. with (unconditional) c.d.f. G(t) = π 0 t + (1π 0 )F (t).

In the above definition, note that the true nulls are automatically assigned to the m 0 (random or not) first coordinates. This can be assumed without loss of generality, since we only consider procedures which only depend on p-values through their reordering in increasing order.

A common additional assumption on F is that F (x) ≥ x or that F is concave.

For instance, these assumptions are both satisfied in the two following standard examples:

-Gaussian location model:

F (t) = Φ(Φ -1 (t) -µ)
, for a given alternative mean µ > 0, where Φ(z) = P(Z ≥ z) for Z ∼ N (0, 1). This corresponds to the alternative distribution of p-values when testing for µ ≤ 0 under a Gaussian location shift model with unit variance.

-Dirac δ 0 distribution: F is identically equal to 1, as introduced by [START_REF] Finner | On the false discovery rate and expected type I errors[END_REF]. The corresponding distribution in the FM model is called Dirac-uniform (DU) configuration (or distribution) and denoted by FM(m, m 0 , F ≡ 1) or simply DU(m, m 0 ). We define similarly RM(m, π 0 , F ≡ 1). Note that the Dirac-uniform configuration can be seen as an instance of the Gaussian c.d.f. for an alternative mean µ = ∞.

In the existing literature, the Dirac-uniform distribution has often be considered as the first candidate for being an LFC of several global type I error rates (with or without a theoretical support) (see, e.g., [START_REF] Finner | Dependency and false discovery rate: asymptotics[END_REF][START_REF] Romano | Control of generalized error rates in multiple testing[END_REF][START_REF] Somerville | Step-up and step-down procedures controlling the number and proportion of false positives[END_REF].

Procedures

In this paper, we consider the particular class of multiple testing procedures called

(k) ≤ t k } if p (λ) > t λ . (2)
In the sequel, for convenience, we identify procedures with their rejection sets, e.g., SUD λ (t) = {1 ≤ i ≤ m : p i ≤ t k}. An important remark is that the cases λ = 1 and λ = m correspond to the traditional step-down and step-up procedures, respectively. An illustration is provided in Figure 1.

A classical choice for the threshold collection t consists of Simes' (1986) critical values t k = αk/m for a pre-specified level α ∈ (0, 1). The corresponding step-up-down procedure is called the linear step-up-down procedure and is denoted by LSUD λ . In particular, for λ = 1 and λ = m, the procedure LSUD λ is simply denoted by LSD and LSU, respectively. LSU corresponds to the famous linear step-up procedure of [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF].

It is common to consider threshold collections of the form

t k = ρ(k/m) for a function ρ : [0, 1] → [0, 1].
This function is generally assumed to satisfy the following assumptions:

ρ : [0, 1] → [0, 1] is continuous and non-decreasing;

(3)

x ∈ (0, 1] → ρ(x)/x is non-decreasing. ( 4
)
The function ρ is called the critical value function (while its inverse is generally called the rejection curve, see e.g. Finner et al., 2009). Observe that assumption

(3) can always be assumed to hold when m is fixed, however it is of interest in the case of an asymptotical analysis as m → ∞ (in which case ρ is assumed to be independent of m). Assumption (4) on the other hand restricts the possible threshold collection also for any fixed m. It will be often used in this paper.

For a fixed finite m, assumptions (3) and ( 4) taken together are equivalent to "k → t k /k is non-decreasing".

False discovery rate and LFCs

As introduced by [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF], the false discovery rate of a multiple testing procedure is defined as the averaged ratio of the number of erroneous rejections to the total number of rejections. In our setting, for a distribution P being either FM(m, m 0 , F ) or RM(m, π 0 , F ), the FDR of a stepup-down procedure can be written as

FDR(SUD λ (t), P ) = E p∼P [FDP(SUD λ (t), m 0 , p)] , (5)
for which the FDP is the false discovery proportion defined by

FDP(SUD λ (t), m 0 , p) = |{1 ≤ i ≤ m 0 : p i ≤ t k}| |{1 ≤ i ≤ m : p i ≤ t k}| ∨ 1 , (6)
where | • | denotes the cardinality function and in which m 0 is either fixed or ran-

dom whether P is FM(m, m 0 , F ) or RM(m, π 0 , F ), respectively. For short, the quantity FDR(SUD λ (t), FM(m, m 0 , F )) is often denoted FDR(SUD λ (t), m 0 , F ),
or FDR(SUD λ (t), F ) when the context makes the interpretation unambiguous.

Similarly, FDR(SUD λ (t), RM(m, π 0 , F )) can be shortened as FDR(SUD λ (t), π 0 , F ), or FDR(SUD λ (t), F ).

Definition 2.2. Any F ′ ∈ F is called a least favorable configuration (LFC) for the FDR of SUD λ (t) in a fixed mixture model with m 0 true hypotheses out of m and over the class

F if ∀F ∈ F, FDR(SUD λ (t), m 0 , F ) ≤ FDR(SUD λ (t), m 0 , F ′ ).
A similar definition holds for a random mixture model with m hypotheses and proportion π 0 of true hypotheses.

The above definition can possibly be restricted to a subclass F ′ ⊂ F (typically, the class of concave c.d.f.s). This will be clearly specified in the context.

Obviously, if F ′ is an LFC for the fixed mixture model for all values of m 0 , then it is also an LFC in the RM(m, π 0 , F ) model for any value of π 0 (by integrating over m 0 ∼ B(m, π 0 )).

Survey of known LFCs under the FDR criterion

Recent results about LFCs for the FDR criterion related to step-up-down type procedures are summarized in Figure 2 of [START_REF] Gontscharuk | Asymptotic and Exact Results on FWER and FDR in Multiple Hypothesis Testing[END_REF], we easily derive the following result:

Theorem 3.1. [START_REF] Gontscharuk | Asymptotic and Exact Results on FWER and FDR in Multiple Hypothesis Testing[END_REF]] Consider a step-up-down procedure using a threshold collection of the form t k = ρ(k/m), where ρ satisfies (3) and (4).

Assume that the step-up-down procedure is performed at an order

λ = λ m such that λ m /m → κ ∈ [0, 1]. Assume that m 0 /m → ζ ∈ [0,
FDR(SUD λm (t), F ) -FDR(SUD λm (t), F ≡ 1) ≤ 0, ( 7
)
either for all ζ ∈ [0, 1] if κ > 0 or for all ζ ∈ [0, 1) if κ = 0.
However, for a finite m, and λ / ∈ {1, m} no result is known about LFC's to our knowledge. This is the point of the paper and is symbolized by the question mark in the middle of Figure 2.

Finally, let us consider the linear SUD procedure, that is, the SUD procedure using the threshold collection t k = αk/m, α ∈ (0, 1) (corresponding to ρ(x) = αx). Since both LSU and LSD procedures satisfy that DU is an LFC and since an SUD procedure can be expressed as a combination of an SU and an SD, we might make the following conjecture, which is the starting point of this paper.

Conjecture 3.2. For any m ≥ 2, the Dirac-uniform configuration (F ≡ 1) is a least favorable configuration for the FDR of the linear step-up-down procedure, in the RM(m, π 0 , F ) and FM(m, m 0 , F ) models.

Obviously, a similar conjecture might be formulated for a (non-linear) stepup-down procedure using ρ satisfying (3) and (4). The exact calculations described in Section 5 allow to compute the value of FDR(LSUD λ (t)) exactly. This shows the following (numerical) result.

Proposition 4.1. For m = 10, consider the linear step-up-down procedure LSUD λ at level α = 0.5 and for λ ∈ {4, 5, 6, 7}. Then, we have

FDR(LSUD λ , F ) > FDR(LSUD λ , F ≡ 1), ( 8
)
in either of the two following cases:

• in the FM(m, m 0 , F ) model, with m 0 = 7 and F (x) = x;

• in the RM(m, π 0 , F ) model, with π 0 = 7/10 and F (x) = x.

This disproves Conjecture 3.2 and shows that finding the LFC for SUD is more difficult than for SU and SD separately. More generally, Figure 3

Nonasymptotic bound

In the present section, we investigate the amplitude of the phenomenon observed in the previous section as a function of the number of hypotheses. In other words, we derive a more explicit and non-asymptotic version of the limit appearing in (7). For this, we consider a perturbation analysis of the SUD procedure as defined In order to state the result in a compact form, we first introduce the following notation for an SUD threshold in a continuous setting.

U (ℓ, G) =    min {u ∈ [ℓ, 1] : G(ρ(u)) ≤ u} if G(ρ(ℓ)) ≥ ℓ; max {u ∈ [0, ℓ] : G(ρ(u)) ≥ u} if G(ρ(ℓ)) < ℓ. (9)
Observe that in the above definition, the infimum and supremum are welldefined since the considered sets are non-empty; using that G is non-decreasing, it can be seen that U (ℓ, G) is a fixed point of the function G • ρ (so that the infimum is indeed a minimum and the supremum, a maximum). Unfortunately, the number of rejections k of the SUD procedure as defined in (2) does not always satisfy k/m = U (λ/m, Ĝm ) (because of the step-down part, see Figure 6 in Section 7.2). Nevertheless, the following lemma is proved in Section 7.2.

Lemma 4.3. With the above notation, if the threshold collection t is defined as

t k = ρ(k/m), we have U (λ/m, Ĝm ) ≤ k/m ≤ U(λ/m, ( Ĝm + m -1 ) ∧ 1), ( 10
)
where k is defined by (2) and Ĝm (x) := m -1 m i=1 1{p i ≤ x} is the empirical c.d.f. of the p-values.

We now state our main result.

Theorem 4.4. Consider a threshold collection t of the form t k = ρ(k/m), where ρ : [0, 1] → [0, 1] satisfies (3) and (4). Let ζ ∈ (0, 1) be an arbitrary constant.

For δ ∈ (0, 1), define 4 for an illustration). Let us define the remaining term: for any y ∈ (0, 1),

u + δ = U (λ, (G DU ζ + δ) ∧ 1) and u - δ = U (λ, (G DU ζ -δ) ∨ 0), where G DU ζ (x) := (1 -ζ) + ζx (see Figure
ε(δ, m, ζ, y) := ρ(u + δ ) -ρ(u - δ ) u + δ + 4 1 -ζ e -2m(δ-y-1 m ) 2 + (1-y/ζ) + . ( 11
)
Then, for any F ∈ F and λ ∈ {1, . . . , m} the following holds.

• In the FM(m, m 0 , F ) model with 0 < m 0 < m and ν = max

k∈{m 0 -1,m 0 } {|k/m- ζ|} ∈ [0, 1], we have FDR(SUD λ (t), m 0 , F ) ≤ FDR(SUD λ (t), m 0 , F ≡ 1) + m 0 m ε(δ, m, ζ, ν); (12)
• In the RM(m, π 0 , F ) model with π 0 = ζ, we have for any γ ∈ (0, 1),

FDR(SUD λ (t), π 0 , F ) ≤ FDR(SUD λ (t), π 0 , F ≡ 1) + π 0 ε(δ, m, ζ, γ) + 4e -2m(γ-1/m) 2 + . ( 13
u - δ = 1-ζ-δ 1-αζ ∨ 0 and u + δ = 1-ζ+δ 1-αζ ∧ 1. Hence, (ρ(u + δ ) -ρ(u - δ ))/u + δ ≤ 2αδ 1 -ζ + δ .
As a result, ( 12) and ( 13) hold by replacing (ρ(u + δ )ρ(u - δ ))/u + δ by 2αδ 1-ζ+δ inside the remaining term ε(δ, m, ζ, y).

Second, for ζ > α, let us use

ρ(u) = αu/(1 -u(1 -α)), that is, ρ -1 (t) = t/(α + t(1 -α)). (14)
The rejection curve ρ -1 , displayed in the right panel of Figure 4, is called the "asymptotically optimal rejection curve" (AORC). It was introduced by Finner et al. (2009) for the purpose of improving the power of the linear critical value function. When applying Theorem 4.4 with this choice of ρ, the calculation of u - δ and u + δ depends on the position of the parameter λ m /m on [0, 1] and ρ(u + δ )-ρ(u - δ ) may not vanish when δ becomes small, see Figure 4. Fortunately, when ρ(λ m /m) (dotted-long-dashed line) is smaller than the quantity ρ(v δ ) (dashed line), the two points ρ(u - δ ) and ρ(u + δ ) (dotted-dashed lines) are expected to be close as δ becomes small and the bound given in Theorem 4.4 will vanish. The exact expressions of u - δ , u + δ and v δ can be easily derived by solving the corresponding quadratic equations. For short, we only report the equivalent as δ tends to 0:

v δ = 1 -δα/(ζ -α) + O(δ 2 ) u + δ = (1 -ζ)/(1 -α) + δζ/(ζ -α) + O(δ 2 ) u - δ = (1 -ζ)/(1 -α) -δζ/(ζ -α) + O(δ 2 ). Since ρ ′ ((1 -ζ)/(1 -α)) = α/ζ 2 , we have ρ(u + δ ) -ρ(u - δ ) = 2αδ/(ζ 2 -αζ) + O(δ 2
). As a result, assuming ζ > α and λ/m < v δ , we can derive that (12) and ( 13) hold and that quantity (ρ(u

+ δ ) -ρ(u - δ ))/u + δ is equivalent to 2α(1-α) (ζ 2 -αζ)(1-ζ)
δ (as δ tends to zero) in the remaining term ε(δ, m, ζ, y).

Convergence rate when m → ∞

We can now use Theorem 4.4 in an asymptotic sense and for specific critical value functions, in order to obtain an explicit bound on the convergence rate of the limit appearing in (7). For ζ m tending to 1 (not too quickly, "fairly" sparse case), the convergence rate is slower.

As a counterpart, assumptions of Corollary 4.5 are more restrictive than those of Theorem 3.1. In particular, they exclude the case where ζ m tends to 1 faster than (log m)/m ("highly" sparse case). This is a limitation of the methodology used to prove the nonasymptotic results. This problem may possibly be fixed by adapting the proof of Lemma 3.7 of [START_REF] Gontscharuk | Asymptotic and Exact Results on FWER and FDR in Multiple Hypothesis Testing[END_REF] to a nonasymptotic setting, but this falls outside of the intended scope of this paper.

Exact formulas

In this section, we gather some of the formulas derived by Roquain and Villers (2011a,b) to calculate the joint distribution of the number of false discoveries and the number of discoveries. Moreover, we complement this work by giving a new recursion that makes these formulas fully usable. These calculations are used to state Proposition 4.1.

A new Steck's recursion

For any k ≥ 0 and any threshold collection t = (t 1 , . . . , t k ), we denote

Ψ k (t) = Ψ k (t 1 , . . . , t k ) = P U (1) ≤ t 1 , . . . , U (k) ≤ t k . ( 16
)
where (U i ) 1≤i≤k is a sequence of k variables i.i.d. uniform on (0, 1) and with the convention Ψ 0 (•) = 1. In practice, quantity ( 16) can be evaluated using standard and Wellner, 1986, p. 366-369).

Steck's recursion Ψ k (t) = (t k ) k -k-2 j=0 k j (t k -t j+1 ) k-j Ψ j (t 1 , . . . , t j ) (Shorack
Next, we generalize the latter to the case of two populations. Define for 0 ≤ k 0 ≤ k and any threshold collection t = (t 1 , . . . , t k ),

Ψ k,k 0 ,F (t 1 , . . . , t k ) = P U (1) ≤ t 1 , . . . , U (k) ≤ t k , (17)
where (U i ) 1≤i≤k is a sequence of k variables such that (U i ) 1≤i≤k 0 are i.i.d. uniform on (0, 1), independently of (U i ) k 0 +1≤i≤k i.i.d. of c.d.f. F and with the convention Ψ 0,0,F (•) = 1. The computation of Ψ k,k 0 ,F is more difficult than Ψ k because it involves non i.i.d. variables. To our knowledge the existing formulas for computing Ψ k,k 0 ,F have a complexity exponential with k [START_REF] Glueck | Fast computation by block permanents of cumulative distribution functions of order statistics from several populations[END_REF]. Here, we propose a substantially less complex computation, by generalizing Steck's recursions as follows.

Proposition 5.1. The following recursion holds: for 0 ≤ k 0 ≤ k,

Ψ k,k 0 ,F (t 1 , . . . , t k ) =(t k ) k 0 F (t k ) k-k 0 - 0≤j 0 ≤j≤k-2 j 0 ≤k 0 j-j 0 ≤k-k 0 k 0 j 0 k -k 0 j -j 0 × (t k -t j+1 ) k 0 -j 0 (F (t k ) -F (t j+1 )) k-k 0 -j+j 0 Ψ j,j 0 ,F (t 1 , . . . , t j ), (18)
with the convention 0 0 = 1. This is proved in Section 7.1. Note that the case k = k 0 reduces to the standard (one population) Steck's recursion.

FDR formulas

Using the Ψ k 's and Ψ k,k 0 's, let us define the following useful quantities: for any

threshold collection t = (t k ) 1≤k≤m , F ∈ F; for any π 0 ∈ [0, 1], k ≥ 0, k ≤ m, 0 ≤ j ≤ k, we let P m,π 0 ,F (t, k, j) = m j m -j k -j π j 0 π k-j 1 (t k ) j (F (t k )) k-j × Ψ m-k (1 -G(t m ), . . . , 1 -G(t k+1 )); (19)
P m,π 0 ,F (t, k, j) = m j m -j k -j π j 0 π k-j 1 (1 -G(t k+1 )) m-k × Ψ k,j,F (t 1 , . . . , t k ), (20)
where

G(t) = π 0 t + (1 -π 0 )F (t). For any m 0 ∈ {0, . . . , m}, k ≥ 0, k ≤ m, j ≤ m 0 , k -j ≤ m -m 0 , we let Q m,m 0 ,F (t, k, j) = m 0 j m -m 0 k -j (t k ) j (F (t k )) k-j × Ψ m-k,m 0 -j,F (1 -t m , . . . , 1 -t k+1 ); (21) Q m,m 0 ,F (t, k, j) = m 0 j m -m 0 k -j (1 -t k+1 ) m 0 -j (1 -F (t k+1 )) m-m 0 -k+j × Ψ k,j,F (t 1 , . . . , t k ), (22)
where

F (t) = 1 -F (1 -t).
The following results have been proved by Roquain and Villers (2011a,b).

Theorem 5.2 (Roquain and Villers (2011)). Consider any threshold collection t and the quantities defined by (19), ( 20), ( 21) and (22). Then the following holds:

(i) In the RM(m, π 0 , F ) model, for any π 0 ∈ [0, 1], F ∈ F, 0 ≤ k ≤ m, 0 ≤ j ≤ k, P(|R ∩ {1, . . . , m 0 }| = j, |R| = k) = P m,π 0 ,F (t, k, j) for R = SU(t), P m,π 0 ,F (t, k, j) for R = SD(t). ( 23
) (ii) In the FM(m, m 0 , F ) model, for any m 0 ∈ {0, . . . , m}, F ∈ F, 0 ≤ k ≤ m, 0 ∨ (k -m + m 0 ) ≤ j ≤ m 0 ∧ k, P(|R ∩ {1, . . . , m 0 }| = j, |R| = k) = Q m,m 0 ,F (t, k, j) for R = SU(t), Q m,m 0 ,F (t, k, j) for R = SD(t). (24)
Classically, any step-up-down procedure can be written as a combination of a step-down and a step-up procedure [START_REF] Sarkar | Some results on false discovery rate in stepwise multiple testing procedures[END_REF]:

SUD λ (t) = SU((t λ ∧ t j ) 1≤j≤m ) if |SU((t λ ∧ t j ) 1≤j≤m )| < λ, SD((t λ ∨ t j ) 1≤j≤m ) if |SD((t λ ∨ t j ) 1≤j≤m )| ≥ λ. (25)
Moreover, since {|SU(( 25) form a partition of the probability space. This yields the following explicit FDR computations:

t λ ∧ t j ) 1≤j≤m )| < λ} = {p (λ) > t λ } and {|SD((t λ ∨ t j ) 1≤j≤m )| ≥ λ} = {p (λ) ≤ t λ } the two cases in (
Corollary 5.3. let λ ∈ {1, . . . , m} and consider any threshold collection t. Then the following

(i) In the model RM(m, π 0 , F ), for any π 0 ∈ [0, 1], F ∈ F, we have FDR(SUD λ (t)) = λ-1 k=1 k j=0 j k P m,π 0 ,F (t ∧ t λ , k, j) + m k=λ k j=0 j k P m,π 0 ,F (t ∨ t λ , k, j). (26)
(ii) In the model FM(m, m 0 , F ), for any m 0 ∈ {0, . . . , m}, F ∈ F, we have

FDR(SUD λ (t)) = λ-1 k=1 m 0 ∧k j=0∨(k-m+m 0 ) j k Q m,m 0 ,F (t ∧ t λ , k, j) + m k=λ m 0 ∧k j=0∨(k-m+m 0 ) j k Q m,m 0 ,F (t ∨ t λ , k, j). (27)
In the model RM(m, π 0 , F ), it turns out that FDR(SUD λ (t)) has an expression only involving the Ψ k s, and not the Ψ k,j,F (Roquain and Villers, 2011a, Section 5.2). Although it has a somewhat less intuitive form, it is better than (26) from a computational point of view.

Discussion

Our aim in this paper was to address the question "is the Dirac-uniform distribution an LFC for an intermediate step-up-down procedure (that uses a standard threshold collection)?" In a nutshell, the answer we found is "no, but almost".

We provided a rigorous quantification of what "almost" means, using an alternative approach to the asymptotic results of [START_REF] Gontscharuk | Asymptotic and Exact Results on FWER and FDR in Multiple Hypothesis Testing[END_REF] that entails nonasymptotic bounds and explicit convergence rates. In practical situations, evaluating such bounds can allow to determine whether we can consider that the FDR is maximum when the signal strength is maximum.

Returning to equations ( 5) and ( 6), an additional question, particularly relevant in practice, is how appropriate it is to base the multiple type I error criterion solely on control of the expectation of the random variable FDP. We notice that Theorem 5.2 may also be used to study this issue by computing exactly the point mass function of the FDP under arbitrary configurations for the alternative, cf.

Section A in the appendix. Based on this, we investigated to what extent the distribution of the FDP concentrates around its expectation for a simple Gaussian location model with parameter µ. (i) The effect size µ is close to zero (weak signal) or (ii) the proportion π 0 of true null hypotheses is close to 1 (sparse signal). Thus, controlling the FDR alone does not guarantee a small FDP for a specific experiment at hand in these cases.

For a well-defined dependency structure induced by exchangeable test statistics, theoretical arguments for m tending to infinity support the observation that the distribution of the FDP often does not degenerate in the limit, see [START_REF] Finner | Dependency and false discovery rate: asymptotics[END_REF][START_REF] Delattre | On the false discovery proportion convergence under Gaussian equi-correlation[END_REF]. For the jointly independent case and in the cases (i) or (ii) above, this phenomenon has not been theoretically studied to the best of our knowledge. The latter can possibly be investigated by extending the asymptotic approach of [START_REF]Asymptotic properties of false discovery rate controlling procedures under independence[END_REF] to the case where µ and π 0 are allowed to depend on m.

Taking these considerations into account, control of the false discovery ex-7 Proofs 7.1 Proof of Proposition 5.1

We follow the proof of the regular Steck's recursion (Shorack and Wellner, 1986, p. 366-369). By using the convention U (0) = t 0 = 0 and by considering the smallest j for which U (j+1) > t j+1 , we can write

P(U (k) ≤ t k ) -P U (1) ≤ t 1 , . . . , U (k) ≤ t k = k-2 j=0 P(∀i ≤ j, U (i) ≤ t i , U (j+1) > t j+1 , U k ≤ t k ) = k-2 j=0 X⊂{1,...,k},|X|=j P(∀i ≤ j, U (i) ≤ t i , ∀i / ∈ X, t j+1 ≤ U i ≤ t k ).
Hence, if U (i:X) denotes the i-th smallest member of the set {U i , i ∈ X}, we obtain

P(U (k) ≤ t k ) -P U (1) ≤ t 1 , . . . , U (k) ≤ t k = k-2 j=0 X⊂{1,...,k},|X|=j P(∀i ≤ j, U (i:X) ≤ t i )P(∀i / ∈ X, t j+1 ≤ U i ≤ t k ) = k-2 j=0 j j 0 =0 X⊂{1,...,k},|X|=j
1{|X ∩ {1, . . . , k 0 }| = j 0 }Ψ j,j 0 ,F (t 1 , . . . , t j )

× P(∀i / ∈ X, t j+1 ≤ U i ≤ t k ).

Proof of Lemma 4.3

In this proof, we denote Ĝ′ m = ( Ĝm + m -1 ) ∧ 1 for short. Let us first note that for any k ∈ {1, . . . , m}, "p (k) ≤ t k " is equivalent to " Ĝm (ρ(k/m)) ≥ k/m". We now distinguish two cases:

U (λ/m, Ĝ′ m ) = min u ∈ [λ/m, 1] : Ĝ′ m (ρ(u)) ≤ u = min u ∈ {λ/m, . . . , m/m} : Ĝ′ m (ρ(u)) ≤ u , (29)
) that if G, G ′ are two nondecreasing functions such that ∀x ∈ [0, 1], G(x) ≥ G ′ (x), then U (λ, G) ≥ U (λ, G ′ ).
Based on the bound (10), we deduce that ∀δ ∈ (0, 1), As a consequence, in the DU(m, k) model, k ∈ {m 0 -1, m 0 }, by using Ĝm (x) -

sup x∈[0,1] | Ĝm (x) -G DU ζ (x)| ≤ δ - 1 m ⊂ u - δ ≤ û ≤ u + δ . λ =
G DU ζ (x) = (k/m -ζ)(1 -Ĝk (x)) -ζ( Ĝk (x) -x), we obtain Ω δ (k) := sup x∈[0,1] | Ĝk (x) -x| ≤ ζ -1 (δ -ν - 1 m ) ⊂ u - δ ≤ û ≤ u + δ . (30)
Remember that p 1 , . . . , p m 0 correspond to true nulls, hence, when k ∈ {m 0 -1, m 0 }, Ĝk involves only variables which are i.i.d. uniform. As a consequence, by using the DKW inequality with [START_REF] Massart | The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality[END_REF] optimal constant, we have in the DU(m, k) model and for k ∈ {m 0 -1, m 0 },

P DU (m,k) [(Ω δ (k)) c ] ≤ 2 exp -2k δ -ν - 1 m 2 + /ζ 2 ≤ 2 exp -2m δ -ν - 1 m 2 + (1 -ν/ζ) + , (31)
because k/m ≥ ζν and ζ ≤ 1.

Upper bound

Let q(x) = ρ(x)/x when x ∈ (0, 1] and q(0) = lim x→0 + ρ(x)/x (the limit exists in R because x ∈ (0, 1] → ρ(x)/x is non-decreasing). Applying Theorem 4.3 of Finner et al. (2009), we have

FDR(SUD λ (t), m 0 , F ) ≤ m 0 m E DU(m,m 0 -1) [q(û)] ≤ m 0 m ρ(u + δ ) u + δ + m 0 m P DU(m,m 0 -1) [(Ω δ (m 0 -1)) c ], ( 32
)
because q is non-decreasing, u + δ is positive and by (30). Next, by (31), we obtain the following upper-bound:

FDR(SUD λ (t), m 0 , F ) ≤ m 0 m ρ(u + δ ) u + δ + m 0 m 2 exp -2m δ -ν - 1 m 2 + (1 -ν/ζ) + .
(33)

Lower bound

In the model DU(m, m 0 ) with m 0 < m, we have û > 0 a.s. and thus

FDR(SUD λ (t), m 0 , F ≡ 1) = m 0 m E DU(m,m 0 ) Ĝm 0 (ρ(û)) û ≥ m 0 m E DU(m,m 0 ) Ĝm 0 (ρ(û)) û 1{Ω δ (m 0 )} ≥ m 0 m E DU(m,m 0 ) Ĝm 0 (ρ(u - δ )) u + δ 1{Ω δ (m 0 )} ≥ m 0 m ρ(u - δ ) u + δ - m 0 m 1 1 -ζ P DU(m,m 0 ) [(Ω δ (m 0 )) c ] ,
by (30) and because u + δ ≥ 1ζ. From (31), we obtain the lower-bound

FDR(SUD λ (t), m 0 , F ≡ 1) ≥ m 0 m ρ(u - δ ) u + δ - m 0 m 1 1 -ζ 2 exp -2m δ -ν - 1 m 2 + (1 -ν/ζ) + . (34)
Finally, (33) and (34) yield the result.

A Formulas for FDP distribution

From Theorem 5.2 and (25), we can compute the exact c.d.f. of the FDP of any SUD procedure in the following way, for each fixed number m ≥ 2 of hypotheses. Q m,m 0 ,F ((t λ ∨ t j ) 1≤j≤m , k, j). In particular, FDR(SD(t)) > FDR(SU(t)) for t 0 > 1 -(1m 0 /m) 1/m 0 .

The proof is straightforward and is left to the reader. As an illustration, for m = 10 and m 0 = 7, 1 -(1m 0 /m) 1/m 0 ≃ 0.158.

C DU is an LFC for the k-FWER

We state here for the sake of completeness a straightforward generalization of Lemma 1 of [START_REF] Finner | Controlling the familywise error rate with plugin estimator for the proportion of true null hypotheses[END_REF] (see also Lemma where V m is defined in (1) and φ is a function from N to R.

Assume the two following conditions are satisfied:

(i) t * is a nonincreasing function of each p-value;

(ii) φ is nondecreasing.

Then it holds that R(P, δ) ≤ R(DU(m, m 0 ), δ), and thus the conclusion.

A straightforward (though less immediately interpretable) extension of this result to procedures that are not necessarily threshold-based is to replace assumption (i) by (i'): p → V m (δ(p)) is a non increasing function of each p-value.

  true for θ and gets rejected by δ}|. (1)

  Given a statistical model, we consider a finite set of m ≥ 2 null hypotheses H 1 , . . . , H m , and a corresponding, fixed collection of tests with associated pvalue family p := (p i , i ∈ {1, ..., m}). For simplicity, we skip somewhat the formal definition of p-values and of the underlying statistical model and consider directly a statistical model for the p-values, that is, a model for the joint distribution of p.

Figure 1 :

 1 Figure 1: Value of k (vertical solid line), defined by (2), for several procedures of the SUD type. The bottom-right SUD procedure (using λ = 2m/3) coincides with the SU procedure for this realization of the p-value family.

Figure 2 :

 2 Figure 2: An arrow "A → B" means "FDR(A) ≤ FDR(B)". These results hold for the fixed mixture (FM) model except when "RM" is written. The brackets are a shortened reference to the corresponding literature, see main text for more details.

  reports some obtained values in the Gaussian case F = Φ(Φ -1 (•)µ) for different values of the alternative mean µ. We observe that the result for FM(m, m 0 , F ) (left panel) or RM(m, π 0 , F ) (right panel) are qualitatively the same: there is a range of values for λ ∈ {2, . . . , m -1} for which the FDR is larger for a smaller µ. However, this phenomenon seems to decay when m becomes larger, see Fig-ure3for m = 100. Also, when α decreases, the phenomenon still occurs but its amplitude decays. To alleviate the concern that this somewhat unexpected phenomenon could be due to numerical inaccuracies in the computation of the exact formulas (which involve several nested recursions), the reported results were double-checked via extensive and independent Monte-Carlo simulations, which confirmed the validity of the reported curves.

Fixed mixture m 0 /Figure 3 :

 03 Figure 3: The LFC of LSUD is not always DU. FDR(LSUD λ ) as a function of the order λ ∈ {1, . . . , m}. Left: fixed mixture; Right: random mixture. One sided Gaussian location model with parameter µ. α = 0.5.

Definition 4. 2 .

 2 Let ρ : [0, 1] → [0, 1] satisfying (3). For any non-decreasing function G : [0, 1] → [0, 1], and ℓ ∈ [0, 1], we define

Figure 4 :

 4 Figure 4: Illustration for u + δ and u - δ in the case of the LSUD and the SUD based on AORC. Here the X-axis is on the "threshold scale" t = ρ(u). ζ = 0.5; δ = 0.03. The area between (G DU ζ δ) ∨ 0 and (G DU ζ + δ) ∧ 1 is displayed in gray.

Corollary 4. 5 .

 5 Let α ∈ (0, 1). Consider a threshold collection of the form t k = ρ(k/m), where ρ is either ρ(x) = αx (linear) or given by (14), that is, associated with the AORC. Consider the SUD procedure of threshold collection t and of order λ = λ m possibly depending on m. Consider either the FM(m, m 0 , F ) model with m 0 = ⌊ζ m m⌋ or RM(m, π 0 , F ) model with π 0 = ζ m , for some ζ m ∈ (α, 1) possibly depending on m. Assume that 1 1-ζm (log m)/m = o(1). For the AORC case, assume moreover 1 1-λm/m (log m)/m = o(1). Then we have for any F ∈ F, (FDR(SUD λm (t), F ) -FDR(SUD λm (t), F ≡ 1)) proved in Section 7.4 and is an easy consequence of Theorem 4.4 when taking δ = δ m (and γ = γ m ) suitably tending to zero. Assuming ζ m > α is not an important restriction because when ζ m ≤ α, controlling the FDR is a trivial problem: the procedure rejecting all the hypotheses has an FDR (and even an FDP) smaller than ζ m ≤ α.While focusing on the linear and AORC rejection curve, the conclusion of Corollary 4.5 is substantially more informative than Theorem 3.1: it evaluates what is the order of the error when considering that the DU is an LFC of an SUD test. For ζ m = ζ ∈ (0, 1) fixed with m, note that the rate of convergence in (15) is equal to the parametric rate, up to a log m factor. Furthermore, the constant in the O(•) can be derived explicitly by using the bound from the previous section.

  Figure5was obtained from these exact formulas for m = 100 and varying values of π 0 and µ. Note that the unrealistically large choice of α = 1/2 has only been used for reasons of readability of the figures; similar plots also obtained when choosing α smaller (the variance of the FDP actually increases with smaller α, because this entails a smaller number of rejections). On inspection of these graphs, it becomes apparent that -even though joint independence of the p-values holds -the distribution of the FDP is not concentrated around the corresponding FDR in the following two situations:

-

  Step-up case: assume p (λ) > t λ , that is, Ĝm (ρ(λ/m)) < λ/m. Let us provethat k/m = U (λ/m, Ĝm ). Since U (ℓ, G) is a fixed point of the function G • ρ, we have U (ℓ, G) ∈ {0, 1/m, . . . , m/m}. Hence, k/m = max{k/m ∈ {0, . . . , λ/m} : Ĝm (ρ(k/m)) ≥ k/m} = max{u ∈ [0, λ/m] : Ĝm (ρ(u)) ≥ u},and we can conclude.-Step-down case: assume p (λ) ≤ t λ , that is, Ĝm (ρ(λ/m)) ≥ λ/m. First assume that k < m and prove that ( k + 1)/m = U (λ/m, Ĝ′ m ). On the one hand,( k + 1)/m = min{k/m ∈ {(λ + 1)/m, . . . , m/m} : Ĝm (ρ(k/m)) < k/m} = min{k/m ∈ {λ/m, . . . , m/m} : Ĝm (ρ(k/m)) < k/m} = min{k/m ∈ {λ/m, . . . , m/m} : Ĝ′ m (ρ(k/m)) ≤ k/m},(28)because m Ĝm (ρ(k/m)) is an integer. On the other hand, since Ĝ′ m (ρ(λ/m)) ≥ λ/m, we have

  because U (λ/m, Ĝ′ m ) ∈ {0, 1/m, . . . , m/m}. Combining (28) and (29) yields the result. Second, in the case where k = m, then for any k/m ∈ {λ/m, . . . , m/m}, we have Ĝm (ρ(k/m)) ≥ k/m. Hence, for all k/m ∈ {λ/m, . . . , (m -1)/m}, Ĝ′ m (ρ(k/m)) > k/m which entails U (λ/m, Ĝ′ m ) = 1. Hence, k/m = U (λ/m, Ĝ′ m ) in that case. Finally, the inequality k/m ≤ U (λ/m, Ĝ′ m ) always holds.7.3 Proof of Theorem 4.4Let us first prove the result in FM(m, m 0 , F ) model. We recall Ĝm (x) :=m -1 m i=1 1{p i ≤ x}, G DU ζ (x) := (1ζ)+ ζx and we put û := k/m, where k is defined by (2). We can easily check from the definition (9

Figure 6 :

 6 Figure 6: Plot of G m (ρ(•)) for ρ(x) = 0.5x and m = 10 p-values. The order λ of the SUD procedure is displayed by the dashed line while the value of k/m is displayed by the dotted line. Left: the SUD procedure of order λ = 8 is such that k/m = U (λ/m, Ĝm ) = 0.7. Right: the SUD procedure of order λ = 4 satisfies k/m = 0.7 but U (λ/m, Ĝm ) = 0.4.

Corollary A. 1 .

 1 Let λ ∈ {1, . . . , m} and consider any threshold collection t. Fix an arbitrary x ∈ (0, 1). Then the following holds:(i) in the model RM(m, π 0 , F ), for any π 0 ∈ [0, 1], F ∈ F, we have P(FDP(SUD λ (t)) ≤ x) = π 0 ,F ((t λ ∧ t j ) 1≤j≤m , k, j) π 0 ,F ((t λ ∨ t j ) 1≤j≤m , k, j);(38)(ii) in the model FM(m, m 0 , F ), for any m 0 ∈ {0, . . . , m}, F ∈ F, we haveP(FDP(SUD λ (t)) ≤ x) = m 0 ,F ((t λ ∧ t j ) 1≤j≤m , k, j)

  figuration

  2.2 of[START_REF] Gontscharuk | Asymptotic and Exact Results on FWER and FDR in Multiple Hypothesis Testing[END_REF] concerning the LFCs of multiple testing procedures under a class of type I criteria containing in particular the k-FWER (but not the FDR, as pointed out in the introduction). This result should be considered as already known by experts in the field, although we failed to locate a precise reference for it. The setting considered assumes independence of p-values corresponding to true nulls, but is more general than the fixed mixture model, since the p-values corresponding to true null hypotheses are only assumed to be stochastically larger than a uniform variable on [0,1]; also, the p-values corresponding to alternatives are not assumed to be identically distributed nor independent. Lemma C.1. Let m ≥ 1 and m 0 ∈ {0, . . . , m} be fixed. Let p = (p 1 , . . . , p m ) be a family p-values with distribution by P such that (p i ) 1≤i≤m 0 form an independent family of variables, each stochastically lower bounded by a uniform variable. Assume that δ is a multiple testing procedure rejecting all hypotheses having pvalue less than a data-dependent threshold t * (p). Let R be a type I error criterion taking the form R(P, δ) = E p∼P [φ(V m (δ(p)))],

  that is, DU(m, m 0 ) is an LFC for δ among the set of distributions satisfying the properties described above.Proof. Using (i) and (ii) together entails that p → φ(V m (δ(p))) is a nonincreasing function of each p-value. Denote p 0 = (p 1 , . . . , p m 0 , 0, . . . , 0) the p-value family obtained by replacing p i by 0 for i > m 0 , and P 0 the distribution of p 0 when p has distribution P . Obviously we haveE p∼P [φ(V m (δ(p)))] ≤ E p∼P [φ(V m (δ(p 0 )))] = E p∼P 0 [φ(V m (δ(p)))]Now applying Lemma A.11 as cited by[START_REF] Gontscharuk | Asymptotic and Exact Results on FWER and FDR in Multiple Hypothesis Testing[END_REF], we obtainE p∼P 0 [φ(V m (δ(p)))] ≤ E p∼DU(m,m 0 ) [φ(V m (δ(p)))] ,

  and explained below. These results hold either under the fixed mixture model or the random mixture model, hence involve a maximization over distributions where the p-values are independent. (While the present paper is focused on this setting, let us mention here briefly, that LFCs for the FDR criterion under arbitrary dependencies have also been studied, see, e.g.[START_REF] Lehmann | Generalizations of the familywise error rate[END_REF][START_REF] Guo | On control of the false discovery rate under no assumption of dependency[END_REF] First, let us consider the problem of the monotonicity of FDR(SUD

	λ = 1) and for a RM(m, π 0 , F ) model, Theorem 4.1 of Roquain and Villers
	(2011a) states that the Dirac-uniform configuration (F ≡ 1) is an LFC under some complex condition on the threshold collection t, that is fulfilled by the
	linear threshold collection t k = αk/m, α ∈ (0, 1) and over the class of concave c.d.f.'s. However, for λ / ∈ {1, m} (i.e., an "intermediate" SUD procedure), finding LFC's is more delicate and the only known result is asymptotic, as m tends to
	infinity. Precisely, combining Theorem 4.3 of Finner et al. (2009) and Lemma 3.7

λ (t)) in λ (vertical arrows). Recently, it was proved that, whenever F is concave, the FDR grows as the rejection set grows

(Zeisel, Zuk, and Domany, 2011, Theorem 4.1)

.

Interestingly, the rejection set R can have a very general form: the only condition is that |R| is a measurable function of the order statistics of the family of p-values under consideration. From (2) and since for any λ ∈ {1, . . . , m -1}, the rejection set of SUD λ (t) is included in the one of SUD λ+1 (t), we obtain that for a concave

F , FDR(SUD λ (t)) ≤ FDR(SUD λ+1 (t)),

both for FM(m, m 0 , F ) and RM(m, π 0 , F ) models. This implies in particular that FDR(SD(t)) ≤ FDR(SU(t)) for a concave F . Other studies establish similar inequalities, but with a condition on the threshold collection t, not on F .

Precisely, Theorem 4.3 of

Finner et al. (2009)

and Theorem 3.10 of Gontscharuk (2010) establish that, when k → t k /k is nondecreasing, for any λ ∈ {1, . . . , m-1}, FDR(SUD λ (t)) ≤ FDR(SU(t)), both for FM(m, m 0 , F ) and RM(m, π 0 , F ) models. In particular, the fact that FDR(SU(t)) dominates the FDR of SD(t) is quite well established in multiple testing literature. Nevertheless, let us stress that this is no longer the case for "atypical" configurations of F and t, as we state in Appendix B.

Secondly, let us consider the monotonicity of FDR(SUD λ (t), F ) in F . In the step-up case (i.e., λ = m), the situation is somewhat simple: Theorem 5.3 of

[START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]

states that F ≤ F ′ implies FDR(SUD λ (t), F ) ≤ FDR(SUD λ (t), F ′ ) whenever k → t k /k is nondecreasing. Moreover, the inequality is reversed whenever k → t k /k is nonincreasing. In the step-down case (i.e.,

Acknowledgements

This work was supported by the French Agence Nationale de la Recherche (ANR grant references: ANR-09-JCJC-0027-01, ANR-PARCIMONIE, ANR-09-JCJC-0101-01) and the French ministry of foreign and european affairs (EGIDE -PRO-COPE project number 21887 NJ).

ceedance (i.e., of the probability that the FDP exceeds a given threshold) has recently been proposed in the literature (see, e.g. [START_REF] Farcomeni | A review of modern multiple hypothesis testing, with particular attention to the false discovery proportion[END_REF], for a review).

Controlling the false discovery exceedance control again brings forth the question of the corresponding LFC: are Dirac-uniform configurations least favorable for, e.g., P(FDP(LSU) > x)? Some non-reported figures show that this is not the case for any x. Hence, finding LFCs for the false discovery exceedance stays an open avenue for future research.

Proof for random mixture model

In the RM(m, π 0 , F ) model with π 0 = ζ, the distribution of m 0 is binomial with parameters (m, ζ). In particular, ν is random. However, we can write

Additionally, using [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF] inequality, we can write

Combining ( 35) and ( 36) with ( 12) finishes the proof.

Proof of Corollary 4.5

First consider the FM(m, m 0 , F ) model with m 0 < m (the case

so that e -2m(δm-νm-1/m) 2 + (1-νm/ζm) + = 1/m for large m. Since ν m ≤ 2/m by assumption, we have δ m ∝ (log m)/m. From Theorem 4.4, it is sufficient to prove

From Section 4.2, this holds for the linear critical value function. This also holds for the AORC as soon as λ m /m < v δm , which is the case for large m by assumption.

The proof in the RM(m, π 0 , F ) model is similar by taking additionally γ m ∝ (log m)/m such that 2e -2m(γm-1/m) 2 + = 1/m.


David Vanderhaeghe

Pascal Barla

William Baxter

Google

Dynamic Stylized Shading Primitives

Keywords: CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Color, shading, shadowing, and texture; non-photorealistic rendering, stylized shading, shape depiction, real-time rendering

Shading appearance in illustrations, comics and graphic novels is designed to convey illumination, material and surface shape characteristics at once. Moreover, shading may vary depending on different configurations of surface distance, lighting, character expressions, timing of the action, to articulate storytelling or draw attention to a part of an object. In this paper, we present a method that imitates such expressive stylized shading techniques in dynamic 3D scenes, and which offers a simple and flexible means for artists to design and tweak the shading appearance and its dynamic behavior. The key contribution of our approach is to seamlessly vary appearance by using a combination of shading primitives that take into account lighting direction, material characteristics and surface features. We demonstrate their flexibility in a number of scenarios: minimal shading, comics or cartoon rendering, glossy and anisotropic material effects; including a variety of dynamic variations based on orientation, timing or depth. Our prototype implementation combines shading primitives with a layered approach and runs in real-time on the GPU.

Introduction

Shading in computer graphics is generally designed to simulate reality by using a combination of complex materials and lighting environments. Although this approach has many applications in architecture or special effects for movies, other applications rather use shading in a stylized manner. Examples of stylized shading choices abound in artistic illustrations and graphic novels [START_REF] Hogarth | Dynamic Light and Shade[END_REF][START_REF] Mccloud | Understanding Comics: The Invisible Art[END_REF]. Even if such illustrations may make use of very different media such as watercolor, oil paint, acrylic or pencil for instance, they all tend to reproduce shading in similar respects. Firstly, they are not constrained by physical accuracy: a few simple shading gradients are enough to produce a convincing appearance, and sharp color transitions are often used to create cartoon effects. More importantly, they convey illumination, material and surface shape characteristics all at once. Shading is thus often used to articulate storytelling, by exaggerating facial expressions through the enhancement of reflections and surface details, or by accentuating character silhouettes via sharp rim lighting.

There is another advantage in relaxing physical constraints though, as illustrated in the work of [START_REF] Hogarth | Dynamic Light and Shade[END_REF]: it allows artists to create expressive compositions, whereby shading takes on different appearances with different configurations of surface distance, * e-mail: vdh@irit.fr † Work completed in part at OLM Digital, Inc. lighting, character expressions, timing of the action, etc. This is especially noticeable in graphic novels where, from one frame to another, shading appearance may change significantly.

Minimal

Turning hand-made artworks that exhibit stylized shading appearance into computer-based animations is not trivial. An important requirement is to provide artists with a direct control over surface shape depiction, material characteristics, lighting directions and stylization effects through time. An obvious approach would be to let users paint and animate shading directly onto 3D objects, but it would require a tremendous amount of hand-tuning and keyframing. The alternative solution that we investigate in this paper is to automate common shading behaviors, while providing sufficient access to artistic control. With this approach, dynamic scene changes raise the additional challenge of seamlessly transforming a shading appearance into another in real-time. As explained in Section 2, existing shading techniques are often limited in the handling of dynamic shading changes; moreover, they are usually restricted to specific stylization effects.

In contrast, our goal is to provide a single approach to create various dynamic stylized shading effects. Our solution is to decompose shading into a set of procedural primitives amenable to dynamic behaviors. Each primitive is controlled by a single light source, and its appearance is determined by material characteristics, stylization choices, and view-dependent surface features. We make use of a single shading primitive model in order to permit continuous changes of material, shape and stylization effects. The design choices that have lead to our model are presented and discussed throughout Section 3. Primitives are evaluated and composited in real-time using the GPU, and a few primitives are sufficient to mimic most common shading styles, as demonstrated in Section 4. We discuss their benefits, limitations and future extensions in Section 5.

Previous work

One of the earliest forms of stylized shading is Toon Shading, which simply segments the result of a diffuse shading function into a few color bands. The Technical Illustration Shader of [START_REF] Gooch | A non-photorealistic lighting model for automatic technical illustration[END_REF] uses a similar approach with an unclamped diffuse term (also called Half-Lambertian), and cool-to-warm color gradients. Alternative color gradients stored in 1D textures have been used in video games [START_REF] Mitchell | Illustrative rendering in Team Fortress 2[END_REF]] to reproduce illustrative shading styles.

More complex effects can be obtained with 2D color gradients. The Lit Sphere of [START_REF] Sloan | The lit sphere: A model for capturing NPR shading from art[END_REF] consists of a single shaded sphere stored in a 2D texture and looked up using screen-space normals.

Variations on this basic idea have been successfully used in different contexts such as volumetric rendering [Bruckner and Gröller 2007] or digital sculpture software (e.g., ZBrush R or MudBox R ).

Although the Lit Sphere provides a very natural approach to the creation of stylized shading, it produces a relatively static shading appearance, since the lighting direction is implicitly "baked" in the 2D texture. The X-Toon method of [START_REF] Barla | X-Toon: An extended toon shader[END_REF] takes a different approach where the 1D toon function is extended to a 2D function, stored in a texture. The additional dimension is used to create effects that vary with depth or surface orientation for instance. Although light direction may vary in this approach, it permits only one effect that is baked in each 2D texture.

An alternative approach to stylized shading is to take advantage of inverse rendering techniques that allow direct material and illumination editing. The Illumination Brush of [START_REF] Okabe | Illumination Brush: Interactive Design of All-Frequency Lighting[END_REF] permit such control: starting from a known BRDF, they infer the environment lighting by directly painting on the object the desired result of either diffuse or specular reflections. Even if the paint strokes may be moved around for stylization purpose, the adaptation of the method to the creation of stylized shading is not straightforward, since it relies on BRDF models and realistic illumination.

We rather seek a method that is not bound by such physical constraints.

One reason for providing additional freedom is to give artists the ability to control how shape is conveyed through shading. Previous techniques have proposed extensions to existing shading models that enhance surface characteristics. The Exaggerated shading technique of [START_REF] Rusinkiewicz | Exaggerated Shading for Depicting Shape and Detail[END_REF] locally aligns light directions with grazing surface angles so that details are revealed through variations of a Half-Lambertian shader. Downsides of this approach are that it tends to produce "flat" appearances where the main illumination direction is lost, and visual artifacts may occur during light motion. The Apparent Relief technique [START_REF] Vergne | Apparent relief: A shape descriptor for stylized shading[END_REF]] rather exploits the additional dimension of the X-Toon shader to convey shape features computed via a combination of object-and image-space measurements. The method is able to depict shape features at varying scales, but exhibits artifacts with distant objects due to the combination of surface measurements. Both methods are restricted to a single shading style though. This is not the case with the Radiance Scaling technique [START_REF] Vergne | Radiance scaling for versatile surface enhancement[END_REF]], which works with arbitrary shading models: it consists in scaling incoming radiance according to surface and material characteristics. However, because incoming light directions are treated independently, stylized shading gradients are likely to be disrupted.

Apart from the simplest stylized shading techniques, the methods presented so far do not permit any dynamic control over stylized material characteristics. The cartoon highlights of [START_REF] Anjyo | Stylized Highlights for Cartoon Rendering and Animation[END_REF]2006] deal with the specific case of the Blinn-Phong specular term. Highlight shape is tweaked via translation, rotation, scaling, splitting, squaring, and Boolean operators applied to the half vector. Artists can control highlight shape deformations through keyframing. [START_REF] Pacanowski | Sketch and Paint-based Interface for Highlight Modeling[END_REF] propose a more direct interface whereby artists draw the 2D profile of a highlight in a plane perpendicular to the reflection direction. In their system, highlight shape can be deformed based on a few light directions to control shading appearance at grazing angles. The two methods focus essentially on highlight shape, whereas other material characteristics are left unchanged, and shape features are ignored.

In contrast to previous work, our approach provides a single solution to the stylized shading of 3D scenes where both surface features and material characteristics are conveyed dynamically.

Shading model

In our approach, stylized shading is designed as a combination of a few procedural shading primitives (Section 3.1) whose parameters are dynamically modified when applied to 3D objects (Section 3.2). Parameters of our shading model are summarized in Table 1.

Symbol Domain Name K(•) [0, π] → [0, 1] 3 color α [0, 1] specularity τ [-π, π] extent f [0, π] intensity fall-off c [f, π] intensity cut-off λ (-1, 1) material anisotropy µ R surface enhancement χ R concave/convex transition Table 1: List of shading parameters.

Shading primitives

Three principal requirements have led to the design of our shading primitives. First, a single primitive model should account for different shading behaviors, from diffuse to specular. Second, continuous changes of primitive parameters should lead to continuous (and ideally perceptually uniform) visual changes. Third, a primitive should be easily manipulated directly on 3D models, and in particular it should not require the tweaking of too many parameters. With these technical goals in mind, we define a shading primitive as composed of three basic ingredients: an angular parametrization u, and a pair of color and intensity profile functions, defined on this parametrization and noted K and I respectively.

The parametrization u :

S 2 × S 2 × S 2 → [0, π] is defined by u(n n n, l l l, v v v) = ⌊acos(d d dα • l l l) -τ ⌋, (1)
where n n n, l l l and v v v correspond to the normal, light and view vectors, d d dα is a reference direction that controls how the parametrization evolves when l l l is rotated, and ⌊x⌋ is a function that clamps x to [0, π]. The user-specified parameter τ ∈ [0, π] permits to control the extent of the shading primitive independently of its color and intensity profiles. The reference direction d d dα is given by

d d dα = (1 -α) n n n + α r r r ||(1 -α) n n n + α r r r|| , (2)
where r r r = 2(v v v •n n n) n n n-v v v is the view vector reflected around the surface normal, and α is a user-controlled parameter that interpolates between diffuse (α = 0) and specular (α = 1) shading behaviors. Different choices of parametrization are illustrated in Figure 2: we set up an orthographic view and render a sphere that represents the entire set of normal directions for a particular light direction l l l. The angular parametrization, displayed with a color code and isolines, evolves continuously when l l l, τ or α are manipulated. Once a parametrization has been chosen, it becomes easier to specify color and intensity profiles. In our system, the color profile K(u) is simply defined as a smooth 1D color gradient (commonly called a color ramp), while intensity is controlled by a more complex profile function. This choice is motivated by the observation that sharp transitions due to either shininess or cartoon effects are more easily controlled through a scalar function.

The intensity profile I : [0, π] → [0, 1] is defined by

I(u) = |β + (1 -β) cos u)| γ + , (3)
where β ∈ [-1/2, 1/2] is a bias parameter that permits to extend the primitive intensity toward the interval [π/2, π], γ ∈ R + is an exponent parameter that determines the intensity fall-off rate, and |x|+ is a function that clamps negative values of x to 0. Figure 3 shows how different values of β and γ allow users to control shading appearance for a fixed parametrization.

The bias and exponent parameters are inspired respectively by the Half-Lambert and Phong specular terms, which have been frequently used in previous work on stylized shading. In our system though, bias and exponent are employed in all situations, not only because we use a single primitive model for both diffuse and specular effects, but also because it provides additional creative freedom. This flexibility has one minor drawback though: specular primitives produce exaggerated rim effects when the light direction comes from the back of the object because our formulation does not incorporate the geometric term (n n n•l l l) in this case. We retain this rim effect because it offers an interesting specular style, but the original Phong-like behavior may optionally be re-introduced by multiplying I(u) by (1 -α) + α(n n n • l l l) (see Supplemental material).

Uniformly changing bias and exponent values leads to sudden changes in shading appearance as shown for γ at the top of Figure 4. This is not only a matter of user interface; it also raises issues during interpolation such as when applying the dynamic shading variations presented in Section 4. A better choice of parameters are the fall-off f ∈ [0, π] and cut-off c ∈ [f, π] angles at which the intensity profile respectively exhibits an inflection and reaches 0. For instance, uniformly changing f values leads to a more uniform interpolation as shown at the bottom of Figure 4. We explain how to compute bias and exponent based on these parameters in Equations 7 and 8 of the Appendix. Classic shading behaviors are then easily reproduced: Phong's diffuse term corresponds to f = c = π/2, and the Half-lambertian term is obtained by further setting c = π; Phong's specular term (with the geometric term re-introduced) is the same as diffuse, but leaves f free to vary in the [0, π/2] range; cartoon shading occurs when c = 0 and τ > 0, and ambient shading when τ = π. In our system, τ , f and c are directly controllable on top of 3D objects via isolines.

Shading variations

When applied to a 3D object, a shading primitive is evaluated at each surface point by computing local light and view vectors as in traditional shading techniques. A single directional or point light source is required for each primitive. Shading variations automatically occur from one surface point to another when using a point light source and/or a perspective camera. However, since primitives are entirely procedural, their parameters may also be dynamically modified across the object to create intentional shading variations. We distinguish between two types of variations: those that only depend on local object characteristics (such as shape or material features), and those that depend on scene or shot characteristics (such as depth, timing of action, etc). The former is presented in the remainder of this Section while the latter is presented in Section 4 when we consider the combination of multiple of our shading primitives.

The main source of real-world local material variation is caused by anisotropic distributions of microfacets. They are typical of materials such as brushed metals or satin, and usually controlled by a user-specified tangent field defined at each surface point. To mimic such material effects, we adapt the solution of [START_REF] Anjyo | Stylized Highlights for Cartoon Rendering and Animation[END_REF]2006] to work with our primitive model. We first compute the half vector h h h between l l l and v v v, then transform it based on the local tangent t t t and bi-tangent b b b vectors as in the following:

h h h = s λ h h ht t t + 1 s λ h h h b b b + h h hn n n ||s λ h h ht t t + 1 s λ h h h b b b + h h hn n n|| (4)
where h h h t t t,b b b,n n n are the component vectors of h h h in tangent space (e.g., h h ht t t = (h h h • t t t)t t t) and s λ is a scaling factor defined by in rotating h h h toward or away from n n n, depending on its projection in the tangent plane: when λ > 0, half vectors are compressed in the t t t direction; when λ < 0, half vectors are stretched toward the t t t direction; and when λ = 0, half-vectors are kept unchanged. We systematically apply the inverse scaling, 1/s λ , in the direction of b b b to roughly maintain highlight size regardless of anisotropy as seen in Equation 4. The back-lighting term η in Equation 5 is needed to revert to an isotropic behavior (s λ = 1) when the light direction is opposite to the view direction, since in this case the half-vector is ill-defined.

s λ = 1 -η + η 1 1-λ if λ ≥ 0 1 -η + η 1 1+λ -1 if λ < 0 (
Our shading model does not produce specular behaviors based on the half-vector h h h but rather makes use of the reflected view vector r r r for interpolation purpose. We thus first reflect l l l around h h h to compute the transformed view vector ṽ

v v = 2(l l l • h h h) h h h -l l l, and then reflect ṽ v v around n n n to obtain r r r = 2(n n n • ṽ v v)n n n -ṽ v v.
As illustrated in Figure 5 where the tangent field is aligned with sphere meridians, anisotropy is ignored in the diffuse case (d d d0 = n n n), and appears progressively when tending toward a specular behavior (d d d1 = r r r).

Another important set of local object characteristics are surface features, and recent techniques have provided solutions for enhancing them through shading. We present a novel approach where we manipulate directly the shape of shading primitives instead of modifying independent shading values as in previous work. We first measure surface curvature κ at each visible surface point. Although our approach is independent of the choice of curvature measure, we have chosen the view-centered mean curvature of [START_REF] Vergne | Radiance scaling for versatile surface enhancement[END_REF] because it provides automatic levels-of-detail and is computed dynamically from screen-space normals. We then offset the parametrization function based on local surface curvature: τ = τ + µ tanh(κχ), (6) where µ ∈ R controls the magnitude of enhancement and χ ∈ R the slope of the transition between concave (κ < 0) and convex (κ > 0) features. The reason why we have chosen to vary τ is that it only modifies the extent of a primitive: for instance, shading may be "attracted" in convexities and "repelled" from concavities. In comparison, other alternatives such as Radiance Scaling modify shading appearance by disrupting existing shading gradients as shown in Figure 6.

Original

Proposed approach Radiance Scaling Figure 6: Shape depiction: primitive extent is modified to better depict surface features. Here, the shading over a blobby object (left) is modified to convey concave features using µ = 0.3 and χ = 30 (middle). Compared to Radiance Scaling (right), our solution better preserves shading gradients while still conveying shape.

We have also tried blending between a pair of color ramps inside a single primitive, using a similar transition function for concave to convex regions. This technique is inspired by the ZBrush R Mat-Caps, where two LitSpheres are blended with a concave-to-convex transition function. An advantage of our approach is that shading primitives may have different blending transitions, while in the case of MatCaps, the entire LitSpheres are blended together. However, this technique strongly modifies material appearance, and is thus only used for specific effects (see supplemental material).

Results

Dynamic stylized shading primitives are currently implemented in GLSL and we obtain real-time performance (> 60 fps on a NVIDIA GTX480) in 1600 × 1200 for all the examples in this paper. Primitive evaluation is performed in a pixel shader, and we have chosen a layered interface for combining primitives in a way similar to classic image-processing software such as Photoshop R or Gimp R . As demonstrated in all our results and the supplemental video, a few primitives combined with alpha or additive blending permit to quickly imitate existing artwork or mimic complex materials. This is not only valuable in interactive applications, but also at the compositing stage every time light-weight re-lighting is preferable to a more heavy-weight solution at the 3D rendering stage. Our shader is given in the Appendix.

The tutorial video illustrates the process of creating a stylized shading in our prototype system. It demonstrates various features: shading behaviors, isoline controls, shape-based primitive extent and color variations, primitive blending, etc. Each primitive position is controlled independently to create ambient, diffuse, highlights and rim effects. We have also found it useful to bind all light sources to a common illumination reference frame. This is the solution we have used for most of the results shown in the gallery video. Some of these shading styles are presented in Figures 1 and7. A few primitives are sufficient to obtain compelling minimal shading results. In the first example, we use a single white specular primitive deformed to convey convex shape features in a quite exaggerated way, which provides a compelling depiction of shape within a single layer. In the second example, we use a pair of opposite specular primitives deformed so that all surface details of the face are represented. The third example makes use of 3 anisotropic primitives that convey a metallic appearance to the model. The comics shading example shows a character rendered with 7 primitives and lit from three different directions. Shape features are used to align shading with anatomical details and to darken colors in concavities. Simpler primitives are used for cartoon shading styles (we have used 6 and 4 primitives for the left and right images respectively). Embossing and contouring effects are created by attaching the light source to the camera, while sharp transitions are obtained when setting c = 0 and τ > 0. More complex glossy effects are obtained by combining multiple primitives. We have used 6 primitives in both images of the dog model; they differ only in color and light source positions. The same shape-based variations have been applied to demonstrate that we obtain similar surface depictions in both cases. The Venus and Buddha images illustrate anisotropic shading styles, where we have combined 4 primitives in both cases to mimic velvet and gold.

Intentional and dynamic variations of primitive parameters prove to be useful in many situations, as shown in Figure 8. For instance, varying primitive profile and extent based on object orientation is useful in scientific illustration scenarios. In the sanguine-like shading, the reddish and whitish shading primitives are modified when seen from behind, while a third ambient primitive is kept constant. In the beetle image (made of 7 primitives), the highlight profiles and their blending are altered to draw attention to the material when on top of the insect, while the bottom surface is depicted with exaggerated shape-based deformations of a diffuse primitive. A more direct shading artistic control is obtained through key-framing, as in the elephant image which shows two frames of an animation where shape-based primitive deformations are activated each time the elephant hits the ground. Although we have applied the effect everywhere on the elephant, it could be easily localized by means of an alpha map. Finally, continuous variations are obtained by mapping any scene parameter to one or more primitive parameters.

An example is given in the landscape image where an aerial perspective effect is produced by mapping depth to shading primitive variations. Here, we alter the color of a diffuse primitive while smoothly suppressing its shape-based deformations with increasing depth; we also let another slightly specular primitive only appear in the foreground to give more volume to nearby relief. Note that using a simpler approach such as X-Toon [START_REF] Barla | X-Toon: An extended toon shader[END_REF]] would only permit one depth-based effect at a time. With our prototype system, it takes a couple minutes to get a first shading result, and up to half an hour for refining more complex styles. Although we have shown quite simple dynamic controls based on depth, orientation and keyframing, our shading primitives may be controled in more complex ways, which ultimately depend on the final application.

Discussion

We have presented an approach for dynamic stylized shading based on the composition of a few procedural primitives. Our solution is well adapted to both interactive applications such as video games, scientific illustration, or digital sculpture; and offline applications such as stylized rendering and compositing for animated films and commercials. Our shading model is designed to be continuously interpolated among a variety of stylized shading effects, such as diffuse or specular behaviors, surface and anisotropy features, etc.

A limitation of our stylized shading model is that it only produces opaque material effects, from diffuse to glossy, in direct illumination settings. We believe another approach is required for producing stylized renderings of objects with highly refractive or reflective materials. Indeed, citing [START_REF] Hogarth | Dynamic Light and Shade[END_REF] (p.120): "The only way to draw [complex refractions] is to observe the seemingly random patterns and record them faithfully." More diffuse effects such as translucency or soft inter-reflections (including color bleeding) may be taken into account by extending the range of shading behaviors. Shadows may also be used to control shading primitives, by mapping occlusion values to primitive parameters. However, we believe that specific shadowing primitives will be necessary to provide an accurate control over shadow appearance.

There are different ways in which our primitive model could be extended in future work. For instance, as demonstrated in sup-plemental material, we have reproduced highlight shapes of [START_REF] Anjyo | Stylized Highlights for Cartoon Rendering and Animation[END_REF]2006] and [START_REF] Pacanowski | Sketch and Paint-based Interface for Highlight Modeling[END_REF]. Although these extensions produce interesting styles, we have preferred retaining simple primitive shapes in our model, mainly because complex shapes do not transfer appropriately to different 3D objects when these have detailed surface features. This is especially noticeable with highlights that are disrupted in the presence of rapidly changing surface curvature. A promising avenue of future work would be to directly control primitive shape on top of 3D objects without any indirection. In particular, such an approach should allow artists to draw highlights that somehow "ignore" small curvature variations to retain their shape during camera, object or lighting motion, while still conveying material properties.

ComicsFigure 1 :

 1 Figure 1: Shading styles. Various styles are created by combining multiple shading primitives.

Figure 2 :

 2 Shading parametrization: Here we illustrate the parametrization u for the light direction shown at top left. The left most color bar presents the color code used. By varying α, different shading behaviors are continuously produced from diffuse (α = 0) to specular (α = 1). We have used an offset of τ = 0 in all cases.

β

  = 0, γ = 1 β = 0.5, γ = 1 β = 0.5, γ = 7 Figure 3: Shading intensity: by varying β and γ, different shading appearances are continuously obtained from sharp to smooth profiles. Here we use a simple diffuse shading parametrization (α = 0, τ = 0) and a smooth color ramp in the right image.

Figure 4 :

 4 Figure4: Profile interpolation: changing the exponent γ from 1 to 30 by uniform increments exhibits non-uniform shading appearance interpolation (top row). In contrast, with uniform increments of f , keeping c = π/2, we obtain an interpolation behavior that is more uniform (bottom row).

Figure 5 :

 5 5)with λ ∈ (-1, 1) a user-controlled anisotropy parameter, and η = (l l l • v v v + 1)/2 a back-lighting term. Intuitively, Equation 4 consists Shading anisotropy: the λ parameter controls anisotropy (here we use λ = -0.5) and is progressivelly introduced when the shading behavior tends toward specular.

Figure 7 :

 7 Figure 7: Shading styles. Various styles are created by combining multiple shading primitives.

Acknowledgments

This work has been sponsored by the Animare (ANR-08-JCJC-0078-01), SeARCH (ANR-09-CORD-019) projects and the INRIA postdoctoral program.

The cutoff parameter c corresponds to the angle at which the profile intensity function reaches zero. The bias function β(c) gives the bias value that satisfies this condition: I(c) = 0. It is defined by

The falloff parameter f locates the angle at which the profile function passes through an inflection point, which corresponds to the zero-crossing of its second derivative. The exponent function γ(f, c) gives the (strictly positive) exponent value that satisfies this condition: I ′′ (f ) = 0. The profile's second derivative is:

There is only one strictly positive zero-crossing for I ′′ (f ):

Replacing β by the formula for the bias function β(c) gives:

GLSL shader to compute a primitive contribution. We use a constant color instead of color ramps for simplicity.