Didier A Girard

ESTIMATING THE ACCURACY OF (LOCAL) CROSS-VALIDATION VIA RANDOMIZED GCV CHOICES IN KERNEL OR SMOOTHING SPLINE REGRESSION

Keywords: Bootstrap, C L criterion, Coverage probability, Generalized cross-validation, Kernel regression, Nonparametric regression, Randomized trace, Regularization, Smoothing spline

In nonparametric regression, it is generally crucial to select "nearly" optimal smoothing parameters for which a given (weighted) average squared error (∆) is "nearly" minimized. The cross-validation (CV) selector, or the GCV selector, are popular for this task but it has been observed by many statisticians that these selectors may happen to be "not sufficiently" accurate in some situations. So a practical matter of great importance is the development of reliable estimates of this accuracy.

The purpose of this paper is to show that the simulation of the randomized GCV selector or a simple general variant using an "augmented-randomized-trace", can provide useful inferences, like consistent estimates of the standard error in the CV selector or of the expected increase of ∆ due to this error. Furthermore this also provides a tool for constructing more parsimonious curve estimates having almost the same asymptotic justification as the CV estimate, namely with similar increase of ∆ up to a given factor.

Rigorous proofs are given in the context of one-dimensional kernel regression. Simulated examples, also in this context, illustrate the usefulness of the methodology even at moderate sample sizes. Some direct extensions (for multi-dimensional kernels, equispaced splines) of the theoretical results are outlined. We give heuristics which indicate that the general methodology proposed in this article should be useful in many curve-, surface-or image-estimation problems when using spline-like smoothers.

Introduction

When regularization methods or nonparametric techniques are used to solve an inverse problem or to estimate the mean curve or surface underlying noisy observations, the choice of the regularization parameter(s), smoothing parameter(s) or bandwidth(s) is generally crucial. Cross-validation (CV), unbiased risk estimate (C L ) and generalized cross-validation (GCV) methods are very popular for this choice. They were studied and applied in a broad variety of contexts: see [START_REF] Craven | Smoothing noisy data with spline functions[END_REF], [START_REF] Rice | Bandwidth choice for nonparametric regression[END_REF], [START_REF] Speckman | Spline smoothing and optimal rates of convergence in nonparametric regression models[END_REF], [START_REF] Li | From Stein's unbiased risk estimates to the method of generalized cross-validation[END_REF][START_REF] Li | Asymptotic optimality of C L and generalized cross-validation in ridge regression with application to spline smoothing[END_REF], [START_REF] Wahba | A comparison of GCV and GML for choosing the smoothing parameter in the generalized spline smoothing problem[END_REF], Härdle and Marron (1985), Härdle, [START_REF] Ärdle | How far are automatically chosen regression smoothing parameters from their optimum (with discussion)?[END_REF], [START_REF] Kneip | Ordered linear smoothers[END_REF], [START_REF] Girard | Asymptotic comparison of (partial) cross-validation, GCV and randomized GCV in nonparametric regression[END_REF], for theoretical results, [START_REF] Hastie | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF], [START_REF] Kohn | The performance of cross-validation and maximum likelihood estimators of spline smoothing parameters[END_REF], [START_REF] Thompson | A cautionary note about crossvalidatory choice[END_REF][START_REF] Thompson | A study of methods of choosing the smoothing parameter in image restoration by regularization[END_REF] for extensive experimentations, [START_REF] Eubank | Nonparametric regression and spline smoothing[END_REF] and [START_REF] Gu | Smoothing spline ANOVA models[END_REF] for recent surveys. However, it is now well known that the resulting best bandwidth estimates may exhibit too large a sample-to-sample variability in certain applications, and thus the computed cross-validated curve or surface (from the single observed sample) should be interpreted with great caution in such a case. Various other bandwidth selectors have been recently proposed (especially for nonparametric techniques of the kernel type) which can give more stable estimates, at the price of some bias. However, their asymptotic justification needs stronger smoothness assumptions and in practice their use is not so 'automatic' (additional parameters). As mentioned by several authors (e.g. [START_REF] Hart | Comment on "Empirical functionals and efficient smoothing parameter selection[END_REF] in the discussion of [START_REF] Hall | Empirical functionals and efficient smoothing parameter selection (with discussion)[END_REF], [START_REF] Loader | Bandwidth selection : Classical or plug-in ?[END_REF]), cross-validation is more "assumptions-robust", and is also generally recognized as an unbiased selector: this makes cross-validation an irreplaceable tool. What is presently missing to any practitioner who uses GCV is a reasonable estimate (or bound) of how far from the optimal smoothing parameter is the GCV choice or at least some feeling about its usefulness, as was noticed by [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF] and others. Confidence interval would be even more appreciated but is a more ambitious goal. Such accuracy estimates would turn out to be very useful to guide the practitioner in his/her choices among curve-estimators of various degrees of sophistication. One instance is the choice between global and local bandwidths: Schukany (1995) aptly points out that the development of an indicator stating whether one has sufficient data to warrant more than a global bandwidth, would be very useful.

In this paper we propose a simulation-based methodology for producing estimates of the accuracy of the CV or GCV selector either in the bandwidth-space or in the curve-or surface-space. These estimators will be obtained as a simple by-product of the repeated use of the fast randomized version of GCV or of variants, where "fast" stands for meaning that one uses a single randomized-trace in place of the exact trace.

In [START_REF] Girard | On the fast Monte-Carlo cross-validation and C L procedures: comments, new results and application to image recovery problems (with discussion)[END_REF], computing several replications of the minimizer of the fast randomized GCV criterion was already advocated, even in contexts where exact GCV is computationally feasible, as a simple methodology for providing a useful indicator of certain "pathological" contexts: if the observable variability in the final solution (curve or surface) due to randomization, is "too large" for the problem at hand then (even exact) GCV should not be considered as "reliable" for this problem (see Sections 3.2 and 5.4 of [START_REF] Girard | On the fast Monte-Carlo cross-validation and C L procedures: comments, new results and application to image recovery problems (with discussion)[END_REF]).

The aim of this paper is to show that the distribution of the randomized GCV selector and a variant using an "augmented-randomized-trace", obtained by such a simulation (such a distribution will be called a randomization-based distribution) can actually provide considerably more precise inferences in contexts where GCV is relevant and where enough regularity is present. Before describing this new methodology, let us recall some definitions and previous results.

Simplified background and general definitions

Let us first consider the simple observational model

y i = m(x i ) + ε i , i = 1, • • • , n (1.1)
where m is an unknown smooth function observed at the points x i = i/n, i = 1, • • • , n and ε i are independent and identically distributed observation errors with mean zero and variance σ 2 . Consider the widely studied kernel regression estimator for m of the form:

mh (x) = 1 nh n j=1 K( x -x j h )y j . (1.2)
where K is a smooth "bell shaped" symmetric function chosen so that it satisfies K(x) dx = 1 and

x 2 K(x) dx > 0, and h > 0 is the smoothing parameter (or bandwidth) to be chosen.

This curve estimator evaluated at the points x i will also be matricially denoted by: mh = A h y (1.3)

where the so-called smoother matrix A h is explicit here [A h ] i,j = 1 nh K( xi-xj h ). This will permit us to concisely formulate some future general definitions and discussions which will be applicable to other common smoothers A h like smoothing splines or "lowess" estimates (local weighted least squares).

A convenient measure for assessing mh as an estimate of m, is the weighted average squared error

∆(h) := n -1 [ mh (x i ) -m(x i )] 2 u(x i ) = n -1 A h y -m 2 U (1.4)
(that we will sometimes simply call the "loss") or its expectation M (h):=E(∆(h)). Here u is a fixed nonnegative function (we denote by U the diagonal matrix diag(u

(x i ), i = 1, • • • , n
) and • U is the weighted l 2 norm associated with the inner product x, y U = x T U y). By assuming u compactly supported on a subinterval of (0,1), the boundary effects are eliminated. By using a "localized" weight function u, like the characteristic function of a "small" subinterval, this also classically permits one to target a "local bandwidth". Let us denote by ĥ0 the minimizer of ∆ and h 0 the minimizer of M . The "optimal bandwidth" will refer to ĥ0 in this paper. (For simplification of notation, the dependence of A h , ∆, M , h 0 , ĥ0 , ĥCV , etc., on n are suppressed.) The popular "leave-one-out" approach (or ordinary cross-validation) for estimating the optimal bandwidth, consists, in our linear framework (1.3)-(1.4), of numerically minimizing the criterion

CV(h) = n -1 D -1 h (I -A h )y 2 U D h := diag(1 -[A h ] i,i , i = 1, • • • , n) (1.5)
which, in the particular setting here, also coincides with the popular generalized cross-validation criterion (e.g. [START_REF] Craven | Smoothing noisy data with spline functions[END_REF], [START_REF] Hastie | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF] Section 3.4.3). Let ĥCV the resulting bandwidth. It is now well known that, for the context (1.1)-(1.2), under standard regularity conditions, we have h 0 ∼ c(K, m, σ, u)n -1/5 as n → ∞ and both ĥ0 /h 0 , ĥCV /h 0 tend to 1 in probability, and moreover the "errors" in ĥCV or h 0 have asymptotic normal distributions:

L n 3/10 (h 0 -ĥ0 ) → N (0, σ 2 1 ), L n 3/10 ( ĥCV -ĥ0 ) → N (0, σ 2 2 ),
where σ 1 , σ 2 are constants which are known functions of K, m, σ, u; see [START_REF] Rice | Bandwidth choice for nonparametric regression[END_REF], Härdle, [START_REF] Ärdle | How far are automatically chosen regression smoothing parameters from their optimum (with discussion)?[END_REF] abbreviated by HHM in the following. In the latter paper, Remark 3.6, it is suggested that, by estimating the unknown terms in σ 2 , this asymptotic result could be used to provide approximate confidence intervals for ĥ0 .

Before going on, we have to recall some other definitions. Let us define

t(h) := trU A h /trU, (1.6)
a trace-term which is simply n -1 h -1 K(0) for the particular estimator (1.2) in the case of an equidistant design, and is thus independent of u in this case. Generalized cross-validation (GCV) is a member of a family of criteria that can be written as

G X (h) := n -1 (I -A h )y 2 U Ξ X (t(h)) (1.7)
where Ξ X is a correction factor (or penalization) satisfying Ξ X (t) = 1 + 2t + O(t 2 ) with Ξ ′′ X bounded on a neighborhood of 0. A list of usual penalizations Ξ X is presented in HHM. GCV which is defined by Ξ GCV (t) := (1t) -2 , is one of the most popular. It has been shown in HHM that, for the setting (1.1)-(1.4), each of these penalized criteria gives a bandwidth equivalent up to second order to ĥCV (that is, with the same normal limiting distribution as the one above). Let ĥ generically denote any one of these GCV-type selectors. Now, the fast randomized version R G of G is obtained by using in place of t(h), in the definition (1.7), the general randomized trace function :

R t(h) := w, A h w U w, w U ,
where w is a simulated unitary "white noise" vector w of size n, i.e. such that Ew = 0 and Var(w) = I (the same w being used for every h). These criteria have been introduced in Girard (1989) (in the unweighted case) as fast Monte-Carlo-type approximations to the exact ones, for all the contexts where computing trU A h is not an easy task: this is not the case in the setting (1.1)-(1.4) above, but typical important examples are smoothing splines or penalized least squares procedures, additive modeling by backfitting, iterative image restorations, etc.; see [START_REF] Girard | On the fast Monte-Carlo cross-validation and C L procedures: comments, new results and application to image recovery problems (with discussion)[END_REF] for references to various applications.

In the setting above, it has been shown in [START_REF] Girard | Asymptotic comparison of (partial) cross-validation, GCV and randomized GCV in nonparametric regression[END_REF] that

L n 3/10 ( ĥR -ĥ0 ) → N (0, σ 2 R ), with σ 2 R < 2σ 2 2
, where ĥR denotes a generic minimizer of any one of these randomized GCV-type criteria.

Outline of results and structure of the article

In this paper, we propose a rather general method for building up inferences concerning the optimally smoothed curve, optimal in terms of the loss ∆(•). It simply consists of repeatedly rerunning, for example, the randomized GCV procedure, and, in order to approximate the desired distribution of ĥGCV -ĥ0 , of using the empirical distribution of the so-obtained randomized choices ĥRGCV which is thus conditional to the data y, centered about its mean or about ĥGCV (these two centerings being asymptotically equivalent).

We claim (see Section 2, Corollary 2.3, and Remark 8.2 for extensions of the theory) that this permits to construct, under classical regularity conditions, a consistent estimate (in probability) of the distribution of ĥGCV -ĥ0 as n → ∞. To be more precise, the obtained conditional distribution slightly underestimates the targeted distribution in the sense L ĥRGCV -ĥGCV y and L κ( ĥGCV -ĥ0 ) have the same limit in probability with respect to y, where the constant κ < 1 is function of K and in general of m and u (and of the density of the x i in the nonequispaced case). We used the word "slightly" because κ will be shown to be often "reasonably close" to 1 (Section 3) : for example, for second order kernel, if u is a characteristic function, then 1/ √ 3 ≤ κ ≤ 1. However, in some cases (e.g. kernel of high order) κ may be too small.

We then show in Sections 4 and 5 that this possible trouble can be completely eliminated via a second simulation method similar as the above one, except that the randomized trace is appropriately modified so that it has an "augmented randomness". Note, that the first simulation study, with second order kernel, that we will describe in Section 8, demonstrates that even without such a correction, these randomizationbased estimates of var ĥGCV -ĥ0 are satisfactory. After some comments on how to use confidence bounds on ĥ0 , we also describe in Section 7 two other useful byproducts of a simulated population of randomized choices. An additional attractive property is that the variation of the GCV criterion over such a population yields a simple consistent estimate of the supplement of risk E ∆( ĥGCV ) -E ∆( ĥ0 ) from estimation of ĥ0 (or "risk regret"). An appealing feature of the methodology is that its theoretical justification does not require any additional smoothness assumption on m.

General heuristics and other definitions

The underlying idea can be more easily described by considering instead the Mallows unbiased risk estimate which is a basic criterion for choosing h when σ 2 is known: [START_REF] Mallows | Some comments on C P[END_REF]). The fast randomized version of CL is

CL(h) := n -1 (I -A h )y 2 U + 2σ 2 n -1 trU A h (
RCL(h) := n -1 (I -A h )y 2 U + 2σ 2 n -1 w, A h w U
The following discussion will remain relevant for cross-validation or GCV-type criteria, because it has been shown that typically CL (respectively RCL) produces a smoothing parameter equivalent up to second order (and essentially equivalent in practice; see e.g. [START_REF] Hastie | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF] Section 3, Girard 1995) to any one of the G-selectors (1.6)-(1.7) (respectively R G-selectors). Now algebraic manipulations show that:

CL(h) -n -1 ε ε ε ε ε ε ε ε ε ε 2 U -∆(h) = -2n -1 ε ε ε ε ε ε ε ε ε ε, A h ε ε ε ε ε ε ε ε ε ε U -σ 2 trU A h + 2n -1 ε ε ε ε ε ε ε ε ε ε, (I -A h )m U = e 1 (ε ε ε ε ε ε ε ε ε ε, h) + e 2 (m, ε ε ε ε ε ε ε ε ε ε, h),
(1.8) say, while the analog intrinsic error of the randomized version is

RCL(h) -n -1 ε ε ε ε ε ε ε ε ε ε 2 U -∆(h) = e 1 (ε ε ε ε ε ε ε ε ε ε, h) + e 2 (m, ε ε ε ε ε ε ε ε ε ε, h) -e 1 (ε ε ε ε ε ε ε ε ε ε * , h) (1.9) where ε ε ε ε ε ε ε ε ε ε * = σw is independent and distributed identically to ε ε ε ε ε ε ε ε ε
ε at least for the first and second moments.

This means that the additional "randomization error" -e 1 (ε

ε ε ε ε ε ε ε ε ε * , h
) is similar to one of the two components of the intrinsic error of CL (see [START_REF] Girard | On the fast Monte-Carlo cross-validation and C L procedures: comments, new results and application to image recovery problems (with discussion)[END_REF] for an extension to the randomized GCV criterion). Now, for the case A h a kernel smoother, it is well known that, under classical regularity conditions, the asymptotic distribution of ĥ is typically obtained by the following linearization of the equation G ′ ( ĥ) = 0, where " ′ " denotes differentiation w.r.t. h:

0 = G ′ ( ĥ) ≈ G ′ (h 0 ) + ( ĥ -h 0 )M ′′ (h 0 )
where the error in "≈" is small enough relatively to the centered normal limit distribution of G ′ (h 0 ); and it is known that this also typically holds (with linearization errors of the same order) for the couple G ′ , ĥ replaced by CL ′ , ĥ or ∆ ′ , ĥ0 [START_REF] Rice | Bandwidth choice for nonparametric regression[END_REF], HHM), and also by RCL ′ , ĥR as may be expected from the above comments [START_REF] Girard | Asymptotic comparison of (partial) cross-validation, GCV and randomized GCV in nonparametric regression[END_REF]. As discussed in HHM, [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF], Girard (1998, Section3) and others, one can anticipate that the so-obtained three linearizations can yield satisfying approximations also for numerous other classes of linear smoothers A h , provided A h is twice continuously differentiable w.r.t. h.

Now by linearly combining these three approximations using the relations (1.8) and (1.9), we obtain

-M ′′ (h 0 )( ĥ -ĥ0 ) ≈ e ′ 1 (ε ε ε ε ε ε ε ε ε ε, h 0 ) + e ′ 2 (m, ε ε ε ε ε ε ε ε ε ε, h 0 ), (1.10) -M ′′ (h 0 )( ĥR -ĥ) ≈ -e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 ). (1.11)
It is important to note that the right hand term of this last approximation is independent of ε

ε ε ε ε ε ε ε ε ε.
This will permit us, under usual regularity conditions, to show the following: the observable fluctuations of ĥR for a given y, about the associated ĥ, tend (in probability) to have the same distribution as the unconditional difference ĥRĥ. Now it is often already known that the two components e

′ 1 (ε ε ε ε ε ε ε ε ε ε, h 0 ) and e ′ 2 (m, ε ε ε ε ε ε ε ε ε ε, h 0 ) are
asymptotically distributed as two independent, centered, normal variables of the same order (identical to that of G ′ (h 0 )). We shall show that moreover, in the second order kernel setting, their asymptotic variances are numerically close for typical weight functions u, independently of m. So that lower bounds on the underestimation factor κ, as the one mentioned above, will be obtained (Section 3).

A second set of heuristics is the following. In order to automatically drop this underestimation factor, the expressions (1.8)-(1.9) suggest that RCL(h) should ideally have a second additional error ±e

2 (m, ε ε ε ε ε ε ε ε ε ε * , h).
If m were known, it would suffice to add such a term to RCL(h). We propose here to add -e 2 ( mg

, ε ε ε ε ε ε ε ε ε ε * , h)
where g, a pilot smoothing parameter, has to be chosen. We choose to use the sign minus and the same

ε ε ε ε ε ε ε ε ε ε *
as the one already used for generating RCL because the whole additional error now "exactly mimics" up to the factor -1 the intrinsic error of RCL, and because the corresponding new randomized criterion, that we call the "augmented randomized" version of CL (or ARCL ), then has the natural expression

ARCL(h) := RCL(h) -e 2 ( mg , ε ε ε ε ε ε ε ε ε ε * , h) = n -1 (I -A h )y 2 U + 2n -1 ε ε ε ε ε ε ε ε ε ε * , A h (A g y + ε ε ε ε ε ε ε ε ε ε * ) -A g y U .
(1.12)

where the adjustment of the residual can be recognized similar to "a bootstrap realization of the optimism" in the work of [START_REF] Efron | How biased is the apparent error rate of a prediction rule?[END_REF]. There, the average of such bootstrap realizations was used to estimate an "expected optimism" whose definition was not limited to linear A h (•) (nor to quadratic ∆) in case of which it is a function of m. As far as we know, there is no published paper touching on the methodology proposed here for kernel or spline-like smoothers where each bootstrap realization is used, via numerical minimization of ARCL, to produce a randomized choice.

The corresponding penalized criterion AR G(h) naturally uses the following augmented-randomizedtrace function in place of t(h) in the definition (1.6) -(1.7) :

AR t(h) := ε ε ε ε ε ε ε ε ε ε * , A h (A g y + ε ε ε ε ε ε ε ε ε ε * ) -A g y U ε ε ε ε ε ε ε ε ε ε * , ε ε ε ε ε ε ε ε ε ε * U . (1.13) Note that contrary to R t(h), this is no longer invariant by a renormalisation of ε ε ε ε ε ε ε ε ε ε *
and an estimate of σ is thus required for generating AR G(h). In practice many good estimates exist for estimating σ. We fell that the arguments which follow should still hold with σ 2 replaced by any one of these variance estimates.

We consider in this paper, the natural automatic choice g = ĥGCV for the pilot bandwidth. Other choices might be "better" but this one is very appealing since it does not add any complications in practice.

Now, with enough regularity assumptions and n large enough, the following stochastic representation typically approximates well the true behavior of the minimizer of ARCL, say ĥAR ,

-M ′′ (h 0 ) ĥAR -ĥ ≈ ARCL ′ -CL ′ (h 0 ) ≈ -e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 ) -e ′ 2 (A h0 y, ε ε ε ε ε ε ε ε ε ε * , h 0 ) ≈ -e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 ) -e ′ 2 (m + A h0 ε ε ε ε ε ε ε ε ε ε, ε ε ε ε ε ε ε ε ε ε * , h 0 ) -e ′ 2 (A h0 m -m, ε ε ε ε ε ε ε ε ε ε * , h 0 ).
(1.14)

The third term is typically negligible relatively to the second. Let us define

v 1 :=8n -2 σ 4 tr((A ′ h0 A h0 ) T U 2 A ′ h0 A h0 ), v 2 :=8n -2 σ 4 tr(A ′ h0 T U 2 A ′ h0 ) b 2 :=4n -2 σ 2 A ′ h0 m 2 U 2
(1.15) (these are analogs of V 1 , V 2 and B 2 of Theorem 2.1; note that Nychka (1991, Section 3), in addition of the bootstrap strategy, also studies an analytic strategy which uses similar formulae for general spline smoothers and replaces the unknowns terms by consistent estimates). Straightforward differentiation of e 1 and e 2 and variance calculus show, assuming for example that ε

ε ε ε ε ε ε ε ε ε and ε ε ε ε ε ε ε ε ε ε * are Gaussian, that, conditionally to
a given y, the first term and the second one are independently distributed with an approximate law for the first typically given by N (0, v 2 ) and an exact law for the second equal to N (0, 4n

-2 σ 2 A ′ h0 m+A ′ h0 A h0 ε ε ε ε ε ε ε ε ε ε 2 U 2
) which, for large n, will be closed to N (0, b 2 + (1/2)v 1 ) with a high probability; and thus the conditional distribution of the sum of the three terms will typically be an inflated version of the "target" N (0, b 2 + v 2 ), i.e. the approximate unconditionned law of -M ′′ (h 0 )( ĥ -ĥ0 ), with inflation (or overestimation) factor given by

κ = (b 2 + v 2 + (1/2)v 1 ) / (b 2 + v 2 ) ≤ 1 + (1/2)(v 1 /v 2 ).
We will see (Sections 5-6) how to easily eliminate this biais. However it is worth to point out that κ will be not much larger than 1 in many context. For example, for spline-like smoother where A h and A ′ h are symmetric matrices with a common base of eigenvectors, and where the eigenvalues of A h are bounded by 1, we have, for example for

U = I, that v 1 /v 2 = tr(A ′ h0 A h0 ) 2 /tr(A ′ h0 )
2 ≤ 1 and thus 1 ≤ κ ≤ 3/2. All this will be rigorously proved in the kernel context, Section 4, and extended to other contexts Remark 8.2, for large sample. Moreover we will prove in Section 5 (and 6) that we can in fact linearly combine the variance estimate obtained from the simulated ĥR 's and the one from the simulated ĥAR 's to obtain an asymptotically unbiased estimate of var ĥ -ĥ0 .

1.4. Comparison with a more classical bootstrap methodology [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF] has proposed a parametric bootstrap approach to approximate the distribution of ĥ -ĥ0 . His approach seems quite natural from a practical point of view: after having estimated m and σ using cross-validation, one can generate new independent data sets by taking these estimates in place of the true m and σ. Since typically these estimates are consistent as n → ∞, one may think that these data sets will have a distribution similar as that of the observed data set for large n. One can then compute the empirical distribution of the differences between cross-validation choice and optimal choice associated to these simulated data sets and use this distribution as an approximation of the desired distribution of ĥCV -ĥ0 . Simulations with some typical test functions m in [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF] demonstrate that this method may be rather effective. However, note that the theoretical results in [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF] state that the mean of the bootstrap distribution actually converges to a constant (< 1) multiply of h 0 . Asymptotic theory suggests that this may be corrected by using an oversmoothed mg for the generation of the pseudo-data; but in practice choosing such a pilot bandwidth g is a difficult matter. This (uncorrected) bootstrap methodology is compared to our proposal in the simulation study, Section 8. The randomized-choices methodology will prove to be more accurate in these experiments.

Statement of results for kernel regression estimates

For the sake of clarity and simplicity, we shall state our results for the case of one-dimensional data, in the widely studied setting of second order kernel as in Section 1.1. However, as in [START_REF] Girard | Asymptotic comparison of (partial) cross-validation, GCV and randomized GCV in nonparametric regression[END_REF] we relax the condition of equispaced design: we assume that there exists a smooth distribution function F over [0, 1] such that x i := F -1 ((i -0.5)/n), i = 1, • • • , n, and, for simplicity, that the density f = F ′ is known.

Then we can still use a simple explicit modified version of (1.2): mh (x) := 1 nhf (x)

K(

x-xj h )y j , for which an asymptotic study remains relatively easy. We shall require the following classical assumptions. a) The errors ε i are iid with mean 0 and all other moments finite. b) K is symmetric, compactly supported and has a Hölder continuous second derivative. c)

m is C 2 [0, 1]. d) f is C 2 [0, 1] and f (x) ≥ c > 0 on the support of u which is assumed C 1 [0, 1].
As is usual in asymptotic studies of kernel estimates, the minimization of the various selection criteria is assumed to be restricted to an interval H n = [n -1+ǫ , n -ǫ ] for an arbitrary (small) constant ǫ > 0. Let us define

J u (m) = ((mf ) ′′ ) 2 f -1 u u.
Then it is known (e.g. Härdle and Marron (1985)), that

h 0 ∼ C 0 n -1/5 with C 0 = (C 1 /C 2 ) 1/5 and C 1 = σ 2 u K 2 , C 2 = x 2 K 2 J u (m) u, (2.1)
and that ĥ h0 → 1 and ĥ0 h0 → 1 in probability. And setting C 3 :=5C 1 /C 3 0 we have

M (h 0 ) ∼ 5 4 C 1 C 0 n -4/5 , M ′′ (h 0 ) ∼ C 3 n -2/5 .
Well known important works by [START_REF] Rice | Bandwidth choice for nonparametric regression[END_REF] and HHM have established the following asymptotic stochastic behaviors of the non-randomized selectors (the extension to nonequidistant designs is analyzed in Girard 1998) that we recall for the sake of completeness. To simplify the statements of properties that are shared by various selectors, in the following, ĥ (respectively ĥR ) also denotes the minimizer of CV(h) or the one of CL(h) (resp. the minimizer of RCL(h)) in addition to the G-selectors (resp. the R G-selectors).

Theorem 2.1. (HHM, Girard (1998)) Under a), b), c) and d) (i.e. the assumptions of HHM except that the design may be non-equidistant, its density f being C 2 and bounded from below on the support of u),

we have

L C 3 n 3/10 (h 0 -ĥ0 ) → N (0, B 2 + V 1 ), L C 3 n 3/10 ( ĥ -ĥ0 ) → N (0, B 2 + V 2 )
with, denoting by L the kernel L(x) = -xK ′ (x),

B 2 = 4C 2 0 σ 2 x 2 K 2 J u 2 (m) u 2 , V 1 = 8 C 3 0 σ 4 (K * K -K * L) 2 u 2 , V 2 = 8 C 3 0 σ 4 L 2 u 2 .
In the following, the convergence of the conditional distribution of a randomized quantity X n given y (i.e. the randomization-based distribution of X n ) toward X will be stated in probability, in the sense that for any fixed t where the distribution of X is continuous, |P (X n ≤ t|y) -P (X ≤ t) | → 0 in probability.

Theorem 2.2. Under the assumptions of Theorem 2.1, assumptions on the generated w identical to those

on σ -1 ε ε ε ε ε ε ε ε ε ε and assuming w independent from ε ε ε ε ε ε ε ε ε ε, L C 3 n 3/10 ( ĥR -ĥ) y → N (0, V 2 )
in probability, where the constants C 3 and V 2 are the same ones as in Theorem 2.1.

Proof: Let X n denotes C 3 n 3/10 ( ĥRĥ). The unconditional version of this result, that is, the convergence in distribution of X n toward N (0, V 2 ), can be classically shown from the approximation

X n = n 7/10 e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 ) + o P (1), (2.2)
where e 1 is defined in (1.18), which is immediately obtained by combining the two linearizations stated in Girard (1998, equations (2.5) and (3.3)), and from the convergence in distribution (Lemma 4 of the Appendix of Girard (1998)):

n 7/10 e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 ) → N (0, V 2 ). (2.3)
Recall that one classical way to prove the well known fact that (2.2) and (2.3) are sufficient for X n → N (0, V 2 ), is to obtain the following inequality which holds for any η > 0 (e.g. section 20.6 of Cramer 1970):

|P (X n ≤ t) -F (t)| ≤ max{|F n (t + η) -F (t)|, |F n (t -η) -F (t)|} + 2P (|r n (h 0 , ε ε ε ε ε ε ε ε ε ε, ε ε ε ε ε ε ε ε ε ε * )| > η) ,
where

F n is the distribution of n 7/10 e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 ), F is the distribution N (0, V 2 ), and r n (h 0 , ε ε ε ε ε ε ε ε ε ε, ε ε ε ε ε ε ε ε ε ε * ) = X n - n 7/10 e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 )
; next, to observe that this bound can be made arbitrarily small with an appropriate choice of η and with n large enough.

Concerning the conditional X n , by noticing that F n is also the distribution of

n 7/10 e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 )|y, it is
easy to see that this inequality actually also holds with the probability P (•) replaced by the conditional probability P (•|y). Now, it suffices to observe that the so-obtained bound of |P (X n ≤ t|y) -F (t)| can also be made arbitrarily small in probability since, of course, α n (y

) = P (|r n (h 0 , ε ε ε ε ε ε ε ε ε ε, ε ε ε ε ε ε ε ε ε ε * )| > η|y) → P 0 (which results from E|α n (y)| = P (|r n (h 0 , ε ε ε ε ε ε ε ε ε ε, ε ε ε ε ε ε ε ε ε ε * )| > η) → 0 and Markov inequality).
Recalling that all the selectors are restricted to a compact set, one can define SD RAND by SD RAND := var 1 2 ĥR y .

(2.4)

We could have used E ( ĥRĥ) y 2 instead of (2.4). For large samples, this would be equivalent.

However, for finite samples this is, of course, always greater than var

1 2
ĥR y , and might then improve it in view of the underestimator factor we discuss in the following. We come back to this point in the simulation study.

By simply combining Theorem 2.1 and Theorem 2.2 and using the expression of C 5 0 given by (2.1), we now obtain the following corollary which describes the behavior of the variances of the normal asymptotic laws (that we simply call asymptotic variances):

Corollary 2.3. Under the assumptions of Theorem 2.2, SD RAND , defined in (2.4), is a consistent estimate of κ var 1 2 ĥ -ĥ0 , in the sense that ĥR y has the same asymptotic variance than κ( ĥ -ĥ0 ), where

κ = V 2 B 2 + V 2 1/2 = 1 2 J u 2 (m) J u (m) K 2 L 2 + 1 -1/2
≤ 1.

We will study in more depth this downward bias of SD RAND in the next section.

3. The underestimation factor κ is never "too small" for second order kernel

At first look, one may think that the underestimation factor κ should be estimated in order to appropriately inflate SD RAND . However, we show in this section that for second order kernel and typical u, V 2 can never be a small part of the whole variance B 2 + V 2 and thus κ is reasonably close to 1. From the expression of κ in Corollary 2.3, we see that an important case can be first considered: the case u 2

proportional to u since κ is then independent of m and of the density f . We thus give in Table 3.1 the values of this factor κ in this case, for four typical kernels K, with a now standard terminology (e.g. [START_REF] Priestley | Non-parametric function fitting[END_REF]). Notice that the first and last columns are extremal cases (see Lemma 3.1 below). It is seen that these values of κ are all reasonably close to 1.

our precise objective here. And this becomes even more acute for kernel of higher order, as detailed in Remark 8.2. So we consider the augmented-randomized-version of CL (i.e. ARCL defined in (1.12)) as well as AR G(h) using the augmented-randomized-trace function defined in (1.13).

To simplify the presentation, we assume σ is known in (1.12)-(1.13). However many estimates with good asymptotic property exist for estimating σ and the following result should remain true when replacing σ by one of them, for example the estimate of [START_REF] Rice | Bandwidth choice for nonparametric regression[END_REF].

As discussed in Section 1.3, we use here the natural automatic choice g = ĥGCV for the pilot bandwidth.

Note that the condition on g in Theorem 4.1 below is then satisfied.

Assuming also the standard conditions of Theorem 2.2, it can be shown that the minimizer of ARCL(h)

and those of any one criterion of the family AR G(h) are all equivalent up to second order. So, let ĥAR denote a generic one of these selectors.

We claim that the simulation of ĥAR given y, provides useful inferences on the underlying distribution of ĥ -ĥ0 (by comparing Theorem 4.1 to Theorem 2.1). It is important to note that we do not require any smoothness condition on the underlying m(•). Note that our approach is quite different from the natural bootstrap approach where generated pseudo-data

y * = mg + ε ε ε ε ε ε ε ε ε ε * would produce "bootstrap
replications" of both ĥ and ĥ0 , as considered by [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF].

Theorem 4.1. Under the assumptions of Theorem 2.2, and assuming that the pilot bandwidth g (random or not) used in ARCL or in AR G satisfies, for some ǫ > 0, |g -

h 0 |/h 0 = O p (n -ǫ ), then L C 3 n 3/10 ( ĥAR -ĥ) y → N (0, V 2 + B 2 + 1 2 V 1 ) in probability,
where the constant C 3 , B, V 1 and V 2 are the same ones as in Theorem 2.1.

Proof: The first step is to show that there still exists a approximate linearized form for expressing ĥAR .

We only sketch its derivation because it uses similar lines of proof as for ĥ in HHM or ĥR in [START_REF] Girard | Asymptotic comparison of (partial) cross-validation, GCV and randomized GCV in nonparametric regression[END_REF].

For example, letting B h denote the smoothing operator associated with the kernel L, what is required in supplement of (A.8) of HHM is now sup

h,g |h-h0|+|g-h0|≤n -1/5-ǫ | 1 nh (A h -B h )A g y, ε ε ε ε ε ε ε ε ε ε * U - 1 nh 0 (A h0 -B h0 )A h0 y, ε ε ε ε ε ε ε ε ε ε * U | = o P (n -7/10 )
which can be shown by similar bounds and partitioning argument to those used in HHM. This permits us to replace both ĥAR and g by h 0 in the right hand term of the approximation ĥARh 0 ≈ -ARCL ′ -M ′ ( ĥAR )/M ′′ (h 0 ) obtained by standard first order Taylor expansion and then, subtracting the known analog Taylor approximation of ĥh 0 , to derive (with e 1 and e 2 defined as in (1.8))

-M ′′ (h 0 ) ĥAR -ĥ = -e ′ 1 (ε ε ε ε ε ε ε ε ε ε * , h 0 ) -e ′ 2 (A h0 y, ε ε ε ε ε ε ε ε ε ε * , h 0 ) + o P (n -7/10 ).
The second step is to establish the limiting normal distribution of e

′ 2 (A h0 y, ε ε ε ε ε ε ε ε ε ε * , h 0 ) conditioned to y. First, we can imitate the proof for L e ′ 2 (m, ε ε ε ε ε ε ε ε ε ε, h)/var 1/2 e ′ 2 (m, ε ε ε ε ε ε ε ε ε ε, h) → N (0, 1
) which uses the Lindeberg condition (e.g. [START_REF] Eubank | Confidence regions in Non-parametric Regression[END_REF]. For this, it suffices to observe that e

′ 2 (A h0 y, ε ε ε ε ε ε ε ε ε ε * , h 0 ) = 2n -1 h 0 -1 (B h0 - A h0 )A h0 y, ε ε ε ε ε ε ε ε ε ε * U where the vector (B h0 -A h0 ))A h0 y has ith element h 0 2 1 2 x 2 (K -L) ( mh0 f ) ′′ (x i )f -1 (x i ) + o P (h 0 2 ),
uniformly in i, and quadratic mean

n -1 ||(B h0 -A h0 ))A h0 y|| 2 U 2 = h 4 0 1 2 x 2 (K -L) 2 J u 2 ( mh0 ) u 2 (1 + o P (1)).
Second, the additional term

1 2 V 1 results from n -1 ||(B h0 -A h0 )A h0 y|| 2 U 2 ∼ n -1 ||(B h0 -A h0 )A h0 m|| 2 U 2 + σ 2 n -1 h 0 -1 (K * K -K * L) 2 u 2
, in probability, where the two terms are of the same order and can be checked to be proportional to B 2 and 1 2 V 1 (for this, note x 2 (L -K) = 2 x 2 K).

Let us now define SD AUG-RAND by SD AUG-RAND := var 1 2 ĥAR y , (4.1) then Corollary 4.2. Under the assumptions of Theorem 4.1, SD AUG-RAND , defined in (4.1), is a consistent estimate of κ var 1/2 ( ĥ -ĥ0 ), in the sense that ĥAR y has the same asymptotic variance than κ( ĥ -ĥ0 ),

where

κ 2 = B 2 + V 2 + 1 2 V 1 B 2 + V 2 = 1 + 1 2 V 1 V 2 κ 2 ≥ 1.
It is important to point out that the overestimation factor κ is also typically reasonably close to 1. For example, for any positive K, it can be shown that V 1 ≤ V 2 (e.g. HHM) and thus κ ≤ 3/2 = 1.2247. If K is the Gaussian kernel, we obtain the numerical value V 1 /V 2 = √ 2/8 and thus κ ≤ 1 + √ 2/16 = 1.0432.

Such an appealing closeness to 1 will be also obtained for the cubic smoothing spline setting, Section 6.

Estimating var( ĥ -ĥ0 ) without asymptotic bias

In the general case where u 2 ∝ u, neither κ nor κ are known, but they are known functions of the functionals J u (m) and J u 2 (m). So, one might be tempted to work with standard estimates of these functionals. In fact, we do not need such estimates to derive a consistent estimate of var( ĥ -ĥ0 ). Indeed, Theorem 5.1. Under the assumptions of Theorem 2.1, and for κ and κ given in Corollaries 2.3 and 4.2, let us assume that n 3/10 SD RAND /κ 2 and n 3/10 SD AUG-RAND /κ 2 both converge in probability toward the asymptotic variance of n 3/10 ( ĥ -ĥ0 ); then

n 3/10 SD AUG-RAND 2 - 1 2 (K * K -K * L) 2 L 2 n 3/10 SD RAND 2 (5.1)
converges in probability toward the asymptotic variance of n 3/10 ( ĥ -ĥ0 ).

Proof: It can be derived from Corollary 4.2 that

κ 2 /κ 2 = 1/κ 2 + (1/2)(V 1 /V 2 )
, where V1 V2 is a constant depending only on K. On the other hand, SD AUG-RAND / SD RAND is a consistent estimate of κ κ . Thus a consistent estimate of κ (or of J u 2 (m)/J u (m)) is deduced of this ratio. Plugging it in n 3/10 SD RAND /κ 2 , it is checked this is equivalent to use the stated linear combination of SD Note that the convergence in probability of n 3/10 SD RAND or n 3/10 SD AUG-RAND is not claimed in Corollaries 2.3 or 4.2; indeed it is well known that convergence in distribution of a random sequence does not imply convergence in moment. A classical technical condition sufficient for this, is a uniform integrability condition that we were not able to obtain here. Nevertheless Theorem 5.1 is already useful for percentile-based variance estimates (such as those in (7.3) below) in place of empirical moments: such estimates are classical in case of asymptotic normality and are even often proposed as robust standard deviation estimates (e.g. [START_REF] Efron | The Bootstrap, the Jacknife and Other Resampling Plans[END_REF]); their consistency holds true in our setting (see Section 7.2).

Extension to other nonparametric regression settings

From the heuristics of Section 1.3 (mathematical rigor as for the theorems above, is clearly possible in some other settings as outlined in Remark 8.2) these accuracy estimators have a large potential but the estimator given Theorem 5.1 is, at first glance, not as easily applicable as SD RAND or SD AUG-RAND since the appropriate general factor -(1/2)(v 1 /v 2 ) in this linear combination (where v 1 and v 2 are defined in (1.5)) is specific to the setting. However, note that for all the settings where v 1 /v 2 is independent of h 0 , a general version of this estimate could be constructed simply be replacing this constant by an approximation that can be obtained, for a given family A h of smoothing operators, by a natural simulation based approach as follows: choose one (or several) realistic "true" function m's and a corresponding σ, simulate several data-sets for each m, and adjust the above linear combination so as to optimize its observed average performance (in a least-squares sense for example) on these simulated problems. This should work in the general smoothing spline setting where A h = (I + hΩ) -1 with Ω a given symmetric ≥ 0 matrix, at least for u ≡ 1. Indeed it can be checked, by using that A ′ h = -h(I -A h )A h in the formulae given Section 1.3, that the appropriate constant in this linear combination becomes

-(1/2)(v 1 /v 2 ) = -(1/2)trA 4 h (I -A h ) 2 /trA 2 h (I -A h ) 2 (6.1)
where, strictly speaking, h should be set to the unknown h 0 ; but it is known in the spline literature that such a ratio of traces is typically asymptotic to a constant independent of h (e.g. [START_REF] Nychka | The average posterior variance of a smoothing spline and a consistent estimate of the average squared error[END_REF] for a proof of this for the one-dimensional polynomial smoothing spline). In fact, v 1 /v 2 can be simply evaluated form the expression (6.1) where the replacement of h 0 by ĥGCV should not degrade the approximation.

Moreover numerical values close to a few percents can often be expected for this ratio, in case of which the correction of SD AUG-RAND is not really mandatory. Indeed for example in the cubic spline setting one obtains (e.g. [START_REF] Kou | On the efficiency of selection criteria in spline regression[END_REF]

) (1/2)trA 4 h (I -A h ) 2 /trA 2 h (I -A h ) 2 ∼ (1/2) [Γ(15/4)/Γ(6)] / [Γ(7/4)/Γ(4)] = 0.118.

Various uses of these accuracy estimates

7.1. The empirical distribution of the randomized GCV choices, and the one for the augmented-randomized version, both obtained by simulation, can be used in several different ways. The first one is classical and was already developed in Nichka (1991) : even if, strictly speaking, such a simulation furnishes a consistent, say 95%, confidence interval for the difference ĥ -ĥ0 , it is natural to very simply transform it (by translating it by ĥ) in a prediction (called "confidence" in Nichka (1991)) interval for ĥ0 which can be interpreted as follows. For the data in hand, the ∆-optimal ĥ0 can be claimed to belong to the produced interval with a confidence level of 95%, that is, the claim being incorrect for only 5% (assuming that the large sample approximation is correct) of the possible replicated data sets from the stochastic structure of the problem at hand. If a data analyst draws the curve estimates for every bandwidth in this interval and observes "little" changes at each design point of the essential support of the weight function u, then he/she may be assured, with a confidence level of 95%, that he/she has found the ∆-optimal fit (up to the mentioned "little changes" in it). The epithet "little" for a change in the curve at a design point may be made more precise by adding the phrase "relatively to the width of a confidence band for the curve around this point"

(we come back to the topic of confidence bands in Section 7.3) since perturbing the center of a such band by "little" changes will then have "little" consequence on its coverage probability.

As advocated by [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF], a good practice is to report the 2 curve-or surface-estimates which correspond to the lower and upper points of the produced interval. When, contrarily to the first situation, two fits associated with the endpoints of such a (1α)100% confidence interval have rather different aspects as α passes from, say, 30% to 5%, the choice of α has a large influence on the report for the data analyst. Note that, concerning two different bandwidths which are in a such 95% confidence interval for ĥ0 , one can say that there is no objective reasons to think that one of the two corresponding fits should be better (for the ∆ criterion) than the other. In the following we state a further interpretation of the two fits reported, which may give a guide for the choice of α in terms of expected loss.

7.2. Consider, again to simplify the presentation, the setting of Sections 2-5. Let τ a given constant and let Ĉ any statistic (i.e. a function of y) assumed to be a consistent estimate of τ (B 2 + V 2 ) 1/2 /C 3 ; and consider ĥGCV (τ ) = ĥGCV +n -3/10 Ĉ a perturbed version of ĥGCV . Then it can be shown that C 3 n 3/10 ( ĥGCV (τ )-ĥ0 ) is asymptotically distributed as (B 2 + V 2 ) 1/2 N (τ, 1) and that the known asymptotic laws for the "excess loss" (i.e. the difference between the (data-driven) resulting ∆ and the best possible ∆)

L n ∆( ĥ) -∆( ĥ0 ) → 1 2C 3 (B 2 + V 2 )χ 2 1 (7.1)
which are stated in HHM and [START_REF] Girard | Asymptotic comparison of (partial) cross-validation, GCV and randomized GCV in nonparametric regression[END_REF], can be easily generalized for ĥGCV (τ ) in place of ĥ : the standard χ 2 1 must simply be replaced by decentred chi-square [N (τ, 1)] 2 (the proofs are very similar to the ones in HHM, Girard (1998)) and are not repeated here).

On the other hand, for α, β ∈ (0, 0.5) let us denote by ĥR (β) and ĥR (α) the (β)100th and (1α)100th

percentiles of the randomization-based distribution, i.e. defined by

P ( ĥR ≤ ĥR (β) y) = β, P ( ĥR ≤ ĥR (α) y) = 1 -α (7.2)
Besides providing a standard simulation-based (1 -(α + β))100% confidence interval, these percentiles also furnish variance estimates whose consistency can be proven by classical techniques, namely

n 3/10 ĥR (α) -ĥR (β) z 1-α -z β , n 3/10 ĥR (α) - ĥ z 1-α or n 3/10 ĥ -ĥR (β) -z β (7.3)
where z β and z 1-α are the (β)100th and (1-α)100th percentiles of the N (0, 1), both converge in probability toward κ(B 2 +V 2 ) 1/2 /C 3 . By combining the above statement, using z 1-α times the second of this standard deviation estimates in place of Ĉ in the above perturbation, which implies the coincidence ĥ + n -3/10 Ĉ = ĥR (α), we obtain Theorem 7.1. Let α ∈ (0, 0.5). Under the assumptions of Theorem 2.2,

L n ∆ ĥR (α) -∆( ĥ0 ) → 1 2C 3 (B 2 + V 2 ) N (z 1-α κ, 1)
2 where the constants C 3 , B, V 2 and κ are the same ones as in Theorem 2.1 and Corollary 2.3.

By combining this with the expression (7.1) for the asymptotic law of ∆( ĥ) -∆( ĥ0 ), one obtain a very simple bound about what is sacrificed, comparatively to the CV choice, in the fit corresponding to the upper point ĥR (α), in terms of relative increase of the excess loss, in average, also called relative "risk regret" :

E ∆ ĥR (α) -∆( ĥ0 ) E ∆( ĥ) -∆( ĥ0 ) = 1 + (z 1-α κ) 2 ≤ 1 + z 2 1-α (7.4)
where E denotes here the asymptotic variance.

For example, a data analyst concerned by data reduction will look after a kind of "parsimonious yet near optimal" fit of his data : Theorem 7.1 implies that the fit corresponding to the upper point ĥR (0.16) has a theoretical justification "almost" as good as the CV fit, namely the relative risk regret is bounded by 1 + z 2 0.84 = 2. As another example, by choosing ĥR (0.32), the bound becomes 1 + z 2 0.68 = 1.21. Recall that modifying the CV choice toward more smoothing is commonly made in practice (e.g. Gu (2002) Section 6.3.2): the standard modification consists of multipying the trace term in CL(h) (or t(h) in the GCV-like criteria) by a factor larger than 1 (e.g. 1.4). This is had hoc but often works well. It would certainly be interesting to compare these two ways of producing "more parsimonious yet near optimal" fits.

It is obvious that an analog of Theorem 7.1 exists for a lower point : the ∆ performance of the fit obtained with ĥR (β) in place of ĥR (α) has the same expression with N (z 1-α κ, 1) replaced by N (-z β κ, 1).

Let us mention how these simple risk-regret estimators might be used in the production of inferences about m. Of course for any selector h(y) the risk E∆(h(y)) measures the (local) performance of mh(y) ; one of the possible uses of an estimate of this risk is the following: A methodology which is now very popular, in the smoothing spline setting, for accompanying the point estimate mĥ (x), with an interval estimate, is the so-called 90% Bayesian confidence interval mĥ (x) ± z 0.05 σ trU A ĥ/trU (7.6) which is well known, at least for the one-dimensional polynomial spline and u ≡ 1, to have good frequentist "across the curve" properties for m deterministic provided ĥ is a good estimate of ĥ0 or of h 0 . Such an interval indeed takes account the fact that mĥ is biaised. However a possibly very serious weakness is that the error in ĥ is not taken in account. Indeed the true (1γ)100% interval for any h(y) should be mh(y) (x) ± z γ/2 E∆(h(y)/trU (7.7)

while the justification of (7.6) stems on a consistency result σ 2 trU A ĥ ≈ E∆( ĥ) which is in fact a consistency toward E∆(h 0 ) and thus neglect the risk-regret term(s). Before to see how this could be attacked, we point out that the theoretical justification of (7.7) developed by [START_REF] Nychka | The average posterior variance of a smoothing spline and a consistent estimate of the average squared error[END_REF] can be easily seen to also holds for a weight function u ≡ 1 and more general smoothing spline or kernel estimates: the only required assumptions are that the associated smoothing operator is symmetric and that the chosen weight function

u is smooth enough so that A h u ≈ u (where u = (u(x 1 , • • • , u(x n )) T ) over the domain of considered h's.
The meaning of "across the curve" is now that the coverage probability of the interval (7.7) is not only relatively to the distribution of ε ε ε ε ε ε ε ε ε ε but also for x drawn with the discrete probability which assigns the mass u(x i )/trU at each x i .

Now by decomposing ∆(h(y) = ∆( ĥ0 ) + ∆(h(y) -∆( ĥ0 ) we see that 1. The risk regret from using h(y), i.e. the expectation of the second term, is a minorant of E∆(h(y). Of course for any one of the selectors h(y) studied theoretically here, the relative size of this minorant, compared to E∆(h(y), tends to 0 as n → ∞, but we must keep in mind that it may take very large sample sizes before this risk regret is a small part of E∆(h(y); indeed for the second order kernel setting of Section 2, this relative size decreases to zero at a rate only of n -1 /n -4/5 = n -1/5 (and n -1/(2k+1) for kernel of order k, see Remark 8.2); in fact this becomes n -1/10 if we come back (by taking square roots)

to the correct scale for the width of the confidence band. Moreover it can be checked that, for any fixed n, one can obtain a setting for which the constant factor c, in the approximation cn -1/10 of the relative size of this minorant, is arbitrarily large simply by taking a more concentrated weight function u. Thus a first use of risk-regret estimates like (7.5) may be as follows, using the fact that, for the trivial no-smoothing curve estimator (or "h(y) = 0"), the value of the risk, E∆(h(y), is known (and equal to σ 2 trU ) : indeed as soon as the data analyst observes that the estimated lower bound (7.5) for E∆( ĥCV ) is not a small part of σ 2 trU , the CV-choice should be considered as likely not better than the raw data in terms of ∆-performance. (this may be called a "first diagnostic for a failure of (local)

CV"). A similar assessment can clearly also be constructed for the toward-more-parcimony selectors ĥR (α) (or ĥAR (α)) proposed section 7.2. For example let us state this for ĥAR (α): a consistent estimate of the associated risk regret E ∆( ĥAR (α) -∆( ĥ0 )

∼ (2nC 3 ) -1 [B 2 + V 2 + z 2 1-α (B 2 + V 2 + (1/2)V 1 ] is RR AUG-RAND - 1 2 V 1 V 2 RR RAND + z 2 1-α RR AUG-RAND .
2. By decomposing further ∆( ĥ0 ) = ∆(h 0 ) -∆(h 0 ) -∆( ĥ0 ) and using the previous asymptotic laws, we see that, for example, again for the choice h(y):= ĥAR (α),

E ∆( ĥAR (α)) ≈ M (h 0 ) - 1 2nC 3 (B 2 + V 1 ) + 1 2nC 3 (B 2 + V 2 )(1 + (z 1-α κ) 2 ),
the error in this approximation being o(n -1 ) we take for granted that the invoked moments at finite size n are asymptotic to the moments of the asymptotic laws. The sum of the second and third terms is easily seen, from the above results, to be consistently (thus up to o P (n -1 ) ) estimated by

(1 - V 1 V 2 ) RR RAND + z 2 1-α RR AUG-RAND . (7.8)
An as much accurate estimate, say M0 , of the minimum M (h 0 ) is available, in theory, provided the well known asymptotic version (nh) -1 C 1 + (1/4)C 2 h 4 of M (h) is accurate enough and, as is well known too, some additional smoothness can be assumed so that √ n-consistent estimates of C 2 and h 0 are available. Unfortunately such accurate M0 requires a good choice of pilot bandwidths which is not easy. This is beyond the scope of this article, and we do not include any experimental study.

3. In some contexts, the M (h) function can be bounded, say by B(h), from a priori bound on the "roughness" of m (like on J u (m) as in [START_REF] Speckman | Spline smoothing and optimal rates of convergence in nonparametric regression models[END_REF]). There, another possible use of the above correction (7.8) of M (h 0 ) may be to apply it to the minimum of B(h) over the considered h (since

M (h 0 ) ≤ M (argmin h B(h)) ≤ min h B(h))
to yield an accurate conservative confidence band centered at mĥ AR(α) .

Further work is certainly useful to develop a methodology for constructing confidence bands or regions ready for practical use.

It is noteworthy that the produced confidence-bands for the toward-more-parcimony selectors ĥR (α)

and ĥAR (α), have a width which increases with z 2 1-α : this is very satisfactory and is in contrast with the standard procedure when one uses the mentioned standard had hoc modification of the criterion to translate ĥ toward a little bit more smoothing and still retains the above formula (7.6) for the width of a confidence-band (indeed trU A h is generally a decreasing function of h).

Simulations, discussion and further remarks

In practice, one may thus propose the following general methodology: First compute SD RAND by repeatedly rerunning the fast randomized version of, for example, GCV. Then, if this estimated lower bound (for the accuracy of ĥGCV ) is reasonably "small" (otherwise the use of GCV is very questionable), one next computes SD AUG-RAND by repeating the whole previous procedure with, this time, the augmented randomized GCV using the cross-validation choice as pilot bandwidth g and one of the usual estimates of σ: one thus obtains a less optimistic estimate of the variance of ĥGCV -ĥ0 which is thus more reliable for providing a conservative confidence interval for ĥ0 . Note that this second stage has the same cost as the first one and requires no further analytic or programming effort (if computing exact GCV is costly, g may be taken as the mean of the randomized choices of the first stage, instead of the exact ĥGCV ).

Simulation study

For our experiments, we used a setting which has already been used in [START_REF] Rice | Bandwidth choice for nonparametric regression[END_REF], HHM, and many others. So, the chosen mean function was m(x) = x 3 (1x) 3 , and we used equispaced designs for simulating data sets. We consider a gaussian noise, i.e. ε

ε ε ε ε ε ε ε ε ε ∼ N (0, σ 2 I).
A periodic version of the kernel estimate (1.2) was considered, as this is described e.g. in HHM. This is appropriate because m is "smoothly" periodic with period 1. We used u ≡ 1 as weight function, which is possible in this periodic setting. The computational cost of the simulation is then greatly reduced by using fast Fourier transforms (the randomized criteria are also easily generated in the Fourier space) as in the studies mentioned above and in [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF].

The kernel function was taken as the biweight

K(x) = (15/8)(1 -4x 2 ) 2 1 [-.5,.5] (x).
In this simulation study, among all the selectors which are asymptotically equivalent to GCV, we only consider the Mallows' C L selector for which we assume σ is known (since this criterion is then the only one which is exactly unbiased).

We report here the results that we obtained for the following two sample sizes: n = 128 and n = 512;

and two noise levels: σ = 0.0015 and σ = 0.011 (the first one corresponds to a peak signal-to-noise equal to 10), these four cases being rather representative of "small" and "large" samples with "small" and "large" noise levels, respectively.

For each of the so-obtained four cases, 100 samples (i.e. 100 data sets) of size n are simulated. A reasonably good approximation of the targeted var 1 2 ( ĥCL -ĥ0 ) is first computed by its empirical version (denoted SD in the following) over these 100 samples. For each of these 100 samples, first, we generated 200 realizations of the randomized selector RCL, each one using a w drawn from N (0, I) to compute SD RAND ; next, using the exact CL selector as pilot bandwidth, we generated 200 realizations of the augmented randomized version (ARCL), also using ε First, it should be observed that, as suggested from the theory, it is reasonable here to neglect the "bias" in quantifying the accuracy of the CL selector: more precisely, one see in the third column of these tables, that the mean of ĥCL -ĥ0 is rather "negligible", in absolute value, as compared to SD.

ε ε ε ε ε ε ε ε ε * 's from N (0, σ 2 I),
Second, one observes in Tables 8.1 and 8.2 that the "sandwich" ( SD RAND , SD AUG-RAND ), is quite satisfactory: the expected value of this interval approximately includes SD and it has reasonably small width in view of the variability of each of these two "bounds". The "noise" in these randomization estimates is of course not negligible : in each of the four cases, their standard deviation is roughly one fourth of SD.

But in view of the rather large range of the different values for the scale SD, even with such a noise, these randomization-based estimates prove to be useful in these experiments.

In this simulation study, we also considered the bootstrap approach mentioned in the Introduction, using, essentially as in [START_REF] Nychka | Choosing a range for the amount of smoothing in nonparametric regression[END_REF], the exact CL choice for each data set as a simple pilot bandwidth for generating 200 pseudo-data sets (each one drawn from N ( mĥ CL , σ 2 I)). A striking fact, summarized in Tables 8.1 and 8.2, is that SD BOOT tends to over-estimate SD, even more than SD AUG-RAND , especially for small samples. The noise in all these accuracy estimates is rather similar, except in the "small sample, small noise" case where SD BOOT is much more variable than SD RAND or SD AUG-RAND .

As another benchmark, we also consider the "plug-in" variance estimate proposal of HHM, mentioned in the Introduction, which is based on the asymptotic analysis of this specific setting: from Theorem 2.1 and the expressions of C 0 and C 3 (given in [START_REF] Girard | Asymptotic comparison of (partial) cross-validation, GCV and randomized GCV in nonparametric regression[END_REF] for the case f ≡ 1), the following expression for the targeted asymptotic variance is easily derived:

var ĥCL -ĥ0 ∼ h 3 0 u 2 K 2 u 2 4 25 J u 2 (m) J u (m) + 2 L 2 K 2 .
So, in the case u 2 ∝ u, a natural and very simple approach is to plug-in ĥCL in place of h 0 in this expression, and to compute the required constant : this yields the estimate (.5796) 2ĥ3 CL for our choice of K which is displayed in the last column of Tables 8.1 and 8.2. It is remarkable that this plug-in estimate, although it uses a detailed analysis of this specific setting, has a noise which is either only comparable to or worse than that of the randomization-based estimates. Note that, in the general case u 2 ∝ u, this plug-in proposal requires estimates of the ratio J u 2 (m)/J u (m) and an additional noise can then be apprehended.

Notice that in these experiments, we observed relatively small differences between SD RAND and SD AUG-RAND , and so computing SD AUG-RAND actually was not very useful. For the same reason, we have also not implemented here the more sophisticated estimate of Section 5. However, since this second stage has the same cost as the first one, it may well be worth implementing it in other applications.

Further remarks

Remark 8.1. Extensions (and limitations) of the proposed randomization-based methodology. Various variants and extensions are possible. For example 1. In this article, we focus on the estimation of var 1/2 ( ĥ -ĥ0 ). Other measures for the "accuracy" of ĥ could have been chosen, or a logarithmic scale (or other monotone transformations, like the one giving the popular equivalent degree of freedom, e.g. [START_REF] Hastie | A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines[END_REF]) could have been taken as well. Of course, a standard deviation is a good summary only for an error distribution which is roughly normal. Note that this is rather typically the case for kernel bandwidths, as shown by the simulations in HHM, Section 4, and e.g. in Hall and Johnstone (1992, Section 5 and Remark 7.2).

2. When considering the generality of the GCV-type criteria and their randomized versions, and the underlying heuristics (Section 1.3), one may guess this methodology can be applied to many other curve-, surface-or image-estimation techniques. However, one must keep in mind that the approximate stochastic representations of ĥ-ĥ0 and ĥRĥ of Section 1.3 are based on local quadratic approximation of the considered criteria at their minima, and thus, at least second order differentiability of A h as a function of h seems to be required as in e.g. Kneip (1995, Section 5.2). Note also that the theory here does not deal with the behavior near zero of the distribution of the GCV selector (since this selector has been appropriately restricted to the interval H n ); for a discussion of this disturbing behavior for small samples, see [START_REF] Wahba | Behavior near zero of the distribution of GCV smoothing parameter estimates[END_REF]. One may also guess that the augmented randomization technique should generally provide a certain correction of the downward bias of SD RAND ; but in situations where such a correction would be mandatory (theory says this might be the case for kernels of very high order, see Remark 8.2 below), the upper-bound property of SD AUG-RAND stems on the consistency properties of the GCV-selector. Note that, if the practitioner desires that SD AUG-RAND be reliably conservative, then a pilot parameter which would undersmooth, compared to the optimal ĥ0 , seems more appropriate than one which oversmoothes.

3. Nonparametric curve or surface estimators using a multidimensional smoothing parameter could be considered as well. One of the popular multidimensional examples is additive modeling with backfitting. There, the randomization-based methodology could naturally produce accuracy estimates, simultaneously for each component of the vector ĥ of GCV bandwidth estimates; an estimated covariance matrix (or even a confidence region for the vector of the optimal bandwidths) might be constructed. This proposal would, of course, need to be investigated carefully.

Remark 8.2. Extensions of the theory. The theoretical results above can be extended to several contexts.

  to compute SD AUG-RAND . Table8.2. σ = .011. Summary (mean and, in brackets, standard deviation) over 100 samples

	n = 128	ĥ0 .209 (.033) (.038) (SD = .054) ĥCL ĥCL -ĥ0 .200 -.009	SD RAND SD AUG-RAND SD BOOT .5796 .052 .060 .107 .053 ĥ3/2 CL (.014) (.014) (.034) (.014)
	n = 512	.156 (.025) (.026) (SD = .040) .151 -.005	.035 (.010)	.040 (.010)	.050 (.013)	.034 (.008)
	Table 8.1. σ = .0015. Summary (mean and, in brackets, standard deviation) over 100 samples
	n = 128	ĥ0 .480 (.093) (.129) (SD = .188) ĥCL ĥCL -ĥ0 .457 -.023	SD RAND SD AUG-RAND SD BOOT .5796 .168 .195 .255 .185 ĥ3/2 CL (.044) (.042) (.044) (.072)
	n = 512	.346 (.072) (.090) (SD = .133) .326 -.020	.109 (.027)	.130 (.028)	.155 (.032)	.111 (.042)

In order to derive general lower bounds for κ (like the last column of Table 3.1), let us state the two following lemmas:

Lemma 3.1. Assume that K is a continuously differentiable function with compact support. Then

Proof: By integration by parts, we have K 2 =x(2K ′ K) whose square is bounded by (2xK

by Cauchy-Scharwz.

Furthermore, it can be easily shown that this inequality is "sharp", i.e. there exists a minimizing sequence of kernels (denoted "K min " in Table 3.1) for which the ratio |xK ′ | 2 / K 2 converges to the lower bound 1/4.

Lemma 3.2. Under the assumptions of Theorem 2.1,

Proof: Note that this ratio can be written

It immediately follows that:

Theorem 3.3. If u 2 ∝ u then the underestimation factor κ which appears in Corollary 2.3 always satisfies:

Otherwise, this still holds, for any m and f , with √ 3 replaced by 1 + 2 u u 2 max u. When the chosen u is a Gaussian density of any width, u u 2 max u = √ 2 if the integrals are taken over (-∞, +∞).

The case u Gaussian considered here is of practical interest since it is a natural weight function often used to define local bandwidths (e.g. [START_REF] Vieu | Nonparametric regression: local optimal bandwidth choice[END_REF]). Note that if the width of a Gaussian u is chosen to be small enough relatively to the local variation of (mf ) ′′ f -1 , then J u 2 (m)/J u (m) will actually be closer to 1 than to the approximate upper bound √ 2.

Augmenting the randomization error in ĥR for another accuracy estimate

We have seen in the previous section, that the intrinsic error (1.9) of the randomized criterion RCL has an additional error (compared to the error (1.8) of CL) which actually may be not large enough for Using such a fit might be useful for a data analyst concerned by discovering features (like peaks) in a data set.

Analogs of these results can also easily be stated for an upper point and a lower point (denoted respectively by ĥAR (α) and ĥAR (β)) of the empirical conditional distribution of repeated bandwidth choices by the augmented-randomized-trace version. The only change is that κ is then replaced by κ.

7.3. One can in fact estimate the absolute value of the risk excess itself, simply from the randomizationbased distribution. For this task one might estimate the curvature of ∆(•) near ĥ0 , estimate E ĥ -ĥ0 2 and invoke the usual second order Taylor approximation ∆( ĥ) -∆( ĥ0 ) ≈ (1/2) ĥ -ĥ0

it can be shown, as a consequence of the previous results, that the variations of, say, GCV(•) over the population of the choices, thus conditional to the observed y, furnishes (at least) as much information :

Theorem 7.2. Let G be any one of the G-selectors defined in (1.6)-(1.7) and ĥ (resp. ĥR and ĥAR ) the minimizer of G (resp. of its randomized-trace version and its augmented-randomized-trace version). Under the assumptions of Theorem 4.1,

in probability, where the constants C 3 , B, V 2 , κ and κ are the same ones as in Theorem 2.1 and Corollaries 2.3 and 4.2.

Proof: Similar steps as in the proof of Theorems 2.2 and 4.1, along with the stochastic approximations for ∆(h) -∆( ĥ0 ) developed in HHM and [START_REF] Girard | Asymptotic comparison of (partial) cross-validation, GCV and randomized GCV in nonparametric regression[END_REF] for h = ĥ, ĥR , give the results.

Either by sample moment estimates or by appropriate percentiles-based estimates, one can then produce a "sandwich" of estimates of lower and upper bounds reasonably (from the "nearness" of κ and κ to 1) closed to the risk regret due to the inaccuracy of CV. One example of percentiles-based estimate, denoted by RR RAND for future references, is simply the median of the empirical conditional distribution of G( ĥR ) -G( ĥ) divided by the median (≈ 0.45) of χ 2 1 . RR AUG-RAND can be defined similarly with ĥR replaced by ĥAR . From Theorem 7.2, these 2 estimates can be easily proven to be consistent up to the bias factor κ 2 and κ 2 respectively. By an argument similar as in the proof of Theorem 5.1, it is seen that

is an asymptotically unbiased estimate of the asymptotic risk regret from using ĥ, E ∆( ĥ) -∆( ĥ0 ) . Note that when u 2 ∝ u, one can check that the ratio (B 2 + V 2 )/C 3 is in fact independent of m (e.g. Section 4.4.5 of Hall and Johnstone 1992 for the case f ≡ 1); but this does not hold any longer when u 2 ∝ u.

We point out that typically the underestimation factor κ still remains reasonably close to 1. 1. Consider for example the setting of Section 2, with only K replaced by a higher order kernel, that is

and assume that m and f have uniformly continuous kth derivatives. Then all the theoretical results of Sections 2, 4, 5 and 6, still hold with appropriate changes in the exponents of convergence (n 3/10 becomes 3/2(2k+1) ) and with modifications in the constants classically calculated. A very simple modification can be easily proved for the expression of κ in Corollary 2.3, which becomes:

, where again

, for which Lemmas 3.1 and 3.2 are also valid, and thus Theorem 3.3 still holds with √ 3 replaced by √ 1 + k, and the factor 2 simply replaced by k in the expression 1 + 2 u u 2 max u. This implies thus reasonable lower bounds for κ when the order k remains small. 2. A justification for this methodology is also possible in the context of smoothing spline estimate.

For one-dimensional cubic smoothing spline, a convincing development in Nychka (1991, Section 3) indicates that the asymptotic formulae agree with those of a fourth-order kernel regression estimate.

In the case of equidistant design, by exploiting that the smoothing spline framework is very well approximated by a known Gaussian white noise framework using Butterworth filter (whose range is parametered by h) in the frequency domain, to which correspond a particular kernel for which compactness condition (a) in Section 3 can be dropped, a rigorous development can be made as in Hall and Johnstone (1992, section 6.2).

3. Extensions to multidimensional settings could also be stated. See HHM and Girard (1998, Remark 3.2) for a two dimensional example using Gasser-Muller kernel estimate for which the result of this paper can be extended with classical modifications in the rate (the power 3/10 replaced everywhere by (d + 2)/(8 + 2d) with d = 2) and in the expressions for the constants; see Girard (1998, Remark 3.2) for the definition of the analog of J u (m).

Remark 8.3. It would be interesting to study the rate of convergence for these accuracy estimates: this deserves further study.

Remark 8.4. Computational cost. Since the theoretical results in this paper only describe the behavior of e.g. SD RAND which uses an infinite number of randomized choices, it is natural to wonder how many repeated minimizations one has to simulate. A complete theoretical analysis of this question would first require an answer to Remark 7.3. In our above experiments where n = 128 or n = 512, sufficiently stabilized estimates were already obtained by simulating 50 randomized choices : Simulating further 150 randomized choices, as was done here, did not actually produce great improvement.