Christophe Giraud

MIXING LEAST-SQUARES ESTIMATORS WHEN THE VARIANCE IS UNKNOWN

Keywords: Mathematics Subject Classification, 62G08, Key words and phrases, Gibbs mixture -shrinkage estimator -oracle inequalities -adaptive estimationlinear regression, 1

We propose a procedure to handle the problem of Gaussian regression when the variance is unknown. We mix least-squares estimators from various models according to a procedure inspired by that of Leung and Barron [17]. We show that in some cases the resulting estimator is a simple shrinkage estimator. We then apply this procedure in various statistical settings such as linear regression or adaptive estimation in Besov spaces. Our results provide non-asymptotic risk bounds for the Euclidean risk of the estimator.

Introduction

We consider the regression framework, where we have noisy observations [START_REF] Akaike | Statistical predictor identification[END_REF] Y i = µ i + σε i , i = 1, . . . , n of an unknown vector µ = (µ 1 , . . . , µ n ) ′ ∈ R n . We assume that the ε i 's are i.i.d standard Gaussian random variables and that the noise level σ > 0 is unknown. Our aim is to estimate µ.

In this direction, we introduce a finite collection {S m , m ∈ M} of linear spaces of R n , which we call henceforth models. To each model S m , we associate the least-squares estimator μm = Π Sm Y of µ on S m , where Π Sm denotes the orthogonal projector onto S m . The L 2 -risk of the estimator μm with respect to the Euclidean norm

• on R n is (2) E µ -μm 2 = µ -Π Sm µ 2 + dim(S m )σ 2 .
Two strategies have emerged to handle the problem of the choice of an estimator of µ in this setting. One strategy is to select a model S m with a data driven criterion and use μ m to estimate µ. In the favorable cases, the risk of this estimator is of order the minimum over M of the risks [START_REF] Baraud | Gaussian model selection with unknown variance[END_REF]. Model selection procedures have received a lot of attention in the literature, starting from the pioneer work of Akaike [START_REF] Akaike | Statistical predictor identification[END_REF] and Mallows [START_REF] Mallows | Some comments on Cp[END_REF]. It is beyond the scope of this paper to make an historical review of the topic. We simply mention in the Gaussian setting the papers of Birgé and Massart [START_REF] Birgé | Gaussian model selection[END_REF][START_REF] Birgé | A generalized Cp criterion for Gaussian model selection[END_REF] (influenced by Barron and Cover [START_REF] Barron | Minimum complexity density estimation[END_REF] and Barron, Birgé and Massart [START_REF] Barron | Risk bounds for model selection via penalization[END_REF]) which give non-asymptotic risk bounds for a selection criterion generalizing Mallows'C p .

An alternative to model selection is mixing. One estimates µ by a convex (or linear) combination of the μm s

(3) μ = m∈M w m μm , with weights w m which are σ(Y )-measurable random variables. This strategy is not suitable when the goal is to select a single model S m, nevertheless it enjoys the nice property that μ may perform better than the best of the μm s. Various choices of weights w m have been proposed, from an information theoretic or Bayesian perspective. Risk bounds have been provided by Catoni [START_REF] Catoni | Universal aggregation rules with exact bias bounds[END_REF], Yang [START_REF] Yang | Combining different procedures for adaptive regression[END_REF][START_REF] Yang | Combining forecasting procedures: some theoretical results[END_REF], Tsybakov [START_REF] Tsybakov | Optimal rates of aggregation[END_REF] and Bunea et al. [START_REF] Bunea | Aggregation for Gaussian regression[END_REF] for regression on a random design and by Barron [3], Catoni [START_REF] Catoni | Mixture approach to universal model selection[END_REF] and Yang [START_REF] Yang | Mixing Strategies for Density Estimation[END_REF] for density estimation. For the Gaussian regression framework we consider here, Leung and Barron [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF] propose a mixing procedure for which they derive a precise non-asymptotic risk bound. When the collection of models is not too complex, this bound shows that the risk of their estimator μ is close to the minimum over M of the risks [START_REF] Baraud | Gaussian model selection with unknown variance[END_REF]. Another nice feature of their mixing procedure is that both the weights w m and the estimators μm are build on the same data set, which enable to handle cases where the law of the data set is not exchangeable. Unfortunately, their choice of weights w m depends on the variance σ 2 , which is usually unknown.

In the present paper, we consider the more practical situation where the variance σ 2 is unknown.

Our mixing strategy is akin to that of Leung and Barron [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF], but is not depending on the variance σ 2 . In addition, we show that both our estimator and the estimator of Leung and Barron are simple shrinkage estimators in some cases. From a theoretical point of view, we relate our weights w m to a Gibbs measure on M and derive a sharp risk bound for the estimator μ. Roughly, this bound says that the risk of μ is close to the minimum over M of the risks [START_REF] Baraud | Gaussian model selection with unknown variance[END_REF] in the favorable cases. We then discuss the choice of the collection of models {S m , m ∈ M} in various situations. Among others, we produce an estimation procedure which is adaptive over a large class of Besov balls.

Before presenting our mixing procedure, we briefly recall that of Leung and Barron [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF]. Assuming that the variance σ 2 is known, they use the weights (4)

w m = π m Z exp -β Y -μm 2 /σ 2 + 2dim(S m ) -n , m ∈ M
where {π m , m ∈ M} is a given prior distribution on M and Z normalizes the sum of the w m s to one. These weights have a Bayesian flavor. Indeed, they appear with β = 1/2 in Hartigan [START_REF] Hartigan | Bayesian regression using Akaike priors[END_REF] which considers the Bayes procedure with the following (improper) prior distribution: pick an m in M according to π m and then sample µ "uniformly" on S m . Nevertheless, in Leung and Barron [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF] the role of the prior distribution {π m , m ∈ M} is to favor models with low complexity. Therefore, the choice of π m is driven by the complexity of the model S m rather than from a prior knowledge on µ. In this sense their approach differs from the classical Bayesian point of view. Note that the term Y -μm 2 /σ 2 + 2dim(S m ) -n appearing in the weights ( 4) is an unbiased estimator of the risk (2) rescaled by σ 2 . The size of the weight w m then depends on the difference between this estimator of the risk (2) and -log(π m ), which can be thought as a complexity-driven penalty (in the spirit of Barron and Cover [5] or Barron et al. [4]). The parameter β tunes the balance between this two terms. For β ≤ 1/4, Theorem 5 in [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF] provides a sharp risk bound for the procedure.

The rest of the paper is organized as follows. We present our mixing strategy in the next section and express in some cases the resulting estimator μ as a shrinkage estimator. In Section 3, we state non-asymptotic risk bounds for the procedure and discuss the choice of the tuning parameters. Finally, we propose in Section 4 some weighting strategies for linear regression or for adaptive regression over Besov balls. Section 5 is devoted to a numerical illustration and Section 6 to the proofs. Additional results are given in the Appendix.

We end this section with some notations we shall use along this paper. We write |m| for the cardinality of a finite set m, and < x, y > for the inner product of two vectors x and y in R n . To any real number x, we denote by (x) + its positive part and by ⌊x⌋ its integer part.

The estimation procedure

We assume henceforth that n ≥ 3.

2.1. The estimator. We start with a finite collection of models {S m , m ∈ M} and to each model S m we associate the least-squares estimator μm = Π Sm Y of µ on S m . We also introduce a probability distribution {π m , m ∈ M} on M, which is meant to take into account the complexity of the family and favor models with low dimension. For example, if the collection {S m , m ∈ M} has (at most) e ad models per dimension d, we suggest to choose π m ∝ e (a+1/2)dim(Sm) , see the example at the end of Section 3.1. As mentioned before, the quantity -log(π m ) can be interpreted as a complexity-driven penalty associated to the model S m (in the sense of Barron et al. [START_REF] Barron | Risk bounds for model selection via penalization[END_REF]). The performance of our estimation procedure depends strongly on the choice of the collection of models {S m , m ∈ M} and the probability distribution {π m , m ∈ M}. We detail in Section 4 some suitable choices of these families for linear regression and estimation of BV or Besov functions.

Hereafter, we assume that there exists a linear space S * ⊂ R n of dimension d * < n, such that S m ⊂ S * for all m ∈ M. We will take advantage of this situation and estimate the variance of the noise by

(5) σ2 = Y -Π S * Y 2 N *
where N * = n -d * . We emphasize that we do not assume that µ ∈ S * and the estimator σ2 is (positively) biased in general. It turns out that our estimation procedure does not need a precise estimation of the variance σ 2 and the choice (5) gives good results. In practice, we may replace the residual estimator σ2 by a difference-based estimator (Rice [START_REF] Rice | Bandwidth choice for nonparametric kernel regression[END_REF], Hall et al. [START_REF] Hall | Asymptotically optimal differencebased estimation of variance in nonparametric regression[END_REF], Munk et al. [START_REF] Munk | On difference based variance estimation in nonparametric regression when the covariate is high dimensional[END_REF], Tong and Wang [START_REF] Tong | Estimating residual variance in nonparametric regression using least squares[END_REF], Wang et al. [START_REF] Wang | Effect of mean on variance function estimation in nonparametric regression[END_REF], etc) or by any non-parametric estimator (e.g. Lenth [START_REF] Lenth | Quick and easy analysis of unreplicated factorials[END_REF]), but we are not able to prove any bound similar to [START_REF] Devore | Constructive approximation. Grundlehren der Mathematischen Wissenschaften[END_REF] or (13) when using one of these estimators.

Finally, we associate to the collection of models {S m , m ∈ M}, a collection {L m , m ∈ M} of non-negative weights. We recommend to set L m = dim(S m )/2, but any (sharp) upper bound of this quantity may also be appropriate, see the discussion after Theorem 1. Then, for a given positive constant β we define the estimator μ by

(6) μ = m∈M w m μm with w m = π m Z exp β μm 2 σ2 -L m ,
where Z is a constant that normalizes the sum of the w m s to one. An alternative formula for w m is

w m = π m exp -β Π S * Y -μm 2 /σ 2 -L m /Z ′ with Z ′ = e -β Π S * Y 2 /σ 2 Z.
We can interpret the term Π S * Y -μm 2 /σ 2 + L m /β appearing in the exponential as a (biased) estimate of the risk (2) rescaled by σ 2 . As in (4), the balance in the weight w m between this estimate of the risk and the penalty -log(π m ) is tuned by β. We refer to the discussion after Theorem 1 for the choice of this parameter. We mention that the weights {w m , m ∈ M} can be viewed as a Gibbs measure on M and we will use this property to assess the performance of the procedure.

We emphasize in the next section, that μ is a simple shrinkage estimator in some cases.

2.2.

A simple shrinkage estimator. In this section, we focus on the case where M consists of all the subsets of {1, . . . , p}, for some p < n and S m = span {v j , j ∈ m} with {v 1 , . . . , v p } an orthonormal family of vectors in R n . We use the convention S ∅ = {0}. An example of such a setting is given in Section 4.2, see also the numerical illustration Section 5. Note that S * corresponds here to S {1,...,p} and d * = p.

To favor models with small dimensions, we choose the probability distribution [START_REF] Birgé | Gaussian model selection[END_REF]

π m = 1 + 1 p α -p p -α|m| , m ∈ M,
with α > 0. We also set L m = b|m| for some b ≥ 0.

Proposition 1. Under the above assumptions, we have the following expression for μ

(8) μ = p j=1 (c j Z j )v j , with Z j =< Y, v j > and c j = exp βZ 2 j /σ 2 p α exp (b) + exp βZ 2 j /σ 2 .
The proof of this proposition is postponed to Section 6.1. The main interest of Formula ( 8) is to allow a fast computation of μ. Indeed, we only need to compute the p coefficients c j instead of the 2 p weights w m of formula [START_REF] Birgé | An adaptive compression algorithm in Besov spaces[END_REF].

The coefficients c j are shrinkage coefficients taking values in [0, 1]. They are close to one when Z j is large and close to zero when Z j is small. The transition from 0 to 1 occurs when Z 2 j ≈ β -1 (b+α log p)σ 2 . The choice of the tuning parameters α, β and b will be discussed in Section 3.2. for which the c j s are given by

c j = 1 0 q k =j q + (1 -q) exp -βZ 2 k /σ 2 + b dq 1 0 p k=1 q + (1 -q) exp -βZ 2 k /σ 2 + b dq , for j = 1, . . . , p.
This formula can be derived from the Appendix of Leung and Barron [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF].

Remark 2. When the variance is known, we can give a formula similar to [START_REF] Birgé | A generalized Cp criterion for Gaussian model selection[END_REF] for the estimator of Leung and Barron [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF]. Let us consider the same setting, with p ≤ n. Then, when the distribution {π m , m ∈ M} is given by ( 7), the estimator (3) with weights w m given by (4) takes the form

(9) μ = p j=1 e βZ 2 j /σ 2
p α e 2β + e βZ 2 j /σ 2 Z j v j .

3. The performance

3.1.

A general risk bound. The next result gives an upper bound on the L 2 -risk of the estimation procedure. We remind the reader that n ≥ 3 and set (10)

φ : ]0, 1[ → ]0, +∞[ x → (x -1 -log x)/2 ,
which is decreasing. with φ defined by [START_REF] Catoni | Mixture approach to universal model selection[END_REF]. Assume also that L m ≥ dim(S m )/2, for all m ∈ M. Then, we have the following upper bounds on the L 2 -risk of the estimator μ

E µ -μ 2 ≤ -(1 + ε n ) σ2 β log m∈M π m e -β[ µ-Π Sm µ 2 -dim(Sm)σ 2 ]/σ 2 -Lm + σ 2 2 log n (12) ≤ (1 + ε n ) inf m∈M µ -Π Sm µ 2 + σ2 β (L m -log π m ) + σ 2 2 log n , ( 13
)
where ε n = (2n log n) -1 and σ2 = σ 2 + µ -Π S * µ 2 /N * .
The proof Theorem 1 is delayed to Section 6.3. Let us comment this result.

To start with, the Bound (12) may look somewhat cumberstone but it improves [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] when there are several good models to estimate µ. For example, we can derive from [START_REF] Devore | Constructive approximation. Grundlehren der Mathematischen Wissenschaften[END_REF] the bound

E µ -μ 2 ≤ (1 + ε n ) inf m∈M µ -Π Sm µ 2 + σ2 β (L m -log π m ) + inf δ≥0 δ - σ2 β log |M δ | + σ 2 2 log n , where M δ is the set made of those m * in M fulfilling µ -Π S m * µ 2 + σ2 β (L m * -log π m * ) ≤ δ + inf m∈M µ -Π Sm µ 2 + σ2 β (L m -log π m ) .
In the extreme case where all the quantities µ - 12) then improves ( 13) by a factor β -1 σ2 log |M|.

Π Sm µ 2 + σ2 β (L m -log π m ) are equal, (
We now discuss the choice of the parameter β and the weights {L m , m ∈ M}. The choice L m = dim(S m )/2 seems to be the more accurate since it satisfies the conditions of Theorem 1 and minimizes the right hand side of ( 12) and ( 13). We shall mostly use this one in the following, but there are some cases where it is easier to use some (sharp) upper bound of the dimension of S m instead of dim(S m ) itself, see for example Section 4.3.

The largest parameter β fulfilling Condition ( 11) is ( 14)

β = 1 4 φ -1 log n N * -2 < 1 4 .
We suggest to use this value, since it minimizes the right hand side of ( 12) and ( 13). Nevertheless, as discussed in Section 3.2 for the situation of Section 2.2, it is sometimes possible to use larger values for β.

Finally, we would like to compare the bounds of Theorem 1 with the minimum over M of the risks given by (2). Roughly, the Bound [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] states that the estimator μ achieves the best trade-off between the bias µ -Π Sm µ 2 /σ 2 and the complexity term C m = L m -log π m . More precisely, we derive from ( 13) the (cruder) bound

(15) E µ -μ 2 ≤ (1 + ε n ) inf m∈M µ -Π Sm µ 2 + 1 β C m σ 2 + R * n σ 2 , with ε n = 1 2n log n and R * n = 1 2 log n + µ -Π S * µ 2 βN * σ 2 sup m∈M C m .
In particular, if C m is of order dim(S m ), then [START_REF] Hartigan | Bayesian regression using Akaike priors[END_REF] allows to compare the risk of μ with the infimum of the risks (2). We discuss this point in the following example.

Example. Assume that the family M has an index of complexity (M, a), as defined in [START_REF] Baraud | Gaussian model selection with unknown variance[END_REF], which means that

| {m ∈ M, dim(S m ) = d} | ≤ M e ad , for all d ≥ 1.
If we choose, example given, ( 16)

π m = e -(a+1/2)dim(Sm)
m ′ ∈M e -(a+1/2)dim(S m ′ ) and L m = dim(S m )/2, we have C m ≤ (a + 1)dim(S m ) + log (3M ). Therefore, when β is given by ( 14) and d * ≤ κn for some κ < 1, we have

E µ -μ 2 ≤ (1 + ε n ) inf m∈M µ -Π Sm µ 2 + a + 1 β dim(S m )σ 2 + R ′ n σ 2 , with R ′ n = log(3M ) β + 1 2 log n + µ -Π S * µ 2 σ 2 × (a + 1)κ + n -1 log(3M ) β(1 -κ) .
In particular, for a given index of complexity (M, a) and a given κ, the previous bound gives an oracle inequality.

3.2. On the choice of the parameter β. The choice of the tuning parameter β is important in practice. Theorem 1 or Theorem 5 in [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF] justify the choice of a β smaller than 1/4. Nevertheless, Bayesian arguments [START_REF] Hartigan | Bayesian regression using Akaike priors[END_REF] suggest to take a larger value for β, namely β = 1/2. In this section, we discuss this issue on the example of Section 2.2.

For the sake of simplicity, we will restrict to the case where the variance is known. We consider the weights (4) proposed by Leung and Barron, with the probability distribution π given by ( 7)

with α = 1, namely 1 π m = 1 + p -1 -p p -|m| .
According to [START_REF] Bunea | Aggregation for Gaussian regression[END_REF], the estimator μ takes the form

(17) μ = p j=1 s β (Z j /σ)Z j v j with Z j =< Y, v j > and s β (z) = e βz 2 pe 2β + e βz 2 .
To start with, we note that a choice β > 1/2 is not to be recommanded. Indeed, we can compare the shrinkage coefficient s β (Z j /σ) to a threshold at level

T = (2 + β -1 log p)σ 2 since s β (Z j /σ) ≥ 1 2
1 {Z 2 j ≥T } . For µ = 0, the risk of μ is then larger than a quarter of the risk of the threshold estimator μT

= p j=1 1 {Z 2 j ≥T } Z j v j , namely E 0 -μ 2 = p j=1 E s β (Z j /σ) 2 Z 2 j ≥ 1 4 p j=1 E 1 {Z 2 j ≥T } Z 2 j = 1 4 E 0 -μT 2 .
1 Note that this choice of α minimizes the rate of growth of -log πm when p goes to infinity since

-log πm = p log 1 + p -α + α|m| log p = p 1-α + α|m| log p + o p 1-α ,
when πm is given by ( 7) with α > 0. Now, when the threshold T is of order 2K log p with K < 1, the threshold estimator is known to behave poorly for µ = 0, see [START_REF] Birgé | Gaussian model selection[END_REF] Section 7.2. Therefore, a choice β > 1/2 would give poor results at least when µ = 0.

On the other hand, next proposition justifies the use of any β ≤ 1/2 by a risk bound similar to [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF]. For p ≥ 1 and β > 0, we introduce the numerical constants γ β (p) = 2 + β -1 log p and

c β (p) = sup x∈[0,4γ β (p)] R (x -(x + z)s β (x + z)) 2 e -z 2 /2 dz/ √ 2π min(x 2 , γ β (p) 2 ) + γ β (p) 2 /p ∨ 0.6.
This constant c β (p) can be numerically computed. For example, c 1/2 (p) ≤ 1 for any 3 ≤ p ≤ 10 6 , see Figure 1.

Proposition 2. For 3 ≤ p ≤ n and β ∈ [1/4, 1/2], the Euclidean risk of the estimator ( 17) is upper bounded by

(18) E µ -μ 2 ≤ µ -Π S * µ 2 + c β (p) inf m∈M Π S * µ -Π Sm µ 2 + 2 + β -1 log p (|m| + 1)σ 2 .
The constant c β (p) is (crudely) bounded by 16 when p ≥ 3 and

β ∈ [1/4, 1/2].
We delayed the proof of Proposition 2 to Section 6.4. We also emphasize that the bound c β (p) ≤ 16 is crude.

In light of the above result, β = 1/2 seems to be a good choice in this case and corresponds to the choice of Hartigan [START_REF] Hartigan | Bayesian regression using Akaike priors[END_REF]. Note that the choice β = 1/2 has no reason to be a good choice in other situations. Indeed, a different choice of α in (7) would give different "good" values for β. For α ≥ 1, one may check that the "good" values for β are β ≤ α/2.

Remark 3. Proposition 2 does not provide an oracle inequality. The Bound [START_REF] Mallows | Some comments on Cp[END_REF] differs from the best trade off between the bias and the variance term by a log p factor. This is unavoidable from a minimax point of view as noticed in Donoho and Johnstone [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF].

Remark 4. A similar analysis can be done for the estimator (8) when the variance is unknown. When the parameters α and b in ( 8) equal 1, we can justify the use of values of β fulfilling

β ≤ 1 2 φ -1 log p n -p , for n > p ≥ 3,
see the Appendix.

Choice of the models and the weights in different settings

In this section, we propose some choices of weights in three situations: the linear regression, the estimation of functions with bounded variation and regression in Besov spaces.

4.1. Linear regression. We consider the case where the signal µ depends linearly on some observed explanatory variables x (1) , . . . , x (p) , namely

µ i = p j=1 θ j x (j) i , i = 1, . . . , n.
The index i usually corresponds to an index of the experiment or to a time of observation. The number p of variables may be large, but we assume here that p is bounded by n -3.

4.1.1. The case of ordered variables. In some favorable situations, the explanatory variables x (1) , . . . , x (p) are naturally ordered. In this case, we will consider the models spanned by the m first explanatory variables, with m ranging from 0 to p. In this direction, we set S 0 = {0} and S m = span x (1) , . . . , x (m) for m ∈ {1, . . . , p}, where x (j) = (x

(j) 1 , . . . , x (j)
n ) ′ . This collection of models is indexed by M = {0, . . . , p} and contains one model per dimension. Note that S * coincides here with S p . We may use in this case the priors π m = e α -1 e α -e -αp e -αm , m = 0, . . . , p, with α > 0, set L m = m/2 and takes the value ( 14) for β, with N * = n -p. Then, according to Theorem 1 the performance of our procedure is controlled by

E µ -μ 2 ≤ (1 + ε n ) inf m∈M µ -Π Sm µ 2 + σ2 β (α + 1/2)m + 1.2 σ2 β log e α e α -1 + σ 2 2 log n .
with ε n = (2n log n) -1 and σ2 = σ 2 + µ -Π Sp µ 2 /(n -p). As mentioned at the end of Section 3.1, the previous bound can be formulated as an oracle inequality when imposing the condition p ≤ κn, for some κ < 1.

4.1.2. The case of unordered variables. When there is no natural order on the explanatory variables, we have to consider a larger collection of models. Set some q ≤ p which represents the maximal number of explanatory variables we want to take into account. Then, we write M for all the subsets of {1, . . . , p} of size less than q and S m = span x (j) , j ∈ m for any nonempty m. We also set S ∅ = {0} and S * = span x (1) , . . . , x (p) . Note that the cardinality of M is of order p q , so when p is large the value q should remain small in practice.

A possible choice for π m is

π m = p |m| (|m| + 1)H q -1
with

H q = q d=0 1 d + 1 ≤ 1 + log(q + 1).
Again, we choose the value ( 14) for β with N * = n -p and L m = |m|/2. With this choice, combining the inequality p |m| ≤ (e|m|/p) |m| with Theorem 1 gives the following bound on the risk of the procedure

E µ -μ 2 ≤ [1 + ε n ] inf m∈M µ -Π Sm µ 2 + σ2 β |m| 3/2 + log p |m| + log(|m| + 1) + σ 2 2 log n + σ2 β log log[(q + 1)e] ,
with ε n = (2n log n) -1 and σ2 = σ 2 + µ -Π S * µ 2 /(n -p).

Remark. When the family x (1) , . . . , x (p) is orthogonal and q = p, we fall into the setting of Section 2.2. An alternative in this case is to use μ given by [START_REF] Birgé | A generalized Cp criterion for Gaussian model selection[END_REF], which is easy to compute numerically.

4.2. Estimation of BV functions. We consider here the functional setting ( 19)

µ i = f (x i ), i = 1, . . . , n
where f : [0, 1] → R is an unknown function and x 1 , . . . , x n are n deterministic points of [0, 1]. We assume for simplicity that 0 =

x 1 < x 2 < • • • < x n < x n+1 = 1 and n = 2 Jn ≥ 8. We set J * = J n -1 and Λ * = ∪ J * j=0 Λ(j) with Λ(0) = {(0, 0)} and Λ(j) = {j} × 0, . . . , 2 j-1 -1 for j ≥ 1. For (j, k) ∈ Λ * we define v j,k ∈ R n by [v j,k ] i = 2 (j-1)/2 1 I + j,k (i) -1 I - j,k (i) , i = 1, . . . , n with I + j,k = 1 + (2k + 1)2 -j n, .
. . , (2k + 2)2 -j n and I - j,k = 1 + 2k2 -j n, . . . , (2k + 1)2 -j n . The family {v j,k , (j, k) ∈ Λ * } corresponds to the image of the points x 1 , . . . , x n by a Haar basis (see Section 6.5) and it is orthonormal for the scalar product

< x, y > n = 1 n n i=1
x i y i .

We use the collection of models S m = span {v j,k , (j, k) ∈ m} indexed by M = P(Λ * ) and fall into the setting of Section 2.2. We choose the distribution π given by ( 7) with p = n/2 and α = 1. We also set b = 1 and take some β fulfilling

β ≤ 1 2 φ -1 2 log(n/2) n .
According to Proposition 1 the estimator (6) then takes the form

(20) μ = J * j=0 k∈Λ(j)   Z j,k exp nβZ 2 j,k /σ 2 en/2 + exp nβZ 2 j,k /σ 2   v j,k , with Z j,k =< Y, v j,k > n and σ2 = 2   < Y, Y > 2 n - J * j=0 k∈Λ(j) Z 2 j,k   .
Next corollary gives the rate of convergence of this estimator when f has bounded variation, in terms of the norm • n induced by the scalar product < •, • > n .

Corollary 1. In the setting described above, there exists a numerical constant C such that for any function f with bounded variation V (f )

E µ -μ 2 n ≤ C max V (f )σ 2 log n n 2/3 , V (f ) 2 n , σ 2 log n n .
The proof is delayed to Section 6.5. The minimax rate in this setting is (V (f )σ 2 /n) 2/3 . So, the rate of convergence of the estimator differs from the minimax rate by a (log n) 2/3 factor. We can actually obtain a rate-minimax estimator by using a smaller collection of models similar to the one introduced in the next section, but we lose then Formula (20).

4.3.

Regression on Besov space B α p,∞ [0, 1]. We consider again the setting [START_REF] Munk | On difference based variance estimation in nonparametric regression when the covariate is high dimensional[END_REF] with f : [0, 1] → R and introduce a L 2 ([0, 1], dx)-orthonormal family φ j,k , j ≥ 0, k = 1 . . . 2 j of compactly support wavelets with regularity r. We will use models generated by finite subsets of wavelets. If we want that our estimator shares some good adaptive properties on Besov spaces, we shall introduce a family of models induced by the compression algorithm of Birgé and Massart [START_REF] Birgé | Gaussian model selection[END_REF]. This collection turns to be slightly more intricate than the family used in the previous section. We start with some κ < 1 and set J * = ⌊log(κn/2)/ log 2⌋. The largest approximation space we will consider is

F * = span φ j,k , j = 0 . . . J * , k = 1 . . . 2 j ,
whose dimension is bounded by κn. For 1 ≤ J ≤ J * , we define

M J =    m = J * j=0 {j} × A j , with A j ∈ Λ j,J    ,
where Λ j,J = 1, . . . , 2 j when j ≤ J -1 and Λ j,J = A ⊂ 1, . . . , 2 j : |A| = ⌊2 J /(j -J + 1) 3 ⌋ , when J ≤ j ≤ J * .

To m in M = ∪ J * J=1 M J , we associate F m = span {φ j,k , (j, k) ∈ m} and define the model S m by

S m = {(f (x 1 ), . . . , f (x n )), f ∈ F m } ⊂ S * = {(f (x 1 ), . . . , f (x n )), f ∈ F * } .
When m ∈ M J , the dimension of S m is bounded from above by

(21) dim(S m ) ≤ J-1 j=0 2 j + J * j=J 2 J (j -J + 1) 3 ≤ 2 J 1 + J * -J+1 k=1 k -3 ≤ 2.2 • 2 J
and dim(S * ) ≤ κn. Note also that the cardinality of M J is

|M J | = J * j=J 2 j ⌊2 J /(j -J + 1) 3 ⌋ .
To estimate µ, we use the estimator μ given by ( 6) with β given by ( 14) and

L m = 1.1 • 2 J and π m =   2 J (1 -2 J * ) J * j=J 2 j ⌊2 J /(j -J + 1) 3 ⌋   -1 , for m ∈ M J .
Next corollary gives the rate of convergence of the estimator μ when f belongs to some Besov ball B α p,∞ (R) with 1/p < α < r (we refer to De Vore and Lorentz [START_REF] Devore | Constructive approximation. Grundlehren der Mathematischen Wissenschaften[END_REF] for a precise definition of Besov spaces). As it is usual in this setting, we express the result in terms of the norm

• 2 n = • 2 /n on R n .
Corollary 2. For any p, R > 0 and α ∈]1/p, r[, there exists some constant C not depending on n and σ 2 such that the estimator μ defined above fulfills

E µ -μ 2 n ≤ C max σ 2 n 2α/(2α+1) , 1 n 2(α-1/p) , σ 2 n
for any µ given by [START_REF] Munk | On difference based variance estimation in nonparametric regression when the covariate is high dimensional[END_REF]

with f ∈ B α p,∞ (R).
The proof is delayed to Section 6.6. We remind that the rate σ 2 /n 2α/(2α+1) is minimax in this framework, see Yang and Barron [START_REF] Yang | Information-theoretic determination of minimax rates of convergence[END_REF].

Figure 2. Recovering a signal from noisy observations (the crosses). The signal is in black, the estimator is in red.

A numerical illustration

We illustrate the use of our procedure on a numerical simulation. We start from a signal f (x) = 0.7 cos(x) + cos(7x) + 1.5 sin(x) + 0.8 sin(5x) + 0.9 sin(8x), x ∈ [0, 1] which is in black in Figure 2. We have n = 60 noisy observations of this signal

Y i = f (x i ) + σε i , i = 1, . . . , 60
, (the crosses in Figure 2)

where x i = i/60 and ε 1 , . . . , ε 60 are 60 i.i.d. standard Gaussian random variables. The noise level σ is not known (σ equals 1 in Figure 2). To estimate f we will expand the observations on the Fourier basis {1, cos(2πx), . . . , cos(40πx), sin(2πx), . . . , sin(40πx)} .

In this direction, we introduce the p = 41 vectors v 1 , . . . , v 41 given by This vectors {v 1 , . . . , v 41 } form an orthonormal family for the usual scalar product of R n and we fall into the setting of Section 2.2.

v j =            2 n sin(2πjx 1 ), . . . , 2 n sin(2πjx n ) ′ when j ∈ {1, . . . , 20} 1
We estimate (f (x 1 ), . . . , f (x n )) ′ with μ given by ( 8

w m = exp β μm 2 /σ 2 -(α log p + b)|m| m ′ ∈M exp (β μm ′ 2 /σ 2 -(α log p + b)|m ′ |) = exp k∈m βZ 2 k /σ 2 -τ m ′ ∈M exp k∈m ′ βZ 2 k /σ 2 -τ .
We write M j for the set of all the subsets of {1, . . . , j -1, j + 1, . . . , p}. Then, for any j ∈ {1, . . . , p}

c j = m∈M ½ j∈m w m = m∈M ½ j∈m exp k∈m βZ 2 k /σ 2 -τ m∈M exp k∈m βZ 2 k /σ 2 -τ .
Note that any subset m ∈ M with j inside can be written as {j} ∪ m ′ with m ′ ∈ M j , so that

c j = e βZ 2 j /σ 2 -τ m ′ ∈M j exp k∈m ′ βZ 2 k /σ 2 -τ e βZ 2 j /σ 2 -τ m ′ ∈M j exp k∈m ′ βZ 2 k /σ 2 -τ + m∈M j exp k∈m βZ 2 k /σ 2 -τ = e βZ 2 j /σ 2 -τ e βZ 2 j /σ 2 -τ + 1
.

Formula (8) follows.

A preliminary lemma.

Lemma 1. Consider an integer N larger than 2 and a random variable X, such that N X is distributed as a χ 2 of dimension N . Then, for any 0 < a < 1,

(22) E (a -X) + ≤ E a X -1 + ≤ 2 (1 -a)(N -2) exp (-N φ(a)) with φ(a) = 1 2 (a -1 -log a) > 1 4 (1 -a) 2 .
Proof: Remind first that when g : R + → R + and f : R + → R + have opposite monotonicity,

E[g(X)f (X)] ≤ E[g(X)] E[f (X)].
Setting g(x) = x and f (x) = (a/x -1) + leads to the first inequality since E(X) = 1.

We now turn to the second inequality. To start with, we note that

E a X -1 + = aE 1 X - 1 a + = a +∞ 1/a P X ≤ 1 t dt .
Markov inequality gives for any λ ≥ 0

P X ≤ 1 t ≤ e λ/t E e -λX = e λ/t 1 + 2λ N -N/2
and choosing λ = N (t -1)/2 leads to ( 23)

P X ≤ 1 t ≤ 1 t N/2 exp N 2 (1 -1/t) ≤ exp (-N φ(1/t))
for any t > 1. Putting pieces together ensures the bound

E a X -1 + ≤ a +∞ 1/a exp N 2 (1 -1/t) dt t N/2 ≤ a a 0 exp N 2 (1 -x) x N/2-2 dx
for any 0 < a < 1. Iterating integrations by parts leads to

E a X -1 + = exp N 2 (1 -a) 2 a N/2 N -2 k≥0 a k (N/2) k k-1 i=0 (N/2 + i) ≤ exp N 2 (1 -a) 2 a N/2 (N -2)(1 -a) ≤ 2 (1 -a)(N -2) exp (-N φ(a))
for 0 < a < 1, and the bound ( 22) follows. It is also convenient to write w m in the form

w m = π m Z ′ exp -β μ * -μm 2 σ2 -L m with Z ′ = e -β μ * 2 /σ 2 Z.
By construction, the estimator μ belongs to S * , so according to Pythagorean equality we have µ -

μ 2 = µ -µ * 2 + µ * -μ 2 .
The first part is non random, and we only need to control the second part.

According to Theorem 1 in Leung and Barron [START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF], the Stein's unbiased estimate of the L 2 -risk E µ * -μ 2 /σ 2 of the estimator μ on S * can be written as

S (μ) = m∈M w m μ * -μm 2 σ 2 + 2d m -d * - μm -μ 2 σ 2 + 2β∇ μ * -μm 2 σ2
. (μ -μm ) .

When expanding the gradient into

∇ μ * -μm 2 σ2 . (μ -μm ) = 2 (μ * -μm ) σ2 . (μ -μm ) + μ * -μm 2 ∇ 1/σ 2 . (μ -μm ) ,
the term μ * -μm 2 ∇ 1/σ 2 . (μ -μm ) turns to be 0, since ∇ 1/σ 2 is orthogonal to S * . Furthermore, the sum m∈M w m (μ -μ * ) . (μ -μm ) also equals 0, so an unbiased estimate of

E µ * -μ 2 /σ 2 is S (μ) = m∈M w m μ * -μm 2 σ 2 + 2d m + 4β - σ2 σ 2 μm -μ 2 σ2 -d * .
We control the last term thanks to the upper bound

m∈M w m μm -μ 2 = m∈M w m μ * -μm 2 -μ -μ * 2 ≤ m∈M w m μ * -μm 2 and get S (μ) ≤ σ2 σ 2 + 4β - σ2 σ 2 + m∈M w m μ * -μm 2 σ2 + 2 m∈M w m d m -d * ≤ σ2 σ 2 + 4β - σ2 σ 2 + m∈M w m μ * -μm 2 σ2 + L m β -d * + m∈M w m 2d m - σ2 σ 2 + 4β - σ2 σ 2 + L m β where (x) + = max(0, x). First note that when L m ≥ d m /2 we have 2d m - σ2 σ 2 + 4β - σ2 σ 2 + L m β ≤ 2 - σ2 2βσ 2 - 1 2β 4β - σ2 σ 2 + d m ≤ min 0, 2 - σ2 2βσ 2 d m ≤ 0. Therefore, setting δβ = 4βσ 2 /σ 2 -1 + we get S(μ) ≤ σ2 σ 2 1 + δβ m∈M w m μ * -μm 2 σ2 + L m β -d * .
Let us introduce the Kullback divergence between two probability distributions {α m , m ∈ M} and {π m , m ∈ M} on M

D(α|π) = m∈M α m log α m π m ≥ 0
and the function

E π β (α) = m∈M α m μ * -μm 2 σ2 + L m β + 1 β D(α|π).
The latter function is convex on the simplex S + M = α ∈ [0, 1] |M| , m∈M α m = 1 and can be interpreted as a free energy function. Therefore, it is minimal for the Gibbs measure {w m , m ∈ M} and for any α ∈ S + M ,

S(μ) ≤ σ2 σ 2 1 + δβ m∈M α m μ * -μm 2 σ2 + L m β + 1 β D(α|π) - 1 β D(w|π) -d * ≤ 1 + δβ m∈M α m μ * -μm 2 σ 2 + σ2 βσ 2 L m + σ2 βσ 2 D(α|π) -d * .
We fix a probability distribution α ∈ S + M and take the expectation in the last inequality to get

E [S(μ)] ≤ 1 + E( δβ ) m∈M α m E μ * -μm 2 σ 2 + E σ2 βσ 2 (1 + δβ ) m∈M α m L m + D(α|π) -d * .
Since σ2 /σ 2 is stochastically larger than a random variable X with χ 2 (N * )/N * distribution, Lemma 1 ensures that the two expectations

E σ2 σ 2 δβ = E 4β - σ2 σ 2 + and E δβ = E 4βσ 2 σ2 -1 + are bounded by 2 (1 -4β)(N * -2) exp (-N * φ(4β)) , with φ(x) = (x -1 -log(x))/2. Furthermore, the condition N * ≥ 2 + (log n)/φ(4β) enforces 2 (1 -4β)(N * -2) exp (-N * φ(4β)) ≤ 2φ(4β)e -2φ(4β) (1 -4β)n log n ≤ 1 2n log n = ε n .
Putting pieces together, we obtain

E µ -μ 2 σ 2 = µ -µ * 2 σ 2 + E[S(μ)] ≤ µ -µ * 2 σ 2 + (1 + ε n ) m∈M α m E μ * -μm 2 σ 2 + σ2 βσ 2 + ε n m∈M α m L m + D(α|π) -d * ≤ µ -µ * 2 σ 2 + (1 + ε n ) m∈M α m µ * -µ m 2 σ 2 + d * -d m + σ2 βσ 2 (L m + D(α|π)) -d * ≤ (1 + ε n ) m∈M α m µ -µ m 2 σ 2 -d m + σ2 βσ 2 L m + σ2 βσ 2 D(α|π) + ε n d * .
This inequality holds for any non-random probability distribution α ∈ S + M , so it holds in particular for the Gibbs measure

α m = π m Z β exp - β σ2 µ -µ m 2 -d m σ 2 -L m , m ∈ M
where Z β normalizes the sum of the α m s to one. For this choice of α m we obtain

E µ -μ 2 σ 2 ≤ - (1 + ε n )σ 2 βσ 2 log m∈M π m exp - β σ2 µ -µ m 2 -d m σ 2 -L m + ε n d *
which ensures (12) since d * ≤ n. To get [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] simply note that

m∈M π m exp - β σ2 µ -µ m 2 -d m σ 2 -L m ≥ exp - β σ2 µ -µ m * 2 -d m * σ 2 -L m * -log π m *
for any m * ∈ M. The risk of μ is given by

E µ -μ 2 = µ -µ * 2 + p j=1 E (< µ, v j > -s β (Z j /σ)Z j ) 2 .
Note that Z j =< µ, v j > +σ < ε, v j >, with < ε, v j > distributed as a standard Gaussian random variable. As a consequence, when | < µ, v j > | ≤ 4γσ we have

(24) E (< µ, v j > -s β (Z j /σ)Z j ) 2 ≤ c β (p) min(< µ, v j > 2 , γ 2 σ 2 ) + γ 2 σ 2 /p .
If we prove the same inequality for

| < µ, v j > | ≥ 4γσ, then E µ -μ 2 ≤ µ -µ * 2 + c β (p) p j=1 min(< µ, v j > 2 , γ 2 σ 2 ) + γ 2 σ 2 ≤ µ -µ * 2 + c β (p) inf m∈M µ * -µ m 2 + γ 2 (|m| + 1)σ 2 .
This last inequality is exactly the Bound [START_REF] Mallows | Some comments on Cp[END_REF].

To conclude the proof of Proposition 2, we need to check that

E (x -s β (x + Z)(x + Z)) 2 ≤ c β (p)γ 2 for x ≥ 4γ,
where Z is distributed as a standard Gaussian random variable. We first note that

E (x -s β (x + Z)(x + Z)) 2 ≤ 2x 2 E (1 -s β (x + Z)) 2 + 2E(s β (x + Z) 2 Z 2 ) ≤ 2x 2 E (1 -s β (x + Z)) 2 + 2.
We can bound (1 -

s β (x + Z)) 2 as follows (1 -s β (x + Z)) 2 = 1 1 + exp (β [(x + Z) 2 -γ 2 ]) 2 ≤ exp -2β (x + Z) 2 -γ 2 + ≤ exp - 1 2 (x + Z) 2 -γ 2 +
where the last inequality comes from β ≥ 1/4. We then obtain

E (x -s β (x + Z)(x + Z)) 2 ≤ 2 + 4x 2 E exp - 1 2 (x -Z) 2 -γ 2 + 1 Z>0 ≤ 2 + 4x 2 P (0 < Z < x/2) exp - 1 2 (x/2) 2 -γ 2 + P (Z > x/2) ≤ 2 + 2x 2 exp -x 2 /8 + γ 2 /2 + 4x 2 exp -x 2 /8 .
When p ≥ 3 and β ≤ 1/2, we have γ ≥ √ 2 + 2 log 3 and then sup x>4γ

x 2 e -x 2 /8 = 16γ 2 exp(-2γ 2 ).

Therefore,

E (x -s β (x + Z)(x + Z)) 2 ≤ 2 + 16γ 2 2e -3γ 2 /2 + 4e -2γ 2 ≤ 0.6 γ 2 ≤ c β (p)γ 2 ,
where we used again the bound γ 2 ≥ 2 + 2 log 3. The proof of Proposition 2 is complete.

6.5. Proof of Corollary 1. We start by proving some results on the approximation of BV functions with the Haar wavelets.

6.5.1. Approximation of BV functions. We stick to the setting of Section 4.2 with 0 = x 1 < x 2 < • • • < x n < x n+1 = 1 and n = 2 Jn . For 0 ≤ j ≤ J n and p ∈ Λ(j) we define t j,p = x p2 -j n+1 .

We also set φ 0,0 = 1 and φ j,k = 2 (j-1)/2 1 [t j,2k+1 ,t j,2k+2 ) -1 [t j,2k ,t j,2k+1 ) , for 1 ≤ j ≤ J n and k ∈ Λ(j).

This family of Haar wavelets is orthonormal for the positive semi-definite quadratic form

(f, g) n = 1 n n i=1 f (x i )g(x i )
on functions mapping [0, 1] into R. For 0 ≤ J ≤ J n , we write f J for the projection of f onto the linear space spanned by {φ j,k , 0 ≤ j ≤ J, k ∈ Λ(j)} with respect to (•, •) n , namely

f J = J j=0 k∈Λ(j) c j,k φ j,k , with c j,k = (f, φ j,k ) n .
We also consider for 1 ≤ J ≤ J n an approximation of f à la Birgé and Massart [START_REF] Birgé | An adaptive compression algorithm in Besov spaces[END_REF] fJ = f J-1

+ Jn j=J k∈Λ ′ J (j) c j,k φ j,k ,
where Λ ′ J (j) ⊂ Λ(j) is the set of indices k we obtain when we select the K j,J largest coefficients |c j,k | amongs {|c j,k |, k ∈ Λ(j)}, with K j,J = ⌊(j -J + 1) -3 2 J-2 ⌋ for 1 ≤ J ≤ j ≤ J n . Note that the number of coefficients c j,k in fJ is bounded from above by

1 + J-1 j=1 2 j-1 + j≥J (j -J + 1) -3 2 J-2 ≤ 2 J-1 + 2 J-2 p≥1 p -3 ≤ 2 J .
Next proposition states approximation bounds for f J and fJ in terms of the (semi-)norm f 2 n = (f, f ) n . Proposition 3. When f has bounded variation V (f ), we have

(25) f -f J n ≤ 2V (f )2 -J/2 , for J ≥ 0 and (26) f -fJ n ≤ c V (f )2 -J , for J ≥ 1 with c = p≥1 p 3 2 -p/2+1 .
Formulaes ( 25) and ( 26) are based on the following fact.

Lemma 2. When f has bounded variation V (f ), we have

k∈Λ(j) |c j,k | ≤ 2 -(j+1)/2 V (f ), for 1 ≤ j ≤ J n .
Proof of the Lemma. We assume for simplicity that f is non-decreasing. Then, we have

c j,k =< f, φ j,k > n = 2 (j-1)/2 n    i∈I + j,k f (x i ) - i∈I - j,k f (x i )    ,
with I + j,k and I - j,k defined in Section 4.2. Since |I + j,k | = 2 -j n and f is non-decreasing

|c j,k | ≤ 2 (j-1)/2 n |I + j,k | f (x (2k+2)2 -j n ) -f (x (2k)2 -j n ) ≤ 2 -(j+1)/2 f (x (2k+2)2 -j n ) -f (x (2k)2 -j n ) ,
and Lemma 2 follows.

We first prove [START_REF] Yang | Combining different procedures for adaptive regression[END_REF]. Since the {φ j,k , k ∈ λ(j)} have disjoint supports, we have for 0

≤ J ≤ J n f -f J n ≤ j>J k∈Λ(j) c j,k φ j,k n ≤ j>J   k∈Λ(j) |c j,k | 2 φ j,k 2 n =1   1/2 ≤ j>J k∈Λ(j) |c j,k |.
Formula (25) then follows from Lemma 2.

To prove [START_REF] Yang | Mixing Strategies for Density Estimation[END_REF] we introduce the set Λ ′′ J (j) = Λ(j) \ Λ ′ J (j). Then, for 1 ≤ J ≤ J n we have

f -fJ n ≤ Jn j=J   k∈Λ ′′ J (j) |c j,k | 2 φ j,k 2 n =1   1/2 ≤ Jn j=J   max k∈Λ ′′ J (j) |c j,k | k∈Λ(j) |c j,k |   1/2
.

The choice of Λ ′ J (j) enforces the inequalities (1 + K j,J ) max

k∈Λ ′′ J (j) |c j,k | ≤ k∈Λ ′′ J (j) |c j,k | + k∈Λ ′ J (j) |c j,k | ≤ k∈Λ(j) |c j,k |.
To complete the proof of Proposition 3, we combine this bound with Lemma 2:

f -fJ n ≤ j≥J 2 -(j+1)/2 V (f )(1 + K j,J ) -1/2 ≤ j≥J 2 -(j+1)/2 V (f )2 -(J-2)/2 (j -J + 1) 3 . ≤ V (f )2 -J p≥1 p 3 2 -p/2+1 .
6.5.2. Proof of Corollary 1. First, note that v j,k = (φ j,k (x 1 ), . . . , φ j,k (x n )) ′ for (j, k) ∈ Λ * . Then, according to ( 25) and ( 26) there exists for any 0

≤ J ≤ J * a model m ∈ M fufilling |m| ≤ 2 J and µ -Π Sm µ 2 n = µ -Π S * µ 2 n + Π S * µ -Π Sm µ 2 n ≤ 2c 2 V (f ) 2 2 -J * ∨ 2 -2J ,
with c = p≥1 p 3 2 -p/2+1 . Putting together this approximation result with Theorem 1 gives

E µ -μ 2 n ≤ C inf 0≤J≤J * V (f ) 2 2 -J * ∨ 2 -2J + 2 J log n n σ 2
for some numerical constant C, when

β ≤ 1 4 φ -1 log n n/2 -2 .
This bound still holds true when

1 4 φ -1 log n n/2 -2 ≤ β ≤ 1 2 φ -1 log(n/2) n/2 ,
see Proposition 4 in the Appendix. To conclude the proof of Corollary 1, we apply the previous bound with J given by the minimum between J * and the smallest integer such that

2 J ≥ V (f )
m ∈ M J -log π m ≤ log 2 J + J * j=J 2 J (j -J + 1) 3 log e2 j-J+1 (j -J + 1) 3 ≤ 2 J   1 + k≥1 k -3 (1 + 3 log k + k log 2)   ≤ 4 • 2 J . ( 27
)
Second, when f belongs to some Besov ball B α p,∞ (R) with 1/p < α < r, Birgé and Massart [START_REF] Birgé | An adaptive compression algorithm in Besov spaces[END_REF] gives the following approximation results. There exists a constant C > 0, such that for any

J ≤ J * and f ∈ B α p,∞ (R), there exists m ∈ M J fulfilling f -ΠFm f ∞ ≤ f -ΠF * f ∞ + ΠF * f -ΠFm f ∞ ≤ C max 2 -J * (α-1/p) , 2 -αJ
where ΠF denotes the orthogonal projector onto F in L 2 ([0, 1]). In particular, under the previous assumptions we have

µ -Π Sm µ 2 n ≤ 1 n n i=1 f (x i ) -ΠFm f (x i ) 2 ≤ f -ΠFm f 2 ∞ ≤ C 2 max 2 -2αJ , 2 -2J * (α-1/p) , (28)
To conclude the proof of Corollary 2, we combine Theorem 1 together with ( 21), ( 27) and ( 28) for Assume also that b is not smaller than 1. Then, we have the following upper bound on the L 2 -risk of the estimator μ defined by (31)

J = min J * , log[max(n/σ 2 ,
(33) E µ -μ 2 ≤ 16 inf m∈M µ -µ m 2 + 1 β (b + log p) (|m| + 1)σ 2 + (2 + b + log p)σ 2 ,
where σ2 = σ 2 + µ -Π S * µ 2 /(n -p).

We postpone the proof of Proposition . So all we need is to prove [START_REF] Wang | Effect of mean on variance function estimation in nonparametric regression[END_REF]. For symmetry reason, we restrict to the case a > 0.

To prove Inequality [START_REF] Wang | Effect of mean on variance function estimation in nonparametric regression[END_REF], we first note that

E   a - exp β(a + ε) 2 p exp(βλ) + exp (β(a + ε) 2 ) (a + ε) 2   = E a 1 + p -1 exp (β [(a + ε) 2 -λ]) - ε 1 + p exp(β [λ -(a + ε) 2 ]) 2 ≤ 2a 2 E 1 1 + p -1 exp (β [(a + ε) 2 -λ]) 2 + 2E ε 1 + p exp(β [λ -(a + ε) 2 ]) 2 . ( 34
)
and then investigate apart the four cases 0 ≤ a ≤ 2γ -1 , 2γ -1 ≤ a ≤ 1, 1 ≤ a ≤ √ 3 γ and a ≥ √ 3 γ. Inequality (29) will follow from Inequalities (36), (37), (38) and (39).

Case 0 ≤ a ≤ 2γ -1 . From (34) we get

E   a - exp β(a + ε) 2 p exp(βλ) + exp (β(a + ε) 2 ) (a + ε) 2   ≤ 2a 2 + 2E ε 1 + p exp(β [λ -(a + ε) 2 ]) 2 ≤ 2a 2 + 4E ε1 ε>0 1 + p exp(β [λ -(a + ε) 2 ]) 2 ≤ 2a 2 + 4E ε 2 1 ε>0 exp -2β γ 2 -(a + ε) 2 + with γ 2 = λ + β -1 log p. Expanding the expectation gives E ε 2 1 ε>0 exp -2β γ 2 -(a + ε) 2 + = e -2βγ 2 γ-a 0 x 2 e 2β(a+x) 2 -x 2 /2 dx √ 2π + ∞ γ-a
x 2 e -x 2 /2 dx √ 2π

Note that when p ≥ 3, λ ≥ 2 and β ≤ 1/2, we have γ ≥ 2. Therefore, when 0 ≤ a ≤ 1, an integration by parts in the second integral gives

∞ γ-a x 2 e -x 2 /2 dx √ 2π = -xe -x 2 /2 √ 2π ∞ γ-a + ∞ γ-a e -x 2 /2 dx √ 2π ≤ 2(γ -a) √ 2π e -(γ-a) 2 /2 .
For the first integral, since β > 1/4, we have the bound

e -2βγ 2 γ-a 0 x 2 e 2β(a+x) 2 -x 2 /2 dx √ 2π ≤ e -(γ-a) 2 /2 γ-a 0 x 2 dx √ 2π ≤ (γ -a) 3 3 √ 2π e -(γ-a) 2 /2 .
Besides an integration by parts also gives

e -2βγ 2 γ-a 0 x 2 e 2β(a+x) 2 -x 2 /2 dx √ 2π ≤ (4β -1) -1 e -2βγ 2 γ-a 0 (4β(a + x) -x)xe 2β(a+x) 2 -x 2 /2 dx √ 2π ≤ e -2βγ 2 (4β -1) √ 2π xe 2β(a+x) 2 -x 2 /2 γ-a 0 = γ -a (4β -1) √ 2π e -(γ-a) 2 /2 .
Putting pieces together, we obtain for 0 ≤ a ≤ 1 and 1/4 The random variable σ2 /σ 2 is stochastically larger than a random variable X distributed as a χ 2 of dimension n -p divided by n -p. The Lemma 1 gives

< β ≤ 1/2 E   a - exp β(a + ε) 2 p exp(βλ) + exp (β(a + ε) 2 ) (a + ε)
P (X ≤ 2β) ≤ (2β) N/2 exp N 2 (1 -2β) ≤ exp (-N φ(2β))
and then P σ2 < 2βσ 2 ≤ exp(-(n -p)φ(2β)), so Condition (32) ensures that p P σ2 < 2βσ 2 ≤ 1.

Finally, since 12(b + log p) + 36 is smaller than 16(2 + b + log p) we get (33).

Université de Nice Sophia-Antipolis, Laboratoire J-A Dieudonné, Parc Valrose, 06108 Nice cedex 02 E-mail address: cgiraud@math.unice.fr

Remark 1 .

 1 Other choices are possible for {π m , m ∈ M} and they lead to different c j s. Let us mention the choice π m = (p + 1) p |m| -1

Theorem 1 .

 1 Assume that β and N * fulfill the condition (11) β < 1/4 and N * ≥ 2 + log n φ(4β) ,

Figure 1 .

 1 Figure 1. Plot of p → c 1/2 (p)

2 n

 2 cos(2π(j -21)x 1 ), . . . , 2 n cos(2π(j -21)x n ) ′ when j ∈ {22, . . . , 41} .

2 . 6 . 6 . 1 .

 2661 ) with the parameter α = 1, b = 1 and β = 1/3. Finally, we estimate f withf n) < μ, v 21 > and âj = (2/n) < μ, v j+21 >, bj = (2/n) < μ, v j > for j = 1, . . . ,20. The plot of f is in red in Figure Proofs Proof of Proposition 1. Let us first express the weights w m in terms of the Z j s and τ = α log p + b. We use below the convention, that the sum of zero term is equal to zero. Then for any m ∈ M, we have

6. 3 .

 3 Proof of Theorem 1. To keep formulas short, we write d m for the dimension of S m , and use the following notations for the various projections μ * = Π S * Y, µ * = Π S * µ and µ m = Π Sm µ, m ∈ M.

6. 4 .

 4 Proof of Proposition 2. We use along the proof the notations µ * = Π S * µ, [x] + = max(x, 0) and γ = γ β (p) = 2 + β -1 log p. We omit the proof of the bound c β (p) ≤ 16 for β ∈ [1/4, 1/2] and p ≥ 3. This proof (of minor interest) can be found in the Appendix.

2 ≤ 2 ≤ 2 ≤ 2  2  ≤ 2a 2 + 5 . 6 γ 2 ≤ 2a 2 + 5 . 6 γ

 222222562256 min (4β -1) -1 , 3 -1 (γ -a) 2 √ 2π (γ -a)e -(γ-a) 2 /2a 2 + 8 + 4 min (4β -1) -1 , 3 -1 (γ -a) 2 e -(1/2-β)(γ-a) 2 √ 2π (γ -a)e -β(γ-a) 2a 2 + 8 + 4 min (4β -1) -1 , [3e(1/2 -β)] -1 √ 2π γe -β(γ-a) 2a 2 + 5.6 γ e -β(γ-a) 2 . (35)Furthermore, when 0 ≤ a ≤ 2γ -1 , we have (γ -a) 2 ≥ γ 2 -4. The inequalities (35) and λ ≥ 2 thus giveE   a -exp β(a + ε) 2 p exp(βλ) + exp (β(a + ε) 2 ) (a + ε) ≤ 2a 2 + 5.6e (4-λ)β γ/p ≤ 2a 2 + 8 γ 2 /p . (36) Case 2γ -1 ≤ a ≤ 1. Starting from (35) we have for 2γ -1 ≤ a ≤ 1 E   a -exp β(a + ε) 2 p exp(βλ) + exp (β(a + ε) 2 )(a + ε) 3 e -β(γ-1) 2 γ 3 e -(γ-1)

  Proof of Corollary 2. First, according to the inequality n k ≤ (en/k) k we have the bound for

	2/3	n σ 2 log n	1/3	.
	6.6.

  Proposition 4. Assume that β and p fulfill the conditions,

	(32)	p ≥ 3, 0 < β < 1/2 and p +	log p φ(2β)	≤ n.
		1)] (2α + 1) log 2	+ 1 .

  4 to Section B.2. Proof of Lemma 3. When β ≤ 1/4 the bound (30) follows from a slight variation of Theorem 5 in[START_REF] Leung | Information Theory and Mixing Least-Squares Regressions[END_REF]. When 1/4 < β ≤ 1/2, it follows from (29) (after a rescaling by σ 2 and a summation)

	Appendix B. Proofs
	B.1.

  2 /4 γ 2 ≤ 2a 2 + 56γ -2 ≤ 16a 2 . -1 exp (β [(a -ε) 2 -λ]) ≤ 4a 2 E exp -2β (a -ε) 2 -γ 2 + 1 ε>0 + 2 ≤ 2a 2 exp -2β (2a/3) 2 -γ 2 + + 4a 2 P (ε > a/3) + 2.Since (2a/3) 2 > a 2 /9 + γ 2 , we finally obtain Proof of Proposition 4. We can express the weights c j appearing in (31) in the following wayc j = e βZ 2 j /σ 2 pe βλ + e βZ 2 j /σ 2 with β = βσ 2 /σ 2 and λ = b/ β. Note that β ≤ 1/2 enforces λ ≥ 2 since b ≥ 1.Since σ2 is independent of the Z j s, we can work conditionally on σ2 . When β is not larger than 1/2 we apply Lemma 3 and getE µ -μ 2 | σ2 1 { β≤1/2} ≤ 16 infWhen β is larger than 1/2, we use the following bound.Lemma 4. Write Z for a Gaussian random variable with mean a and variance σ 2 . Then for any λ ≥ 0 and β > 1/2 pe βλ + e βZ 2 /σ 2 Z e βZ 2 /σ 2 pe βλ + e βZ 2 /σ 2 Z + σ 2 ). Whenever a 2 is smaller than 3 λ + β-1 log p σ 2 , this quantity remains smaller than 6 λ + β-1 log p σ 2 + 2σ 2 . When a 2 is larger than 3 λ + β-1 log p σ 2 we follow the same lines as in the last case of Section B.1 and getpe β λ + e βZ 2 /σ 2 Z Since (2a/3) 2 ≥ a 2 /9 + λ + β-1 log p σ 2 , we obtain E   a -e βZ 2 /σ 2 pe βλ + e βZ 2 /σ 2 Z ≤ 2σ 2 + 4a 2 exp -a 2 /(18σ 2 ) + 2a 2 exp -a 2 /(9σ 2 ) ≤ 36σ 2 . E µ -μ 2 | σ2 1 { β>1/2} ≤ p 6 β (b + log p) σ2 + 36σ 2 1 { β>1/2} + µ -Π S * µ 2 1 { β>1/2} .Futhermore σ2 is smaller than 2βσ 2 when β is larger than 1/2, so taking the expectation of (42) and (40) givesE µ -μ 2 ≤ 16 inf

	(37) B.2. m∈M From (41), we obtain after summation ≤ 16 inf m∈M (40) (42) m∈M	µ -µ m µ -µ m µ -µ m	2 + λ + β-1 log p (|m| + 1)σ 2 1 { β≤1/2} 2 + 1 β (b + log p) (|m| + 1)σ 2 1 { β≤1/2} . 2 + 1 β (b + log p) (|m| + 1)σ 2
	Case 1 ≤ a ≤ E   a -√ 3 γ. From (34) we get p exp(βλ) + exp (β(a + ε) 2 ) exp β(a + ε) 2 (41) E   a -e βZ 2 /σ 2	(a + ε)
	(38)		2	
	Case a >	√	3γ. From (34) we get	
	E 1 + p -1 exp β (a/σ + ε) 2 -λ   a -exp β(a + ε) 2 p exp(βλ) + exp (β(a + ε) 2 ) ≤ 4a 2 E 1 1 + p 2 (a + ε) 2 E   a -exp β(a + ε) 2 p exp(βλ) + exp (β(a + ε) 2 ) (a + ε)   1 ε>0 + 2 2   a -σε 1 + p exp β λ -(a/σ + ε) 2 2  e βZ 2 /σ 2   ≤ 2(a 2 E  = E    a -	 	2  
	(39)		≤ 2a 2 exp -2βa 2 /9 + 4a 2 exp -a 2 /18 + 2 ≤ 42 ≤ 11γ 2 . ≤ 2σ 2 + 4a 2 P (ε > |a|/(3σ)) + 2a 2 p 2 exp -2 β	2a 3σ	2	-λ	.

2   ≤ 2(a 2 + 1) ≤ 4a 2 ≤ 12 min(a 2 , γ 2 ). 2   ≤ 6 λ + β-1 log p σ 2 + 36σ 2 .

Proof. First, we write ε for a standard Gaussian random variable and obtain

E   a 2   + p [12 (b + log p) + 36] P σ2 < 2βσ 2 σ 2 .

APPENDIX

The Appendix is devoted to the proof of the bound c β (p) ≤ 16 and gives further explanations on the Remark 4, Section 3.2. The results are stated in Appendix A and the proofs (of minor interest) can be found in Appendix B.

Appendix A. Two bounds

A.1. When the variance is known. In this section, we assume that the noise level σ is known. We remind that 3 ≤ p ≤ n and {v 1 , . . . , v p } is an orthonormal family of vectors in R n . Next lemma gives an upper bound on the Euclidean risk of the estimator

Lemma 3. For any λ ≥ 2 and β > 0, we set γ 2 = λ + β -1 log p.

1. For 1/4 ≤ β ≤ 1/2 and a ∈ R we have the bound

where ε is distributed as a standard Gaussian random variable. 2. As a consequence, for any 0 < β ≤ 1/2 and λ ≥ 2 the risk of the estimator μ is upper bounded by

Note that (29) enforces the bound c β (p) ≤ 16. This constant 16 is certainly far from optimal. Indeed, when β = 1/2, λ = 2 and 3 ≤ p ≤ 10 6 , Figure 1 in Section 3.2 shows that the bounds ( 29) and (30) hold with 16 replaced by 1.

A.2. When the variance is unknown. We consider the same framework, except that the noise level σ is not known. Next proposition provides a risk bound for the estimator We remind that φ(x) = (x -1 -log x)/2 .